package edu.stanford.rsl.tutorial.atract;
import ij.ImageJ;
import ij.ImagePlus;
import ij.io.Opener;
import edu.stanford.rsl.conrad.data.numeric.Grid2D;
import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators;
import edu.stanford.rsl.conrad.utils.ImageUtil;
import edu.stanford.rsl.tutorial.fan.CosineFilter;
import edu.stanford.rsl.tutorial.fan.FanBeamBackprojector2D;
import edu.stanford.rsl.tutorial.fan.FanBeamProjector2D;
import edu.stanford.rsl.tutorial.fan.dynamicCollimation.copyRedundantData;
import edu.stanford.rsl.tutorial.fan.redundancy.BinaryWeights;
import edu.stanford.rsl.tutorial.fan.redundancy.ParkerWeights;
import edu.stanford.rsl.tutorial.filters.RamLakKernel;
import edu.stanford.rsl.tutorial.phantoms.DotsGrid2D;
import edu.stanford.rsl.tutorial.phantoms.FilePhantom;
import edu.stanford.rsl.tutorial.phantoms.MTFphantom;
import edu.stanford.rsl.tutorial.phantoms.MickeyMouseGrid2D;
import edu.stanford.rsl.tutorial.phantoms.Phantom;
import edu.stanford.rsl.tutorial.phantoms.SheppLogan;
import edu.stanford.rsl.tutorial.phantoms.TestObject1;
import edu.stanford.rsl.tutorial.phantoms.UniformCircleGrid2D;
public class AtractWithBinaryWeight {
/**
* @param args
*/
public static void main(String[] args) {
double speedup = 1;
// image params
int imgSzXMM = (int)(1024/speedup), // [mm]
imgSzYMM = imgSzXMM; // [mm]
float pxSzXMM = 1.0f, // [mm]
pxSzYMM = pxSzXMM; // [mm]
// sinogram params
@SuppressWarnings("unused")
double focalLength = 4416/speedup,
maxT = 1472/speedup, // Detector length according to image dimensions to avoid truncation
deltaT = 1.0,
//gammaM = Math.atan((maxT / 2.f - 0.5) / focalLength),
maxBeta = 200*Math.PI/180,
//deltaBeta = maxBeta / 165,
deltaBeta = maxBeta / 165,
fanAngle = Math.atan((maxT/2.0-0.5*deltaT)/focalLength);
/*double fanAngle = 10*Math.PI/180,//Math.atan((maxT/2.0-0.5*deltaT)/focalLength);
maxT = 1472/speedup, // Detector length according to image dimensions to avoid truncation
deltaT = 1.0,
focalLength = (maxT-0.5f*deltaT)*0.5f/Math.tan(fanAngle),//4416/speedup,
//gammaM = Math.atan((maxT / 2.f - 0.5) / focalLength),
maxBeta = 200*Math.PI/180,
deltaBeta = maxBeta / 165;*/
int phantomType = 4; // 0 = circle, 1 = MickeyMouse, 2 = TestObject1,
// 3=DotsGrid, 4=SheppLogan, 5=MTF, 6=File
// size in grid units
int imgSzXGU = (int) Math.floor(imgSzXMM / pxSzXMM), // [GU]
imgSzYGU = (int) Math.floor(imgSzYMM / pxSzYMM); // [GU]
new ImageJ();
// image object
Phantom phantom;
Grid2D Sinogram = null;
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;
case 4:
phantom = new SheppLogan(imgSzXGU);
break;
case 5:
phantom = new MTFphantom(imgSzXGU, imgSzYGU, 0.95, 0.7895, 1.f, 4.f);
break;
case 6:
//phantom = new FilePhantom(imgSzXGU, imgSzYGU, "D:\\!Stanford\\RecoSVN\\berger\\Data\\PigSlice.zip");
//phantom = new FilePhantom(imgSzXGU, imgSzYGU, "C:\\Users\\Martin\\Desktop\\lena.zip");
phantom = new FilePhantom(imgSzXGU, imgSzYGU, "C:\\Users\\Martin\\Desktop\\donald.zip");
break;
case 7:
// Load your own saved image here!!
Opener op = new Opener();
ImagePlus ipl = op.openZip("D:\\!Stanford\\RecoSVN\\berger\\MICCAI\\Results\\MTF\\MTFPhantom_1024_1024_r2_384.zip");
phantom = new DotsGrid2D(ipl.getWidth(),ipl.getHeight());
imgSzXMM = phantom.getWidth();
imgSzYMM = phantom.getHeight();
imgSzXGU = (int) Math.floor(imgSzXMM / pxSzXMM); // [GU]
imgSzYGU = (int) Math.floor(imgSzYMM / pxSzYMM); // [GU]
Sinogram = ImageUtil.wrapImagePlusSlice(op.openZip("D:\\!Stanford\\RecoSVN\\berger\\MICCAI\\Results\\MTF\\Sinogram_1024_1024.zip"),0, false);
Sinogram.setSpacing(deltaBeta, deltaT);
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("Phantom");
Grid2D grid = phantom;
if (phantomType != 7)
{
// Fan Beam Projection
FanBeamProjector2D fanBeamProjector = new FanBeamProjector2D(focalLength, maxBeta, deltaBeta, maxT, deltaT);
Sinogram = fanBeamProjector.projectRayDrivenCL(grid);
}
Sinogram.show("Sinogram");
// ****************** From here on we have a valid sinogram ******************************
// Binary Weights
BinaryWeights binWeights = new BinaryWeights(focalLength, maxT, deltaT, maxBeta, deltaBeta);
//binWeights.show("Binary Weights");
// Parker Weights
ParkerWeights parkWeights = new ParkerWeights(focalLength, maxT, deltaT, maxBeta, deltaBeta);
//parkWeights.show("Parker Weights");
// apply parker weights
Grid2D parkerSinogram = new Grid2D(Sinogram);
NumericPointwiseOperators.multiplyBy(parkerSinogram, parkWeights);
//parkerSinogram.show("Parker Weights Sinogram");
// apply binary Weights (artificial dynamic collimator)
Grid2D binSinogram = new Grid2D(Sinogram);
NumericPointwiseOperators.multiplyBy(binSinogram, binWeights);
//binSinogram.show("Binary Weights Sinogram");
// try to recover missing (collimated) data by known redundant part and use parker weighting
copyRedundantData copyClass = new copyRedundantData(focalLength, maxT, deltaT, maxBeta, deltaBeta);
Grid2D copiedFullSinogram = new Grid2D(binSinogram);
copyClass.applyToGrid(copiedFullSinogram);
NumericPointwiseOperators.multiplyBy(copiedFullSinogram, parkWeights);
// In this case we apply two weighting schemes: binary weights and parker weights!
// therefore the scale correction was done twice, so we need to divide once by the same scaling factor
NumericPointwiseOperators.multiplyBy(copiedFullSinogram, (float)( (Math.PI) / maxBeta));
//copiedFullSinogram.show("Copied Full Sinogram");
// show difference between parker weighting on full data and parker on copied redundancies
Grid2D nonColl_park_Sino_minus_Coll_copied_park_Sino = (Grid2D)NumericPointwiseOperators.subtractedBy(parkerSinogram, copiedFullSinogram);
nonColl_park_Sino_minus_Coll_copied_park_Sino.show("NonColl_Park_Sino minus Coll_Copied_Park_Sino");
// Cosine Filter
CosineFilter cosFilt = new CosineFilter(focalLength, maxT, deltaT);
// Filtering
RamLakKernel ramLak = new RamLakKernel((int) (maxT / deltaT), deltaT);
// Create ATRACT filter
AtractFilter1D atractFilter = new AtractFilter1D();
// Filtering for Binary ATRACT
Grid2D atractSinogram = new Grid2D(binSinogram);
Grid2D parkerAtractSinogram = new Grid2D(binSinogram);
//Grid2D atract2DSinogram = new Grid2D(binSinogram);
for (int theta = 0; theta < Sinogram.getSize()[0]; ++theta) {
cosFilt.applyToGrid(atractSinogram.getSubGrid(theta));
cosFilt.applyToGrid(parkerAtractSinogram.getSubGrid(theta));
//cosFilt.applyToGrid(atract2DSinogram.getSubGrid(theta));
}
// apply the Atract kernel
atractFilter.applyToGrid(atractSinogram, binWeights.getBinaryMask());
boolean checkMask = false;
if (checkMask) {
Grid2D mask = new Grid2D(
new float[parkerAtractSinogram.getSize()[0]*parkerAtractSinogram.getSize()[1]],
parkerAtractSinogram.getSize()[0],
parkerAtractSinogram.getSize()[1]
);
NumericPointwiseOperators.fill(mask, 1.0f);
atractFilter.applyToGrid(parkerAtractSinogram, mask);
}
/*
// apply 2D atract
LaplaceKernel2D laplace2D = new LaplaceKernel2D();
laplace2D.applyToGrid(atract2DSinogram);
AtractKernel2D atract2D = new AtractKernel2D(atract2DSinogram.getSize()[0], atract2DSinogram.getSize()[1]);
atract2D.applyToGrid(atract2DSinogram);
*/
// Filtering for Binary and Parker FBP
for (int theta = 0; theta < Sinogram.getSize()[0]; ++theta) {
cosFilt.applyToGrid(binSinogram.getSubGrid(theta));
cosFilt.applyToGrid(parkerSinogram.getSubGrid(theta));
cosFilt.applyToGrid(copiedFullSinogram.getSubGrid(theta));
ramLak.applyToGrid(binSinogram.getSubGrid(theta));
ramLak.applyToGrid(parkerSinogram.getSubGrid(theta));
ramLak.applyToGrid(copiedFullSinogram.getSubGrid(theta));
}
/*
FanBeamBackprojector2D fbp_beads = new FanBeamBackprojector2D(focalLength, deltaT, deltaBeta, 64, 64);
Grid2D[] beadSinogram = {parkerSinogram, copiedFullSinogram, atractSinogram};
String[] sinoNames = {"Parker_RamLak_CentralBead","Binary_Copied_Parker_RamLak_CentralBead","Binary_ATRACT_CentralBead"};
for (int i=0; i<beadSinogram.length; ++i)
{
Grid2D parkerResBeads = fbp_beads.backprojectPixelDriven_HighResRegion(beadSinogram[i], 1/16.f, new double[] {0, 0} );
parkerResBeads.show(sinoNames[i]);
}
*/
// Do the backprojections
FanBeamBackprojector2D fbp = new FanBeamBackprojector2D(focalLength, deltaT, deltaBeta, imgSzXMM, imgSzYMM);
Grid2D parkerRes = fbp.backprojectPixelDrivenCL(parkerSinogram);
parkerRes.show("Parker_RamLak");
Grid2D binaryRes = fbp.backprojectPixelDrivenCL(binSinogram);
binaryRes.show("Binary_RamLak");
Grid2D copiedParkerRes = fbp.backprojectPixelDrivenCL(copiedFullSinogram);
copiedParkerRes.show("Binary_Copied_Parker_RamLak");
Grid2D binaryAtractRes = fbp.backprojectPixelDrivenCL(atractSinogram);
binaryAtractRes.show("Binary_ATRACT");
//Grid2D Atract2DRes = fbp.backprojectPixelDrivenCL(atract2DSinogram);
//Atract2DRes.show("2D_ATRACT");
if (checkMask) {Grid2D parkerAtractRes = fbp.backprojectPixelDrivenCL(parkerAtractSinogram);
parkerAtractRes.show("Parker_ATRACT");
}
}
}
/*
* Copyright (C) 2010-2014 Martin Berger
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/