content: package edu.stanford.rsl.tutorial.atract;
  
  import ij.ImageJ;
  import edu.stanford.rsl.conrad.data.numeric.Grid2D;
  import edu.stanford.rsl.conrad.data.numeric.Grid3D;
  import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators;
  import edu.stanford.rsl.conrad.geometry.trajectories.Trajectory;
  import edu.stanford.rsl.conrad.phantom.NumericalSheppLogan3D;
  import edu.stanford.rsl.conrad.utils.Configuration;
  import edu.stanford.rsl.tutorial.cone.ConeBeamBackprojector;
  import edu.stanford.rsl.tutorial.cone.ConeBeamCosineFilter;
  import edu.stanford.rsl.tutorial.cone.ConeBeamProjector;
  import edu.stanford.rsl.tutorial.fan.CosineFilter;
  import edu.stanford.rsl.tutorial.filters.RamLakKernel;
  import edu.stanford.rsl.tutorial.motion.estimation.CylinderVolumeMask;
  import edu.stanford.rsl.tutorial.parallel.ParallelBackprojector2D;
  import edu.stanford.rsl.tutorial.parallel.ParallelProjector2D;
  import edu.stanford.rsl.tutorial.phantoms.DotsGrid2D;
  import edu.stanford.rsl.tutorial.phantoms.MickeyMouseGrid2D;
  import edu.stanford.rsl.tutorial.phantoms.Phantom;
  import edu.stanford.rsl.tutorial.phantoms.Phantom3D;
  import edu.stanford.rsl.tutorial.phantoms.SimpleCubes3D;
  import edu.stanford.rsl.tutorial.phantoms.Sphere3D;
  import edu.stanford.rsl.tutorial.phantoms.TestObject1;
  import edu.stanford.rsl.tutorial.phantoms.UniformCircleGrid2D;
  
  /**
   * This class is an example implementation to show the usage of the atract filter.
   * The 1D atract implementation is run on a 2D image.
   * The 2D atract implementation is run on a 2D image.
   * the 2D atract implementation is run on a 3D volume.
   * 
   * @author Marco Boegel (Reco Project 2012 - Individual Project)
   *
   */
  public class DisplayAtract {
  
  	/**
  	 * @param args
  	 */
  	public static void main(String[] args) {
  		// sinogram params
  		double maxTheta=Math.PI,       // [rad]
  			   deltaTheta=Math.PI/360, // [rad]
  			   maxS=150,               // [mm]
  			   deltaS=1.0;             // [mm]
  		// image params
  		int imgSzXMM = 150,            // [mm]
  			imgSzYMM = imgSzXMM;       // [mm]
  		float pxSzXMM = 1.0f,          // [mm]
  			pxSzYMM = pxSzXMM;         // [mm]
  
  		//float focalLength = 400, maxBeta = (float) Math.PI*2, deltaBeta = maxBeta / 200, maxT = 200, deltaT = 1;
  
  		int phantomType = 3; // 0 = circle, 1 = MickeyMouse, 2 = TestObject1,
  							// 3=DotsGrid
  		// size in grid units
  		int imgSzXGU = (int) Math.floor(imgSzXMM / pxSzXMM), // [GU]
  			imgSzYGU = (int) Math.floor(imgSzYMM / pxSzYMM); // [GU]
  		new ImageJ();
  
  		ParallelProjector2D projector = new ParallelProjector2D(maxTheta, deltaTheta, maxS, deltaS);
  
  		// image object
  		Phantom phantom;
  		switch (phantomType) {
  		case 0:
  			phantom = new UniformCircleGrid2D(imgSzXGU, imgSzYGU);
  			break;
  		case 1:
  			phantom = new MickeyMouseGrid2D(imgSzXGU, imgSzYGU);
  			break;
  		case 2:
  			phantom = new TestObject1(imgSzXGU, imgSzYGU);
  			break;
  		case 3:
  			phantom = new DotsGrid2D(imgSzXGU, imgSzYGU);
  			break;
  		default:
  			phantom = new UniformCircleGrid2D(imgSzXGU, imgSzYGU);
  			break;
  		}
  
  		phantom.setSpacing(pxSzXMM, pxSzYMM);
  		// origin is given in (negative) world coordinates
  		phantom.setOrigin(-(imgSzXGU * phantom.getSpacing()[0]) / 2.0,
  				-(imgSzYGU * phantom.getSpacing()[1]) / 2.0);
  		phantom.show();
  		Grid2D grid = phantom;
  		
  		// create projections
  		Grid2D sinoRay = projector.projectRayDriven(grid);
  		sinoRay.show("Sino CPU Ray");
  		
  		int maxSIndex = (int) (maxS / deltaS + 1);
  		int maxThetaIndex = (int) (maxTheta / deltaTheta + 1);
  		Kollimator koll = new Kollimator(maxThetaIndex, maxSIndex);
  		
  		koll.applyToGrid(sinoRay, 50);
  		sinoRay.show("Kolli");
  	
  		AtractFilter1D at = new AtractFilter1D();
  		at.applyToGrid(sinoRay);
  		
  		ParallelBackprojector2D bp = new ParallelBackprojector2D(imgSzXMM, imgSzYMM, pxSzXMM, pxSzYMM);
  		Grid2D recon = bp.backprojectRayDriven(sinoRay);
  
  		recon.show("Reconstruction");
  
  		
  		Grid2D sinoRay2 = projector.projectRayDriven(grid);
  		sinoRay2.show("Sino CPU Ray");
  		RamLakKernel ramLak = new RamLakKernel((int) (maxS / deltaS),
  				deltaS);
  		for (int theta = 0; theta < sinoRay2.getSize()[0]; ++theta) {
  			 ramLak.applyToGrid(sinoRay2.getSubGrid(theta));
  		}
  		Grid2D recon2 = bp.backprojectRayDriven(sinoRay2);
  		recon2.show("Normal Reconstruction");
  		
  		Grid2D sinoRay3 = projector.projectRayDriven(grid);
  		sinoRay3.show("Sino CPU Ray");
  		
  		
  		koll.applyToGrid(sinoRay3, 50);
  		sinoRay3.show("Kolli");
  	
  		AtractFilter2D at2 = new AtractFilter2D();
  		at2.applyToGrid2D(sinoRay3);
  		
  		Grid2D recon3 = bp.backprojectRayDriven(sinoRay3);
  
  		recon3.show("Reconstruction");
  		
  		
  		
  		
  		@SuppressWarnings("unused")
  		double focalLength = 800, maxT = maxS, deltaT = deltaS, gammaM = Math
  				.atan2(maxT / 2.f /*- 0.5*/, focalLength), maxBeta = Math.PI * 2, deltaBeta = maxBeta / 180;
  
  Configuration.loadConfiguration();
  		
  		Configuration conf = Configuration.getGlobalConfiguration();
  		
  		Trajectory geo = conf.getGeometry();
  		int maxU = geo.getDetectorWidth();
  		int maxV = geo.getDetectorHeight();
  		double deltaU = geo.getPixelDimensionX();
  		double deltaV = geo.getPixelDimensionY();
  		int imgSizeX = geo.getReconDimensionX();
  		int imgSizeY = geo.getReconDimensionY();
  		int imgSizeZ = geo.getReconDimensionZ();
  		Phantom3D test3D = new SimpleCubes3D(imgSizeX, imgSizeY, imgSizeZ);
  
  		Grid3D grid3 = test3D;
  		ConeBeamProjector cbp = new ConeBeamProjector();
  		grid3.show("object");
  		Grid3D sino = cbp.projectRayDrivenCL(grid3);
  
  		sino.show("sinoCL");
  		ConeBeamCosineFilter cbFilter = new ConeBeamCosineFilter(conf
  				.getGeometry().getSourceToDetectorDistance(), conf
  				.getGeometry().getDetectorWidth(), conf.getGeometry().getDetectorHeight(),1.0, 1.0);
  
  		for (int i = 0; i < conf.getGeometry().getProjectionStackSize(); ++i) {
  			cbFilter.applyToGrid(sino.getSubGrid(i));
  
  			float D = (float) conf
  					.getGeometry().getSourceToDetectorDistance();
  			NumericPointwiseOperators.multiplyBy(sino.getSubGrid(i), (float) (D*D * Math.PI / geo.getNumProjectionMatrices()));
  		}
  		sino.show("sinoFilt");
  		
  	
  		
  		Kollimator koll3 = new Kollimator();
  		
  		koll3.applyToGrid(sino, 100, 100);
  		sino.show("Kolli");
  		
  		AtractFilter2D af = new AtractFilter2D();
  		af.applyToGrid2D(sino);
  		
  		
  		
  		
  		ConeBeamBackprojector cbbp = new ConeBeamBackprojector();
  		Grid3D recImage = cbbp.backprojectPixelDrivenCL(sino);
  		CylinderVolumeMask mask = new CylinderVolumeMask(imgSizeX, imgSizeY, imgSizeX/2, imgSizeY/2, 20);
  		mask.applyToGrid(recImage);
  		recImage.show("recImage");
  
  	
  		
  		Grid3D sino2 = cbp.projectRayDrivenCL(grid3);
  
  		
  		koll3.applyToGrid(sino2, 100, 100);
  		sino2.show("Kolli");
  		
  		RamLakKernel ramK = new RamLakKernel(maxU, deltaU);
  		for (int i = 0; i < conf.getGeometry().getProjectionStackSize(); ++i) {
  			cbFilter.applyToGrid(sino2.getSubGrid(i));
  			for (int j = 0;j <maxV; ++j)
  				ramK.applyToGrid(sino2.getSubGrid(i).getSubGrid(j));
  			float D = (float) conf
  					.getGeometry().getSourceToDetectorDistance();
  			NumericPointwiseOperators.multiplyBy(sino2.getSubGrid(i), (float) (D*D * Math.PI / geo.getNumProjectionMatrices()));
  		}
  		sino2.show("sinoFilt");
  		
  	
  		Grid3D recImage2 = cbbp.backprojectPixelDrivenCL(sino2);
  		
  		mask.applyToGrid(recImage2);
  		recImage2.show("recImage");
  
  
  	
  	}
  
  }
  /*
   * Copyright (C) 2010-2014  Marco B�gel
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
  */