package edu.stanford.rsl.conrad.reconstruction;
import edu.stanford.rsl.conrad.data.numeric.Grid2D;
import edu.stanford.rsl.conrad.numerics.SimpleMatrix;
import ij.process.FloatProcessor;
/**
* Implementation of the backprojector which is required for Noo's super-short-scan reconstruction. Note that the distance weight in the weighted sum of the backprojection is different.
* @author akmaier
*
*/
public class CPUSuperShortScanBackprojector extends VOIBasedReconstructionFilter {
/**
*
*/
private static final long serialVersionUID = -1248364198978750688L;
protected boolean mapped = false;
@Override
public void prepareForSerialization(){
super.prepareForSerialization();
mapped = false;
}
@Override
public void backproject(Grid2D projection, int projectionNumber) {
int count = 0;
//System.out.println(projectionVolume);
if (!init){
initialize(projection);
}
FloatProcessor currentProjection = new FloatProcessor(projection.getWidth(), projection.getHeight(), projection.getBuffer(), null);
//ImageProcessor currentProjection = projection;
int p = projectionNumber;
double[] voxel = new double [4];
double[] homogeniousPointi = new double[3];
double[] homogeniousPointj = new double[3];
double[] homogeniousPointk = new double[3];
double[][] updateMatrix = new double [3][4];
SimpleMatrix mat = getGeometry().getProjectionMatrix(p).computeP();
if (mat !=null) {
//mat.print(NumberFormat.getInstance(), 6);
voxel[3] = 1;
updateMatrix[0][3] = mat.getElement(0,3);
updateMatrix[1][3] = mat.getElement(1,3);
updateMatrix[2][3] = mat.getElement(2,3);
boolean nanHappened = false;
double min = 20000;
double max = 0;
for (int k = 0; k < maxK ; k++){ // for all slices
if (debug) System.out.println("here: " + " " + k);
voxel[2] = (this.getGeometry().getVoxelSpacingZ() * k) - offsetZ;
updateMatrix[0][2] = mat.getElement(0,2) * voxel[2];
updateMatrix[1][2] = mat.getElement(1,2) * voxel[2];
updateMatrix[2][2] = mat.getElement(2,2) * voxel[2];
homogeniousPointk[0] = updateMatrix[0][3] + updateMatrix[0][2];
homogeniousPointk[1] = updateMatrix[1][3] + updateMatrix[1][2];
homogeniousPointk[2] = updateMatrix[2][3] + updateMatrix[2][2];
for (int i=0; i < maxI; i++){ // for all lines
voxel[0] = (this.getGeometry().getVoxelSpacingX() * i) - offsetX;
updateMatrix[0][0] = mat.getElement(0,0) * voxel[0];
updateMatrix[1][0] = mat.getElement(1,0) * voxel[0];
updateMatrix[2][0] = mat.getElement(2,0) * voxel[0];
homogeniousPointi[0] = homogeniousPointk[0] + updateMatrix[0][0];
homogeniousPointi[1] = homogeniousPointk[1] + updateMatrix[1][0];
homogeniousPointi[2] = homogeniousPointk[2] + updateMatrix[2][0];
for (int j = 0; j < maxJ; j++){ // for all voxels
// compute real world coordinates in homogenious coordinates;
boolean project = true;
if (useVOImap){
if (voiMap != null){
project = voiMap[i][j][k];
}
}
if (project){
voxel[1] = (this.getGeometry().getVoxelSpacingY() * j) - offsetY;
updateMatrix[0][1] = mat.getElement(0,1) * voxel[1];
updateMatrix[1][1] = mat.getElement(1,1) * voxel[1];
updateMatrix[2][1] = mat.getElement(2,1) * voxel[1];
homogeniousPointj[0] = homogeniousPointi[0] + updateMatrix[0][1];
homogeniousPointj[1] = homogeniousPointi[1] + updateMatrix[1][1];
homogeniousPointj[2] = homogeniousPointi[2] + updateMatrix[2][1];
//Matrix point3d = new Matrix(voxel);
//Compute coordinates in projection data.
//Matrix point2d = geometry.getProjectionMatrix(p).times(point3d);
double coordX = homogeniousPointj[0] / homogeniousPointj[2];
double coordY = homogeniousPointj[1] / homogeniousPointj[2];
// back project
double increment = currentProjection.getInterpolatedValue(coordX + lineOffset, coordY);//*homogeniousPointj[2]);
if (Double.isNaN(increment)){
nanHappened = true;
if (count < 10) System.out.println("NAN Happened at i = " + i + " j = " + j + " k = " + k + " projection = " + projectionNumber + " x = " + coordX + " y = " + coordY );
increment = 0;
count ++;
}
// Weighting according to super-short-scan formula for cone beam geometry.
// homogeniousPointj[2] => ||a- a(lambda)
// conversion of coordX (the element index) back to metric detector coordinates (\tilde{u})required:
// \tilde{u} = (coordX-(this.detectorWidth/2))*this.detectorElementSizeX
//homogeniousPointj[2] /= this.geometry.getSourceToDetectorDistance();
//homogeniousPointj[2] *= homogeniousPointj[2];
//homogeniousPointj[2] = geometry.getSourceToDetectorDistance() - homogeniousPointj[2];
min = (homogeniousPointj[2]< min)? homogeniousPointj[2]: min;
max = (homogeniousPointj[2]> max)? homogeniousPointj[2]: max;
double utilde = (coordX-(this.getGeometry().getDetectorWidth()/2))*this.getGeometry().getPixelDimensionX();
//double weight = Math.sqrt(Math.pow(utilde,2) + Math.pow(geometry.getSourceToDetectorDistance(), 2))/ (Math.pow(geometry.getSourceToDetectorDistance(),1)*Math.abs(homogeniousPointj[2])*Math.abs(homogeniousPointj[2]));
double weight = Math.sqrt(Math.pow(utilde,2) + Math.pow(getGeometry().getSourceToDetectorDistance(), 2))/ (Math.pow(getGeometry().getSourceToDetectorDistance(),1) * homogeniousPointj[2]);
// Not bad:
// double weight = 1.0/ (200+ (Math.sqrt(Math.pow(((maxI/2)-i)*this.geometry.getVoxelSpacingX(),2)+ Math.pow(((maxJ/2)-j)*this.geometry.getVoxelSpacingY(), 2))));
//double weight = 1.0 / (Math.pow(homogeniousPointj[2],1));
increment *= weight ;//* this.parkerLikeWeight(primaryAngles[projectionNumber] /180.0 * Math.PI, utilde);
updateVolume(i, j, k, increment);
}
}
}
System.out.println("Min max: " + min + " " + max);
}
if (nanHappened) {
throw new RuntimeException("Encountered NaN in projection!");
}
if (debug) System.out.println("done with projection");
}
}
@Override
public String getName(){
return "Specialized Cone-beam Super Short Scan Backprojector";
}
@Override
public String getBibtexCitation() {
String bibtex = "@ARTICLE{Noo02-IRF,\n" +
" author = {Noo, F. and Defrise, M. and Clackdoyle, R. and Kudo, H.},\n" +
" title = \"{{Image reconstruction from fan-beam projections on less than a short scan}}\",\n" +
" journal = {Physics in Medicine and Biology},\n" +
" year = 2002,\n" +
" volume = 47,\n"+
" number = 14,\n" +
" pages = {2525-2546}\n" +
"}";
return bibtex;
}
@Override
public String getMedlineCitation(){
return "Noo F, Defrise M, Clackdoyle R, Kudo H. Image reconstruction from fan-beam projections on less than a short scan. " +
"Phys Med Biol 47(14):2525-46. 2002.";
}
@Override
public String getToolName(){
return "Super Short Scan Backprojector";
}
}
/*
* Copyright (C) 2010-2014 Andreas Maier
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/