package edu.stanford.rsl.tutorial.dmip;
import edu.stanford.rsl.conrad.data.numeric.Grid1D;
import edu.stanford.rsl.conrad.data.numeric.Grid2D;
import edu.stanford.rsl.conrad.data.numeric.InterpolationOperators;
import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators;
import edu.stanford.rsl.conrad.filtering.redundancy.ParkerWeightingTool;
import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND;
import edu.stanford.rsl.conrad.geometry.shapes.simple.StraightLine;
import edu.stanford.rsl.conrad.numerics.SimpleOperators;
import edu.stanford.rsl.conrad.numerics.SimpleVector;
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;
import edu.stanford.rsl.tutorial.filters.SheppLoganKernel;
import ij.IJ;
import ij.ImageJ;
/**
* Exercise 7 of Diagnostic Medical Image Processing (DMIP)
* @author Marco Boegel
*
*/
public class DMIP_FanBeamBackProjector2D {
//source to detector distance
private double focalLength;
//increment of rotation angle beta between two projections in [rad]
private double betaIncrement;
//number of pixels on the detector
private int detectorPixels;
//number of acquired projections
private int numProjs;
//spacing - mm per pixel
private double detectorSpacing;
//detector length in [mm]
private double detectorLength;
//max rotation angle beta in [rad]
private double maxBeta;
public enum RampFilterType {NONE, RAMLAK, SHEPPLOGAN};
public DMIP_FanBeamBackProjector2D( )
{
}
/**
* Initialize all relevant parameters for the reconstruction
* @param sino the sinogram
* @param focalLength source to detector distance
* @param maxRot maximum rotation angle in [rad] at which the sinogram was acquired
*/
public void setSinogramParams(Grid2D sino, double focalLength, double maxRot)
{
this.focalLength = focalLength;
this.numProjs = sino.getHeight();
this.detectorPixels = sino.getWidth();
this.detectorSpacing = sino.getSpacing()[0];
this.detectorLength = detectorSpacing*detectorPixels;
double halfFanAngle = 0;//TODO
System.out.println("Half fan angle: " + halfFanAngle*180.0/Math.PI);
//TODO
this.betaIncrement = maxBeta /(double) numProjs;
System.out.println("Short-scan range: " + maxBeta*180/Math.PI);
}
/**
* A pixel driven backprojection algorithm. Cosine, Redundancy and Ramp filters need to be applied separately beforehand
* @param sino the filtered sinogram
* @param recoSize the dimension of the output image
* @param spacing the spacing of the output image
* @return the reconstruction
*/
public Grid2D backproject(Grid2D sino, int[] recoSize, double[] spacing)
{
Grid2D result = new Grid2D(recoSize[0], recoSize[1]);
result.setSpacing(spacing[0], spacing[1]);
for(int p = 0; p < numProjs; p++)
{
//First, compute the rotation angle beta and pre-compute cos(beta), sin(beta)
float beta = (float) (betaIncrement * p);
float cosBeta = (float) Math.cos(beta);
float sinBeta = (float) Math.sin(beta);
//Compute direction and normal of the detector at the current rotation angle
final PointND detBorder = new PointND();//TODO
final SimpleVector dirDet = new SimpleVector();//TODO
final StraightLine detLine = new StraightLine(detBorder, dirDet);
//Compute rotated source point
final PointND source = new PointND(focalLength * cosBeta, focalLength*sinBeta, 0.d);
//pick current projection
Grid1D currProj = sino.getSubGrid(p);
//pixel driven BP: iterate over all output image pixels
for(int x = 0; x < recoSize[0]; x++)
{
//transform the image pixel coordinates to world coordinates
float wx =0; //TODO
for(int y = 0; y < recoSize[1]; y++)
{
float wy = 0;//TODO
final PointND reconstructionPointWorld = new PointND(wx, wy, 0.d);
//intersect the projection ray with the detector
//TODO
final PointND detPixel = new PointND();//TODO
float valueTemp;
if(detPixel != null)
{
//calculate position of the projected point on the 1D detector
final SimpleVector pt = SimpleOperators.subtract(detPixel.getAbstractVector(), detBorder.getAbstractVector());
double length = pt.normL2();
//deal with truncation
if(pt.getElement(0)*dirDet.getElement(0)+pt.getElement(1)*dirDet.getElement(1)<0)
{
length -= length;
}
double t = (length - 0.5d)/detectorSpacing;
//check if projection of this world point hit the detector
if(currProj.getSize()[0] <= t+1 || t< 0)
continue;
float value = InterpolationOperators.interpolateLinear(currProj, t);
//Apply distance weighting
//see Fig 1a) exercise sheet
//TODO
//TODO
float dWeight = 0;//TODO
valueTemp = (float) (value / (dWeight*dWeight));
}
else
{
//projected pixel lies outside the detector
valueTemp = 0.f;
}
result.addAtIndex(x, y, valueTemp);
}
}
}
//adjust scale. required because of shortscan
float normalizationFactor = (float) (numProjs / Math.PI);
NumericPointwiseOperators.divideBy(result, normalizationFactor);
return result;
}
/**
* Cosine weighting for 2D sinograms
* @param sino the sinogram
* @return the cosine weighted sinogram
*/
public Grid2D applyCosineWeights(Grid2D sino)
{
Grid2D result = new Grid2D(sino);
//Create 1D kernel (identical for each projection)
Grid1D cosineKernel = new Grid1D(detectorPixels);
for(int i=0; i<detectorPixels; ++i){
//TODO
cosineKernel.setAtIndex(i, 9999);//TODO
}
//apply cosine weights to each projection
for(int p = 0; p < numProjs; p++)
{
NumericPointwiseOperators.multiplyBy(result.getSubGrid(p), cosineKernel);
}
return result;
}
/**
* Parker redundancy weights for a 2D sinogram
* @param sino the sinogram
* @return parker weighted sinogram
*/
public Grid2D applyParkerWeights(Grid2D sino)
{
Grid2D result = new Grid2D(sino);
Grid2D parker = new Grid2D(sino.getWidth(), sino.getHeight());
// Initialize parameters
double delta = (double) (Math.atan((detectorLength / 2.d)/ focalLength) );
double beta, gamma;
// iterate over the detector elements
for (int t = 0; t < detectorPixels; ++t) {
// compute alpha of the current ray (detector element)
gamma = Math.atan((t * detectorSpacing - detectorLength / 2.d + 0.5*detectorSpacing) / focalLength);
// iterate over the projection angles
for (int b = 0; b < numProjs; ++b) {
beta = b * betaIncrement;
// Shift weights such that they are centered (Important for maxBeta < pi + 2 * gammaM)
beta += (Math.PI+2*delta-maxBeta)/2.0;
// Adjust beta if out of range [0, 2*pi]
if (beta < 0) {
continue;
}
if (beta > Math.PI *2.d) {
continue;
}
// implement the conditions as described in Parker's paper
if (beta <= 2 * (delta - gamma)) {
float val = 0; //TODO
if (Double.isNaN(val)){
continue;
}
parker.setAtIndex(t, b , val);
} else if (beta < Math.PI - 2.d * gamma) {
parker.setAtIndex(t, b , 1);
}
else if (beta <= (Math.PI + 2.d * delta) + 1e-12) {
float val = 0;//TODO
if (Double.isNaN(val)){
continue;
}
parker.setAtIndex(t, b , val);
}
}
}
// Correct for scaling due to varying angle
NumericPointwiseOperators.multiplyBy(parker, (float)( maxBeta / (Math.PI)));
parker.show();
//apply the parker weights to the sinogram
NumericPointwiseOperators.multiplyBy(result, parker);
return result;
}
public static void main(String arg[])
{
ImageJ ij = new ImageJ();
double focalLength = 600.f;
//maximum rotation angle in [rad]
double maxRot = Math.PI;
//choose the ramp filter (none, ramlak, shepplogan)
RampFilterType filter = RampFilterType.RAMLAK;
DMIP_FanBeamBackProjector2D fbp = new DMIP_FanBeamBackProjector2D();
//Load and visualize the projection image data
String filename = "D:/02_lectures/DMIP/exercises/2014/6/Sinogram0.tif";
Grid2D sino = ImageUtil.wrapImagePlus(IJ.openImage(filename)).getSubGrid(0);
sino.show("Sinogram");
//getSubGrid() only yields rows. -> Transpose so that each row is a projection
sino = sino.getGridOperator().transpose(sino);
sino.show();
//Initialize parameters
double detectorSpacing = 0.5;
sino.setSpacing(detectorSpacing, 0);
fbp.setSinogramParams(sino, focalLength, maxRot);
sino.setSpacing(fbp.detectorSpacing, fbp.betaIncrement);
//apply cosine weights
Grid2D cosineWeighted = fbp.applyCosineWeights(sino);
cosineWeighted.show("After Cosine Weighting");
//apply parker redundancy weights
Grid2D parkerWeighted = fbp.applyParkerWeights(cosineWeighted);
parkerWeighted.show("After Parker Weighting");
//apply ramp filter
switch(filter){
case RAMLAK:
RamLakKernel ramLak = new RamLakKernel(fbp.detectorPixels, fbp.detectorSpacing);
for(int i = 0; i < fbp.numProjs; i++)
{
ramLak.applyToGrid(parkerWeighted.getSubGrid(i));
}
parkerWeighted.show("After RamLak Filter");
break;
case SHEPPLOGAN:
SheppLoganKernel sheppLogan = new SheppLoganKernel(fbp.detectorPixels, fbp.detectorSpacing);
for(int i = 0; i < fbp.numProjs; i++)
{
sheppLogan.applyToGrid(parkerWeighted.getSubGrid(i));
}
parkerWeighted.show("After Shepp-Logan Filter");
break;
case NONE:
default:
}
//setup reconstruction image size
int[] recoSize = new int[]{128, 128};
double[] spacing = new double[]{1.0, 1.0};
Grid2D reconstruction = fbp.backproject(parkerWeighted, recoSize, spacing);
reconstruction.show("Reconstructed image");
//Compare our backprojection algorithm to the tutorial code
FanBeamBackprojector2D fbp2 = new FanBeamBackprojector2D(focalLength, fbp.detectorSpacing, fbp.betaIncrement,recoSize[0], recoSize[1]);
fbp2.backprojectPixelDriven(parkerWeighted).show("Tutorial reconstruction");
}
}