content: package edu.stanford.rsl.tutorial.phantoms;
  
  import ij.ImageJ;
  import edu.stanford.rsl.conrad.numerics.SimpleMatrix;
  import edu.stanford.rsl.conrad.numerics.SimpleVector;
  
  public class SheppLogan extends Phantom{
  	
  	final private SimpleMatrix Ellipses; 
  	
  	
  	public SheppLogan(int xy, boolean type) {
  		super(xy,xy,"Shepp-Logan-Phantom");
  		
  		if (type)
  			Ellipses = ShepOrig();
  		else
  			Ellipses = ShepMod();
  		CreatePhantom();
  	}
  	
  	
  	public SheppLogan(int xy) {
  		super(xy,xy,"Shepp-Logan-Phantom");
  		Ellipses = ShepMod();
  		CreatePhantom();
  	}
  	
  	private void CreatePhantom()
  	{
  	
  		// Iterate over all pixels and sum up intensity values of all corresponding ellipses
  		double sizeX = (double)super.getSize()[0];
  		double sizeY = (double)super.getSize()[1];
  		
  		for (int i=0; i < super.getSize()[0]; ++i)
  		{
  			double x = ((double)i-(sizeX-1)/2.0) / ((sizeX-1)/2.0);
  			for (int j=0; j < super.getSize()[1]; ++j)
  			{
  				double y = ((double)j-(sizeY-1)/2.0) / ((sizeY-1)/2.0);
  				super.setAtIndex(i, super.getSize()[1]-j-1, 0.f);
  				for (int k=0; k < Ellipses.getRows(); ++k)
  				{
  					// Extract the ellipse properties here
  					double xc = x - Ellipses.getElement(k, 3);
  					double yc = y - Ellipses.getElement(k, 4);
  					double phi = Ellipses.getElement(k, 5)*Math.PI/180.0;
  					double cos = Math.cos(phi);
  					double sin = Math.sin(phi);
  					double asq = Ellipses.getElement(k, 1)*Ellipses.getElement(k, 1);
  					double bsq = Ellipses.getElement(k, 2)*Ellipses.getElement(k, 2);
  					double Val = Ellipses.getElement(k, 0);
  					
  					// Check if this pixel is part of the ellipse, if yes, add the given intensity value to it
  					double help = Math.pow((xc*cos + yc*sin),2.0);
  					double help2 = Math.pow((yc*cos - xc*sin),2.0);
  					if ( help/asq + help2/bsq <= 1.0 )
  						super.setAtIndex(i, super.getSize()[1]-j-1, super.getAtIndex(i, super.getSize()[1]-j-1) + (float)Val);
  				}
  				
  			}
  		}
  	}
  	
  	private SimpleMatrix ShepOrig()
  	{
  		// Original Shepp-Logan Phantom according to:
  		
  		// Shepp, L. A., & Logan, B. F. (1974). 
  		// The Fourier reconstruction of a head section-LA Shepp. 
  		// IEEE Transactions on Nuclear Science, NS-21, 21�43.
  		
  		// One row describes properties for a single ellipse
  		// Colum Values: A    a     b    x0    y0    phi
  		
  		SimpleMatrix Shep = new SimpleMatrix(10,6);
  		
  		Shep.setRowValue(0, new SimpleVector(new double[] {2.0, 0.69, 0.92, 0, 0, 0}));
  		Shep.setRowValue(1, new SimpleVector(new double[] {-0.98, 0.6624, 0.8740, 0, -0.0184, 0}));
  		Shep.setRowValue(2, new SimpleVector(new double[] {-0.02, 0.1100, 0.3100, 0.22, 0.0, -18.0}));
  		Shep.setRowValue(3, new SimpleVector(new double[] {-0.02, 0.1600, 0.4100, -0.22, 0.0, 18.0}));
  		Shep.setRowValue(4, new SimpleVector(new double[] {0.01, 0.2100, 0.2500, 0, 0.35, 0}));
  		Shep.setRowValue(5, new SimpleVector(new double[] {0.01, 0.0460, 0.0460, 0, 0.1, 0}));
  		Shep.setRowValue(6, new SimpleVector(new double[] {0.01, 0.0460, 0.0460, 0, -0.1, 0}));
  		Shep.setRowValue(7, new SimpleVector(new double[] {0.01, 0.0460, 0.0230, -0.08, -0.605, 0}));
  		Shep.setRowValue(8, new SimpleVector(new double[] {0.01, 0.0230, 0.0230, 0, -0.606, 0}));
  		Shep.setRowValue(9, new SimpleVector(new double[] {0.01, 0.0230, 0.0460,  0.06, -0.605, 0}));
  		return Shep;
  	}
  	
  	private SimpleMatrix ShepMod()
  	{
  	// Modified (better contrast) Shepp-Logan Phantom according 
  	// to P. A. Toft, "The Radon Transform, Theory and Implementation" 
  	// (unpublished dissertation), p. 199.
  	
  	// One row describes properties for a single ellipse
  	// Colum Values: A    a     b    x0    y0    phi
  	
  	SimpleMatrix Shep = new SimpleMatrix(10,6);
  	
  	Shep.setRowValue(0, new SimpleVector(new double[] {1.0, 0.69, 0.92, 0, 0, 0}));
  	Shep.setRowValue(1, new SimpleVector(new double[] {-0.8, 0.6624, 0.8740, 0, -0.0184, 0}));
  	Shep.setRowValue(2, new SimpleVector(new double[] {-0.2, 0.1100, 0.3100, 0.22, 0.0, -18.0}));
  	Shep.setRowValue(3, new SimpleVector(new double[] {-0.2, 0.1600, 0.4100, -0.22, 0.0, 18.0}));
  	Shep.setRowValue(4, new SimpleVector(new double[] {0.1, 0.2100, 0.2500, 0, 0.35, 0}));
  	Shep.setRowValue(5, new SimpleVector(new double[] {0.1, 0.0460, 0.0460, 0, 0.1, 0}));
  	Shep.setRowValue(6, new SimpleVector(new double[] {0.1, 0.0460, 0.0460, 0, -0.1, 0}));
  	Shep.setRowValue(7, new SimpleVector(new double[] {0.1, 0.0460, 0.0230, -0.08, -0.605, 0}));
  	Shep.setRowValue(8, new SimpleVector(new double[] {0.1, 0.0230, 0.0230, 0, -0.606, 0}));
  	Shep.setRowValue(9, new SimpleVector(new double[] {0.1, 0.0230, 0.0460,  0.06, -0.605, 0}));
  	return Shep;
  	}
  	
  	
  	public static void main (String [] args){
  		new ImageJ();
  		SheppLogan test = new SheppLogan(256);
  		test.show("Shepp Logan");
  		
  		SheppLogan test2 = new SheppLogan(256,true);
  		test2.show("Shepp Logan");
  		
  		}
  
  }
  /*
   * Copyright (C) 2010-2014 Andreas Maier
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
  */