content: package edu.stanford.rsl.tutorial.fan;
  
  import ij.ImageJ;
  import ij.io.Opener;
  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.utils.Configuration;
  import edu.stanford.rsl.conrad.utils.ImageUtil;
  import edu.stanford.rsl.tutorial.fan.FanBeamBackprojector2D;
  import edu.stanford.rsl.tutorial.fan.redundancy.ParkerWeights;
  import edu.stanford.rsl.tutorial.filters.RamLakKernel;
  
  /**
   * This class is used to compare the tutorial fan beam reconstruction pipeline 
   * with the pipeline that is implemented in CONRAD. Therefore, the forward projected 
   * phantoms are imported from a saved file and the central line of each projection
   * is extracted to build a fan beam sinogram.
   * 
   * The tutorial code and the CONRAD code are then compared after several steps, e.g.
   * Cosine weighting, ramp filtering and backprojection.
   * 
   * 
   * Note, that we do not compare any forward projection of the tutorial code and CONRAD!
   * 
   * @author berger
   * 
   */
  public class FanBeamComparisonToCONRAD {
  
  
  	/**
  	 * @param args
  	 */
  	public static void main(String[] args) {
  
  		
  		// get geometry details from global configuration ("Conrad.xml")
  		Configuration.loadConfiguration();
  		Configuration config = Configuration.getGlobalConfiguration();
  		Trajectory traj = config.getGeometry();
  		// image params
  		int imgSzXGU = traj.getReconDimensionX(), // [GU]
  				imgSzYGU = traj.getReconDimensionY(); // [GU]
  		double pxSzX = traj.getReconVoxelSizes()[0], // [mm]
  				pxSzY = traj.getReconVoxelSizes()[1]; // [mm]
  
  		// fan beam bp parameters
  		double 	maxT = traj.getDetectorWidth(), 
  				deltaT = traj.getPixelDimensionX(), 
  				focalLength = traj.getSourceToDetectorDistance(),
  				maxBeta = traj.getNumProjectionMatrices()*traj.getAverageAngularIncrement(), 
  				deltaBeta = traj.getAverageAngularIncrement(),
  				gammaM = Math.atan((maxT/2.0)/focalLength);
  
  		System.out.println(gammaM*180/Math.PI);
  
  		new ImageJ();
  		
  		// Load projection data and extract sinogram
  		Grid2D Sinogram = extractSinogram("C:\\Users\\berger\\Desktop\\WaterCylProjDatShort.zip", traj, false, "ProjectionData");
  
  		// show sinogram
  		Sinogram.clone().show("Sinogram");
  
  		// define cosine weights and RamLak kernel and Parker weights
  		RamLakKernel ramLak = new RamLakKernel((int) maxT, deltaT);
  		CosineFilter cKern = new CosineFilter(focalLength, maxT*deltaT, deltaT);
  		ParkerWeights pWeights = new ParkerWeights(focalLength, maxT, deltaT, maxBeta*Math.PI/180, deltaBeta*Math.PI/180);
  		// Apply cosine filtering
  		for (int theta = 0; theta < Sinogram.getSize()[1]; ++theta) {
  			//cKern.applyToGrid(Sinogram.getSubGrid(theta));
  		}
  		//Sinogram.clone().show("After Cosine Weighting");
  		
  		// Load cosine weighted projection data and extract sinogram
  		Grid2D SinogramCosine = extractSinogram("C:\\Users\\berger\\Desktop\\WaterCylProjDatCosineShort.zip", traj, false, "CosineData");
  		
  		Grid2D diffCosineFanCONRAD = (Grid2D)NumericPointwiseOperators.subtractedBy(Sinogram,SinogramCosine);
  		NumericPointwiseOperators.abs(diffCosineFanCONRAD);
  		//diffCosineFanCONRAD.show("Difference Cosine Weighting Tutorial vs. CONRAD");
  		
  		// Apply Parker Weights
  		pWeights.applyToGrid(Sinogram);
  		Sinogram.clone().show("After Parker Weighting");
  		
  		Grid2D SinogramParker = extractSinogram("C:\\Users\\berger\\Desktop\\oneMaskParker.zip", traj, false, "ParkerData");
  		SinogramParker.show("CONRAD after Parker");
  		
  		Grid2D diffParkerFanCONRAD = (Grid2D)NumericPointwiseOperators.subtractedBy(Sinogram,SinogramParker);
  		NumericPointwiseOperators.abs(diffParkerFanCONRAD);
  		diffParkerFanCONRAD.show("Difference Parker Weighting Tutorial vs. CONRAD");
  		
  		
  		
  		// Apply RamLak filter
  		for (int theta = 0; theta < Sinogram.getSize()[1]; ++theta) {
  			ramLak.applyToGrid(Sinogram.getSubGrid(theta));
  		}
  		Sinogram.clone().show("After RamLak Filtering Tutorial");
  
  		// Load cosine weighted projection data and extract sinogram
  		Grid2D SinogramRamLak = extractSinogram("C:\\Users\\berger\\Desktop\\WaterCylProjDatRamLakShort.zip", traj, false, "RamLakData");
  		SinogramRamLak.clone().show("After RamLak Filtering Conrad");
  		
  		Grid2D diffRamLakFanCONRAD = (Grid2D)NumericPointwiseOperators.subtractedBy(Sinogram,SinogramRamLak);
  		NumericPointwiseOperators.abs(diffRamLakFanCONRAD);
  		diffRamLakFanCONRAD.show("Difference after RamLak Tutorial vs. CONRAD");
  		
  
  		// Do the backprojection
  		FanBeamBackprojector2D fbp = new FanBeamBackprojector2D(focalLength,
  				deltaT, deltaBeta, imgSzXGU, imgSzYGU);
  		Grid2D FanbeamRecon = fbp.backprojectPixelDrivenCL(Sinogram);
  		FanbeamRecon.show("Reconstruction Result");
  		
  		
  		Opener op = new Opener();		
  		Grid3D CONRADRecon = ImageUtil.wrapImagePlus(op.openZip("C:\\Users\\berger\\Desktop\\WaterCylProjDatReconShort.zip"));		
  		CONRADRecon.show("CONRAD Reconstruction Result");
  		
  		Grid2D diffReconFanCONRAD = (Grid2D)NumericPointwiseOperators.subtractedBy(FanbeamRecon, CONRADRecon);
  		diffReconFanCONRAD.show("Difference after Recon Tutorial vs. CONRAD");
  
  	}
  	
  	
  	
  	static Grid2D extractSinogram(String filename, Trajectory traj, boolean showData, String name)
  	{
  		Opener op = new Opener();
  		Grid3D projData = ImageUtil.wrapImagePlus(op.openZip(filename));
  		if (showData)
  			projData.show(name);
  		
  		// extract sinogram out of cone beam projection data (central detector line in v direction)
  		Grid2D Sinogram = new Grid2D(projData.getSize()[0], projData.getSize()[2]);
  		for (int i=0; i < projData.getSize()[2]; ++i)
  		{
  			NumericPointwiseOperators.copy(Sinogram.getSubGrid(i), 
  					projData.getSubGrid(i).getSubGrid(traj.getDetectorHeight()/2));
  		}
  		Sinogram.setSpacing(traj.getPixelDimensionX(), traj.getAverageAngularIncrement());
  		
  		return Sinogram;
  		
  	}
  }
  /*
   * Copyright (C) 2010-2014 Martin Berger
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
  */