package edu.stanford.rsl.conrad.reconstruction.iterative;
import edu.stanford.rsl.conrad.data.numeric.Grid2D;
import edu.stanford.rsl.conrad.data.numeric.Grid3D;
import edu.stanford.rsl.conrad.numerics.SimpleMatrix;
import edu.stanford.rsl.conrad.reconstruction.iterative.IterativeReconstruction;
public class RayDrivenBasedReconstruction extends IterativeReconstruction {
private static final long serialVersionUID = 1L;
protected Grid3D projectionImage;
protected int lineOffset = 0;
protected final int maxU = getGeometry().getDetectorWidth(); //or it should be projection.getWidth();
protected final int maxV = getGeometry().getDetectorHeight();
protected final int maxK = getGeometry().getReconDimensionZ();
protected final int maxI = getGeometry().getReconDimensionX();
protected final int maxJ = getGeometry().getReconDimensionY();
protected final double dx = getGeometry().getVoxelSpacingX();
protected final double dy = getGeometry().getVoxelSpacingY();
protected final double dz = getGeometry().getVoxelSpacingZ();
int umin;
int umax;
int vmin;
int vmax;
double []weights;
SimpleMatrix mat;
@Override
public void prepareForSerialization(){
super.prepareForSerialization();
init = false;
}
public synchronized void Initialize(){
if (!init){
super.init();
//configuration?????
InitializeProjectionImage();
Grid2D projection = inputQueue.get(0);
//if projectionWidth does not equal to detectorWidth
//this needs to be further discussed.
if (getGeometry().getDetectorWidth() != -1){
System.out.println("row size projection: " + projection.getWidth() + "\nrow size detector: " + getGeometry().getDetectorWidth());
lineOffset = (projection.getWidth() - getGeometry().getDetectorWidth()) / 2; //have not considered if the lineOffset is not integer
}
copyProjectionImage();
}
}
public void InitializeProjectionImage() {
projectionImage = new Grid3D(nImages,maxU,maxV);
if (debug) System.out.println("Created Projection Image");
computeOffsets();
}
public void copyProjectionImage(){
Grid2D projection;
for ( int ip = 0; ip < nImages; ip++ ){
try {
projection = inputQueue.get(ip);
for (int iu = 0; iu <= maxU ; iu ++ ){
for (int iv = 0; iv <= maxV ; iv++){
// there may be a problem
projectionImage.setAtIndex(ip, iu, iv, projection.getPixelValue(iu + lineOffset, iv));
}
}
} catch (Exception e){
System.out.println("An error occured during backprojection of projection " + ip);
}
}
}
public void backproject() throws Exception {
for ( int p = 0; p < nImages; p++ ){
mat = getGeometry().getProjectionMatrix(p).computeP();
for( int k = 0; k < maxK; k++){
for( int i = 0; i < maxI; i++ ){
for( int j = 0; j < maxJ; j++){
if ( getSystemMatrixEntries( i, j, k) ){
//Add values to projectionImage
}
}
}
}
}
}
public void forwardproject() throws Exception {
for ( int p = 0; p < nImages; p++ ){
mat = getGeometry().getProjectionMatrix(p).computeP();
for( int k = 0; k < maxK; k++){
for( int i = 0; i < maxI; i++ ){
for( int j = 0; j < maxJ; j++){
if ( getSystemMatrixEntries( i, j, k) ){
//Add values to projectionImage
}
}
}
}
}
}
public void getSystemMatrixRow(){
//TODO
}
public void getSystemMatrixColumn(){
//TODO
}
// column/voxel index -> list of pixel indexes with weights
public boolean getSystemMatrixEntries( int i, int j, int k){
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];
if (mat != null){
updateMatrix[0][3] = mat.getElement(0,3);
updateMatrix[1][3] = mat.getElement(1,3);
updateMatrix[2][3] = mat.getElement(2,3);
voxel[3] = 1;
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];
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];
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];
double coordX = homogeniousPointj[0] / homogeniousPointj[2];
double coordY = homogeniousPointj[1] / homogeniousPointj[2];
umin = (int) Math.floor (coordX);
umax = (int) Math.ceil (coordX);
vmin = (int) Math.floor (coordY);
vmax = (int) Math.ceil (coordY);
umin = Math.max(umin, maxU);
umax = Math.min(umax, 0);
vmin = Math.max(vmin, maxV);
vmax = Math.min(vmax, 0);
if ( umax > umin && vmax > vmin){
return false;
}
//compute the weights
}
return true;
}
@Override
public String getBibtexCitation() {
// TODO Auto-generated method stub
return null;
}
@Override
public String getMedlineCitation() {
// TODO Auto-generated method stub
return null;
}
@Override
public String getToolName() {
// TODO Auto-generated method stub
return null;
}
@Override
public void iterativeReconstruct() throws Exception {
// TODO Auto-generated method stub
}
}
/*
* Copyright (C) 2010-2014 Meng Wu
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/