package edu.stanford.rsl.tutorial.parallel;
import java.util.ArrayList;
import java.util.Arrays;
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.geometry.shapes.simple.Box;
import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND;
import edu.stanford.rsl.conrad.geometry.shapes.simple.StraightLine;
import edu.stanford.rsl.conrad.geometry.transforms.Transform;
import edu.stanford.rsl.conrad.geometry.transforms.Translation;
import edu.stanford.rsl.conrad.numerics.SimpleOperators;
import edu.stanford.rsl.conrad.numerics.SimpleVector;
/**
* Implementation of a simple parallel backprojector. In order to create the
* image, the angular range and the angular sampling and the detector element
* size and detector element sampling has to be defined. We show both, a ray
* driven and a pixel driven backprojector.
*
* See L. Zeng. "Medical Image Reconstruction: A Conceptual tutorial". 2009, page 8
*
* @author Recopra Seminar Summer 2012
*
*/
public class ParallelBackprojector2D {
final double samplingRate = 3.d;
int maxThetaIndex, maxSIndex; // image dimensions [GU]
double maxTheta, deltaTheta, // detector angles [rad]
maxS, deltaS; // detector (pixel) size [mm]
int imgSizeX, imgSizeY; // image dimensions [GU]
float pxSzXMM, pxSzYMM; // pixel size [mm]
/**
* Sampling of projections is defined in the constructor.
*
* @param imageSizeX
* @param imageSizeY
* @param pxSzXMM
* @param pxSzYMM */
public ParallelBackprojector2D(int imageSizeX, int imageSizeY, float pxSzXMM, float pxSzYMM){
this.imgSizeX = imageSizeX;
this.imgSizeY = imageSizeY;
this.pxSzXMM = pxSzXMM;
this.pxSzYMM = pxSzYMM;
}
void initSinogramParams(Grid2D sino) {
this.maxThetaIndex = sino.getSize()[1];
this.deltaTheta = sino.getSpacing()[1];
this.maxTheta = (maxThetaIndex -1) * deltaTheta;
this.maxSIndex = sino.getSize()[0];
this.deltaS = sino.getSpacing()[0];
this.maxS = (maxSIndex-1) * deltaS;
}
/**
* The ray driven solution.
*
* @param sino
* the sinogram
* @return the image
*/
public Grid2D backprojectRayDriven(Grid2D sino) {
this.initSinogramParams(sino);
Grid2D grid = new Grid2D(this.imgSizeX, this.imgSizeY);
grid.setSpacing(this.pxSzXMM, this.pxSzYMM);
grid.setOrigin(-(grid.getSize()[0]*grid.getSpacing()[0]) / 2.0, -(grid.getSize()[1]*grid.getSpacing()[1]) / 2.0);
// set up image bounding box in WC
Translation trans = new Translation(
-(grid.getSize()[0] * grid.getSpacing()[0])/2, -(grid.getSize()[1] * grid.getSpacing()[1])/2, -1
);
Transform inverse = trans.inverse();
Box b = new Box((grid.getSize()[0] * grid.getSpacing()[0]), (grid.getSize()[1] * grid.getSpacing()[1]), 2);
b.applyTransform(trans);
for(int e=0; e<maxThetaIndex; ++e){
// compute theta [rad] and angular functions.
double theta = deltaTheta * e;
double cosTheta = Math.cos(theta);
double sinTheta = Math.sin(theta);
for (int i = 0; i < maxSIndex; ++i) {
// compute s, the distance from the detector edge in WC [mm]
double s = deltaS * i - maxS / 2;
// compute two points on the line through s and theta
// We use PointND for Points in 3D space and SimpleVector for directions.
PointND p1 = new PointND(s * cosTheta, s * sinTheta, .0d);
PointND p2 = new PointND(-sinTheta + (s * cosTheta),
(s * sinTheta) + cosTheta, .0d);
// set up line equation
StraightLine line = new StraightLine(p1, p2);
// compute intersections between bounding box and intersection line.
ArrayList<PointND> points = b.intersect(line);
// only if we have intersections
if (2 != points.size()){
if(points.size() == 0) {
line.getDirection().multiplyBy(-1.d);
points = b.intersect(line);
}
if(points.size() == 0)
continue;
}
PointND start = points.get(0); // [mm]
PointND end = points.get(1); // [mm]
// get the normalized increment
SimpleVector increment = new SimpleVector(end.getAbstractVector());
increment.subtract(start.getAbstractVector());
double distance = increment.normL2();
increment.divideBy(distance * samplingRate);
float val = sino.getAtIndex(i, e);
start = inverse.transform(start);
// compute the integral along the line.
for (double t = 0.0; t < distance * samplingRate; ++t) {
PointND current = new PointND(start);
current.getAbstractVector().add(increment.multipliedBy(t));
double x = current.get(0) / grid.getSpacing()[0],
y = current.get(1) / grid.getSpacing()[1];
if (x + 1 > grid.getSize()[0]
|| y + 1 >= grid.getSize()[1]
|| x < 0 || y < 0)
continue;
InterpolationOperators.addInterpolateLinear(grid, x, y, val);
}
}
}
// float correctionFactor = (float) ((float) samplingRate * maxThetaIndex / maxTheta * Math.PI);
float normalizationFactor = (float) ((float) samplingRate * maxThetaIndex / deltaS / Math.PI);
NumericPointwiseOperators.divideBy(grid, normalizationFactor);
return grid;
}
/**
* The pixel driven solution.
*
* @param sino
* the sinogram
* @return the image
*/
public Grid2D backprojectPixelDriven(Grid2D sino) {
this.initSinogramParams(sino);
Grid2D grid = new Grid2D(this.imgSizeX, this.imgSizeY);
grid.setSpacing(pxSzXMM, pxSzYMM);
grid.setOrigin(-(grid.getSize()[0]*grid.getSpacing()[0])/2, -(grid.getSize()[1]*grid.getSpacing()[1])/2);
// loop over the projection angles
for (int i = 0; i < maxThetaIndex; i++) {
// compute actual value for theta
double theta = deltaTheta * i;
// precompute sine and cosines for faster computation
double cosTheta = Math.cos(theta);
double sinTheta = Math.sin(theta);
// get detector direction vector
SimpleVector dirDetector = new SimpleVector(sinTheta,cosTheta);
// loops over the image grid
for (int x = 0; x < grid.getSize()[0]; x++) {
for (int y = 0; y < grid.getSize()[1]; y++) {
// compute world coordinate of current pixel
double[] w = grid.indexToPhysical(x, y);
// wrap into vector
SimpleVector pixel = new SimpleVector(w[0], w[1]);
// project pixel onto detector
double s = SimpleOperators.multiplyInnerProd(pixel, dirDetector);
// compute detector element index from world coordinates
s += maxS/2; // [mm]
s /= deltaS; // [GU]
// get detector grid
Grid1D subgrid = sino.getSubGrid(i);
// check detector bounds, continue if out of array
if (subgrid.getSize()[0] <= s + 1
|| s < 0)
continue;
// get interpolated value
float val = InterpolationOperators.interpolateLinear(subgrid, s);
// sum value to sinogram
grid.addAtIndex(x, y, val);
}
}
}
// apply correct scaling
NumericPointwiseOperators.divideBy(grid, (float) (maxThetaIndex / Math.PI));
return grid;
}
/**
* The pixel driven solution. for one angle
*
* @param sino
* the sinogram
* @return the image
*/
public Grid2D backprojectPixelDriven(Grid2D sino, int i) {
this.initSinogramParams(sino);
Grid2D grid = new Grid2D(this.imgSizeX, this.imgSizeY);
grid.setSpacing(pxSzXMM, pxSzYMM);
grid.setOrigin(-(grid.getSize()[0]*grid.getSpacing()[0])/2, -(grid.getSize()[1]*grid.getSpacing()[1])/2);
// compute actual value for theta
double theta = deltaTheta * i;
// precompute sine and cosines for faster computation
double cosTheta = Math.cos(theta);
double sinTheta = Math.sin(theta);
// get detector direction vector
SimpleVector dirDetector = new SimpleVector(cosTheta,sinTheta);
// loops over the image grid
for (int x = 0; x < grid.getSize()[0]; x++) {
for (int y = 0; y < grid.getSize()[1]; y++) {
// compute world coordinate of current pixel
double[] w = grid.indexToPhysical(x, y);
// wrap into vector
SimpleVector pixel = new SimpleVector(w[0], w[1]);
// project pixel onto detector
double s = SimpleOperators.multiplyInnerProd(pixel, dirDetector);
// compute detector element index from world coordinates
s += maxS/2; // [mm]
s /= deltaS; // [GU]
// get detector grid
Grid1D subgrid = sino.getSubGrid(i);
// check detector bounds, continue if out of array
if (subgrid.getSize()[0] <= s + 1
|| s < 0)
continue;
// get interpolated value
float val = InterpolationOperators.interpolateLinear(subgrid, s);
// sum value to sinogram
grid.addAtIndex(x, y, val);
}
}
// apply correct scaling
NumericPointwiseOperators.divideBy(grid, (float) (maxThetaIndex / Math.PI));
return grid;
}
/**
* The pixel driven solution with median filter on the backprojection valeus.
*
* @param sino
* the sinogram
* @return the image
*/
public Grid2D backprojectPixelDrivenPercentile(Grid2D sino, double percentile) {
this.initSinogramParams(sino);
Grid2D grid = new Grid2D(this.imgSizeX, this.imgSizeY);
grid.setSpacing(pxSzXMM, pxSzYMM);
grid.setOrigin(-(grid.getSize()[0]*grid.getSpacing()[0])/2, -(grid.getSize()[1]*grid.getSpacing()[1])/2);
SimpleVector [] directionVectors = new SimpleVector[maxThetaIndex];
int start = (int) (maxThetaIndex *percentile);
int end = maxThetaIndex -start;
float range = end - start;
// precompute the direction vectors for faster computation
for (int i = 0; i < maxThetaIndex; i++) {
// compute actual value for theta
double theta = deltaTheta * i;
// precompute sine and cosines for faster computation
double cosTheta = Math.cos(theta);
double sinTheta = Math.sin(theta);
// get detector direction vector
directionVectors[i] = new SimpleVector(sinTheta,cosTheta);
}
// loop over the reconstruction grid
for (int x = 0; x < grid.getSize()[0]; x++) {
for (int y = 0; y < grid.getSize()[1]; y++) {
// loop over the projection angles
double [] components = new double[maxThetaIndex];
for (int i = 0; i < maxThetaIndex; i++) {
// compute world coordinate of current pixel
double[] w = grid.indexToPhysical(x, y);
// wrap into vector
SimpleVector pixel = new SimpleVector(w[0], w[1]);
// project pixel onto detector
double s = SimpleOperators.multiplyInnerProd(pixel, directionVectors[i]);
// compute detector element index from world coordinates
s += maxS/2; // [mm]
s /= deltaS; // [GU]
// get detector grid
Grid1D subgrid = sino.getSubGrid(i);
// check detector bounds, continue if out of array
if (subgrid.getSize()[0] <= s + 1
|| s < 0)
continue;
// get interpolated value
components[i] = InterpolationOperators.interpolateLinear(subgrid, s);
// sum value to sinogram
}
// sort array
Arrays.sort(components);
// backproject median value
for (int i = start; i <= end; i++) {
grid.addAtIndex(x, y, (float)components[i]);
}
}
}
// apply correct scaling
NumericPointwiseOperators.divideBy(grid, (float) (range / Math.PI));
return grid;
}
}
/*
* Copyright (C) 2010-2014 Andreas Maier
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/