package edu.stanford.rsl.conrad.reconstruction.iterative;
import edu.stanford.rsl.conrad.data.numeric.Grid3D;
import edu.stanford.rsl.conrad.data.numeric.NumericGrid;
import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators;
import edu.stanford.rsl.conrad.geometry.trajectories.Trajectory;
public class LeastSquaresCG extends DistanceDrivenBasedReconstruction {
/**
* Model-based iterative reconstruction using least-squares conjugate gradient solver.
* by Meng Wu
*/
private static final long serialVersionUID = 1L;
static final int maxNumOfIterations = 20;
public LeastSquaresCG(Trajectory dataTrajectory) {
super(dataTrajectory);
// TODO Auto-generated constructor stub
}
public void initialize( Grid3D proj){
maxI = getGeometry().getReconDimensionX();
maxJ = getGeometry().getReconDimensionY();
maxK = getGeometry().getReconDimensionZ();
maxU = getGeometry().getDetectorWidth(); //or it should be projection.getWidth();
maxV = getGeometry().getDetectorHeight();
dx = getGeometry().getVoxelSpacingX();
dy = getGeometry().getVoxelSpacingY();
dz = getGeometry().getVoxelSpacingZ();
time = System.currentTimeMillis();
projectionViews = proj;
//projectionViews = InitializeProjectionViews();
//volumeImage = InitializeVolumeImage();
}
public void iterativeReconstruct() throws Exception{
Grid3D r = InitializeProjectionViews();
Grid3D q = InitializeProjectionViews();
Grid3D f = InitializeVolumeImage();
Grid3D d = InitializeVolumeImage();
Grid3D g_new = InitializeVolumeImage();
Grid3D g_old = InitializeVolumeImage();
double gamma;
double alpha;
NumericPointwiseOperators.copy(projectionViews, r );
System.out.print("Model-based iterative reconstruction using least-squares conjugate gradient solver: \n");
System.out.print("itn \t||r|| \t\t x(0) \t\n" );
System.out.print("---------------------------------------- \n");
for ( int itr = 1; itr <= maxNumOfIterations; itr++ ){
if ( itr <= 5 || itr%10 == 0 || itr >= maxNumOfIterations - 5 )
System.out.print( itr + " \t" + String.format( "%8.3g", NumericPointwiseOperators.dotProduct( r )) + " \t"
+ String.format( "%8.4f", f.getAtIndex(maxJ/2, maxK/2, maxK/2)) + " \n" );
backproject( r, g_new );
if ( itr == 1 ){
gamma = 0.0;
NumericPointwiseOperators.copy( g_new, d );
}else{
gamma = NumericPointwiseOperators.dotProduct( g_new ) / NumericPointwiseOperators.dotProduct( g_old );
NumericPointwiseOperators.multiplyBy(d, (float) gamma);
NumericPointwiseOperators.addBy(d, g_new);
}
forwardproject( q, d );
alpha = NumericPointwiseOperators.dotProduct( g_new, d ) / NumericPointwiseOperators.dotProduct( q );
NumericGrid tmp = d.clone();
NumericPointwiseOperators.multiplyBy(tmp, (float) alpha);
NumericPointwiseOperators.addBy(f, tmp);
tmp = q.clone();
NumericPointwiseOperators.multiplyBy(tmp, (float) -alpha);
NumericPointwiseOperators.addBy(r, tmp);
NumericPointwiseOperators.copy( g_new, g_old);
}
volumeImage = f;
}
}
/*
* Copyright (C) 2010-2014 Meng Wu
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/