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). */