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