package edu.stanford.rsl.tutorial;
import ij.ImageJ;
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.FanBeamBackprojector2D;
import edu.stanford.rsl.tutorial.fan.FanBeamProjector2D;
import edu.stanford.rsl.tutorial.filters.RamLakKernel;
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.Sphere3D;
import edu.stanford.rsl.tutorial.phantoms.TestObject1;
import edu.stanford.rsl.tutorial.phantoms.UniformCircleGrid2D;
/**
* Simple example that computes and displays a reconstruction.
*
* @author Recopra Seminar Summer 2012
*
*/
public class DisplayReconstruction {
/**
* @param args
*/
public static void main(String[] args) {
// sinogram params
double maxS = 200, // [mm]
deltaS = 1.0, // [mm]
maxTheta = Math.PI * 2, // [rad]
deltaTheta = maxTheta / 200; // [rad]
// image params
int imgSzXMM = 100, // [mm]
imgSzYMM = imgSzXMM; // [mm]
float pxSzXMM = 1.0f, // [mm]
pxSzYMM = pxSzXMM; // [mm]
// fan beam bp parameters
@SuppressWarnings("unused")
double focalLength = 600, maxT = maxS, deltaT = deltaS, gammaM = Math
.atan2(maxT / 2.f /*- 0.5*/, focalLength), maxBeta = Math.PI * 2, deltaBeta = maxBeta / 180;
int phantomType = 1; // 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();
@SuppressWarnings("unused")
ParallelProjector2D projector = new ParallelProjector2D(maxTheta,
deltaTheta, maxS, deltaS);
@SuppressWarnings("unused")
ParallelBackprojector2D backprojector = new ParallelBackprojector2D(
imgSzXGU, imgSzYGU, pxSzXMM, pxSzYMM);
@SuppressWarnings("unused")
FanBeamProjector2D fanBeamProjector = new FanBeamProjector2D(
focalLength, maxBeta, deltaBeta, maxT, deltaT);
@SuppressWarnings("unused")
FanBeamBackprojector2D fanBeamBackprojector = new FanBeamBackprojector2D(
focalLength, deltaT, deltaBeta, imgSzXGU, imgSzYGU);
// image object
// MickeyMouseGrid2D phantom = new MickeyMouseGrid2D(imgSzXGU,
// imgSzYGU);
// DotsGrid2D grid = new DotsGrid2D(imgSzXGU, imgSzYGU);
// UniformCircleGrid2D phantom = new UniformCircleGrid2D(imgSzXGU,
// imgSzYGU);
@SuppressWarnings("unused")
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;
}
Configuration.loadConfiguration();
Configuration conf = Configuration.getGlobalConfiguration();
Trajectory geo = conf.getGeometry();
focalLength = geo.getSourceToDetectorDistance();
//int maxU = geo.getDetectorWidth();
//int maxV = geo.getDetectorHeight();
int maxU_PX = geo.getDetectorWidth();
int maxV_PX = geo.getDetectorHeight();
double deltaU = geo.getPixelDimensionX();
double deltaV = geo.getPixelDimensionY();
double maxU = (maxU_PX) * deltaU;
double maxV = (maxV_PX) * deltaV;
int imgSizeX = geo.getReconDimensionX();
int imgSizeY = geo.getReconDimensionY();
int imgSizeZ = geo.getReconDimensionZ();
Phantom3D test3D = new Sphere3D(imgSizeX, imgSizeY, imgSizeZ);
// test3D.show("Sphere");
@SuppressWarnings("unused")
NumericalSheppLogan3D shepp3d = new NumericalSheppLogan3D(imgSizeX,
imgSizeY, imgSizeZ);
// 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();
Grid3D grid = test3D;
// Grid3D grid =shepp3d.getNumericalSheppLoganPhantom();
grid.show("object");
ConeBeamProjector cbp = new ConeBeamProjector();
//Grid3D sino = cbp.projectPixelDriven(grid);
Grid3D sino = cbp.projectRayDrivenCL(grid);
sino.show("sinoCL");
ConeBeamCosineFilter cbFilter = new ConeBeamCosineFilter(focalLength, maxU, maxV, deltaU, deltaV);
RamLakKernel ramK = new RamLakKernel(maxU_PX, deltaU);
for (int i = 0; i < conf.getGeometry().getProjectionStackSize(); ++i) {
cbFilter.applyToGrid(sino.getSubGrid(i));
//ramp
for (int j = 0;j <maxV_PX; ++j)
ramK.applyToGrid(sino.getSubGrid(i).getSubGrid(j));
float D = (float) conf
.getGeometry().getSourceToDetectorDistance();
NumericPointwiseOperators.multiplyBy(sino.getSubGrid(i), (float) (D*D * Math.PI / geo.getNumProjectionMatrices()));
}
sino.show("sinoFilt");
ConeBeamBackprojector cbbp = new ConeBeamBackprojector();
Grid3D recImage = cbbp.backprojectPixelDrivenCL(sino);
recImage.show("recImage");
if (true)
return;
// create projections
/*
* Grid2D sinoRay = projector.projectRayDriven(grid); GridImage
* sinoImage = new GridImage(sinoRay); sinoImage.show("Sinogram - Ray");
*/
// Grid2D sinoPx = projector.projectRayDriven(grid);
// sinoPx.show("P Sinogram");
// Grid2D fanBeamSinoRay = fanBeamProjector.projectRayDrivenCL(grid);
// fanBeamSinoRay.show("FB CL Sinogram");
// // if (true)
// // return;
// RamLakKernel ramLak = new RamLakKernel((int) (maxS / deltaS),
// deltaS);
//ConeBeamCosineFilter cKern = new ConeBeamCosineFilter(focalLength,
// maxU, maxV, deltaU, deltaV);
// // ParkerWeights pWeights = new ParkerWeights(focalLength, maxT,
// deltaT,
// // maxBeta, deltaBeta);
// // pWeights.applyToGrid(fanBeamSinoRay);
// // pWeights.show();
//for (int theta = 0; theta < sino.getSize()[0]; ++theta) {
//cKern.applyToGrid(sino.getSubGrid(theta));
// // ramLak.applyToGrid(fanBeamSinoRay.getSubGrid(theta));
// // // ramLak.applyToGrid(sinoPx.getSubGrid(theta));
//}
//
// FanBeamBackprojector2D fbp = new FanBeamBackprojector2D(focalLength,
// deltaT, deltaBeta, imgSzXMM, imgSzYMM);
// System.out.print("Pixel:");
// long t0 = System.currentTimeMillis();
// Grid2D CLPixel = fbp.backprojectPixelDrivenCL(fanBeamSinoRay);
// System.out.println(System.currentTimeMillis() - t0);
// t0 = System.currentTimeMillis();
// Grid2D CLRay = fbp.backprojectRayDrivenCL(fanBeamSinoRay);
// //fbpCL.show("FB P/CPU Reconstruction");
// System.out.println("Ray: " + (System.currentTimeMillis() - t0));
// CLRay.show("Ray");
// CLPixel.show("Pixel");
}
}
/*
* Copyright (C) 2010-2014 Andreas Maier
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/