content: package edu.stanford.rsl.conrad.reconstruction;
  
  /**
   * @author febmeng
   *
   */
  
  import edu.stanford.rsl.conrad.data.numeric.Grid3D;
  import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators;
  import edu.stanford.rsl.conrad.geometry.Projection;
  import edu.stanford.rsl.conrad.geometry.trajectories.Trajectory;
  import edu.stanford.rsl.conrad.numerics.SimpleMatrix;
  import edu.stanford.rsl.conrad.numerics.SimpleOperators;
  import edu.stanford.rsl.conrad.numerics.SimpleVector;
  
  
  public class DistanceDrivenBasedReconstruction  extends ModelBasedIterativeReconstruction {
  
  
  	private static final long serialVersionUID = 1L;
  	public boolean Debug = true;
  	public boolean Debug1 = false;
  	public boolean Debug2 = false;
  	protected static final int MAX_WEIGHT_LENGTH = 6;
  
  	public DistanceDrivenBasedReconstruction( Trajectory dataTrajectory ) {
  		super(dataTrajectory);
  		// TODO Auto-generated constructor stub
  	}
  
  
  	@Override
  	public void iterativeReconstruct() throws Exception {
  
  
  		forwardproject( projectionViews, volumeImage );
  	}
  
  	@Override
  	protected void forwardproject(Grid3D projImage, Grid3D volImage) throws Exception {
  
  		//zero out whole volume image		
  		NumericPointwiseOperators.fill(projImage, 0.0f);
  		NumericPointwiseOperators.fill(volImage, 1.0f);
  		//if (Debug) volImage.printOneSlice(3);
  		Projection proj;
  
  		for ( int p = 0; p < nImages ; p++ ){
  			proj = getGeometry().getProjectionMatrix(p);
  			distanceDrivenProjView( projImage, volImage, proj, p );
  		}
  
  		//if (Debug) projImage.printOneSlice(8);
  	}
  
  
  	@Override
  	protected void backproject(Grid3D projImage, Grid3D volImage) throws Exception {
  
  		//zero out whole volume image
  		NumericPointwiseOperators.fill(volImage, 0.0f);
  
  		//if (Debug) volImage.printOneSlice(3);
  		Projection proj;
  
  		for ( int p = 0; p < nImages; p++ ){
  			proj = getGeometry().getProjectionMatrix(p);
  			distanceDrivenBackView( projImage, volImage, proj, p );
  		}
  
  	}
  
  	protected void distanceDrivenProjView( Grid3D projImage, Grid3D volImage, Projection proj, final int ip ){
  
  
  		SimpleMatrix mat = proj.computeP();
  		SimpleVector cameraCenter = proj.computeCameraCenter();
  
  		// buffer have move of voxel
  		SimpleVector halfVoxelMoveX = mat.getCol(0);
  		halfVoxelMoveX.multiplyBy(dx/2);
  		SimpleVector halfVoxelMoveY = mat.getCol(1);
  		halfVoxelMoveY.multiplyBy(dy/2);
  		SimpleVector fullVoxelMoveZ = mat.getCol(2);
  		fullVoxelMoveZ.multiplyBy(dz);
  
  		SimpleVector point3d = new SimpleVector(4);
  		point3d.setElementValue(3, 1);
  
  		SimpleVector point2d, point2dMinus, point2dPlus;
  
  		double cx = cameraCenter.getElement(0);
  		double cy = cameraCenter.getElement(1);
  		double cz = cameraCenter.getElement(2);
  		double dsx0, dsy0, dsxy0_sqr;
  		double coordLeft, coordRight, coordBottom, coordTop, coordStep;
  		double ds0, amplitude;
  		int iumin, iumax, ivmin, ivmax;
  
  		if (Debug1) {
  			System.out.println( "Camera Center: " + cx + ", " + cy);
  			System.out.println( "Voxel Spaceing: " + dx + ", " + dy);
  			System.out.println( "Projection matrix: ");
  			printSimpleMatrix(mat);
  			System.out.println( "Half voxel moves: ");
  			printSimpleVector(halfVoxelMoveX);
  			printSimpleVector(halfVoxelMoveY);
  		} 
  
  		for (int i=1; i < maxI; i++){
  
  			point3d.setElementValue(0, i*dx - offsetX  );
  
  			for (int j = 0; j < maxJ; j++){
  				point3d.setElementValue(1, j*dy - offsetY );
  				point3d.setElementValue(2, -dz/2-offsetZ);
  				point2d = SimpleOperators.multiply(mat, point3d);
  
  				dsx0 = Math.abs(point3d.getElement(0) - cx);
  				dsy0 = Math.abs(point3d.getElement(1) - cy);
  				dsxy0_sqr = dsx0*dsx0 + dsy0*dsy0;
  
  				if (Debug1){
  					System.out.println("dsx0="+ dsx0 + ", dsy0=" + dsy0 + ", dsxy0_sqr= " + dsxy0_sqr);
  				}
  
  				if ( dsy0 > dsx0){
  					ds0 = dsy0;
  					point2dMinus = SimpleOperators.subtract(point2d, halfVoxelMoveX);
  					point2dPlus = SimpleOperators.add(point2d, halfVoxelMoveX);					
  				}else{
  					ds0 = dsx0;
  					point2dMinus = SimpleOperators.subtract(point2d, halfVoxelMoveY);
  					point2dPlus = SimpleOperators.add(point2d, halfVoxelMoveY);
  				}
  
  				coordLeft =  point2dMinus.getElement(0) / point2dMinus.getElement(2) + 0.5;
  				coordRight =  point2dPlus.getElement(0) / point2dPlus.getElement(2) + 0.5;
  
  				if ( coordLeft > coordRight ){
  					double temp = coordRight;
  					coordRight = coordLeft;
  					coordLeft = temp;	
  				}
  
  				coordLeft = Math.max(coordLeft, 0);
  				coordRight = Math.min(coordRight, maxU);
  
  				if (coordLeft >= coordRight){
  					continue;
  				}
  
  				iumin = (int)coordLeft;
  				iumax = (int)coordRight;
  				iumax = Math.min(iumax, maxU-1);
  
  				iumax = Math.min(iumax, iumin + MAX_WEIGHT_LENGTH - 1);
  
  				double[] weightU = new double[]{1,1,1,1,1,1};
  
  				if ( iumax == iumin){
  					weightU[0] = coordRight - coordLeft;
  				}else{
  					weightU[0] = iumin + 1 - coordLeft;
  					weightU[iumax - iumin] = coordRight - iumax;
  				};
  
  
  				coordBottom = point2d.getElement(1) / point2d.getElement(2) + 0.5;
  				point2d.add(fullVoxelMoveZ);
  				coordTop = point2d.getElement(1) / point2d.getElement(2) + 0.5;
  				coordStep = coordTop - coordBottom;
  
  
  				for ( int k = 0; k < maxK; k++ ){
  
  					double dsz0 = Math.abs(k*dz-offsetZ-cz); 
  
  					if ( coordBottom <= 0 ){
  						coordBottom = coordTop;
  						coordTop = coordTop + coordStep;
  						continue;
  					}
  
  					if ( coordTop >= maxV ){
  						break;
  					}
  
  					ivmin = (int)coordBottom;
  					ivmax = (int)coordTop;
  
  					if (Debug2) {
  						System.out.println( "iu: iumin= " + iumin + ", \t iumax=" + iumax);
  						System.out.println("iv: ivmin= " + ivmin + ", \t ivmax=" + ivmax );
  					}
  
  					double[] weightV = new double[]{1,1,1,1,1,1};
  
  					//ivmax = Math.min(ivmax, ivmin + MAX_WEIGHT_LENGTH - 1);			
  					if ( ivmin == ivmax ){
  						weightV[0] = coordTop - coordBottom;
  					}else{
  						weightV[0] = ivmin + 1 - coordBottom;
  						weightV[ivmax - ivmin] = coordTop - ivmax;
  					};
  
  					amplitude = Math.sqrt(dsz0*dsz0 + dsxy0_sqr) / ds0;
  
  					if (Debug2){
  						int[] detSize = projImage.getSize();
  						int[] volSize = volImage.getSize();
  						//int[] detBound = projImage.getBoundarySize();
  						System.out.println("Detector Size: " + detSize[0] + " X " + detSize[1] + " X " + detSize[2] );
  						//System.out.println("Detector Bound: " + detBound[0] + " X " + detBound[1] + " X " + detBound[2] );
  						System.out.println("Volume Size: "+ volSize[0] + " X " + volSize[1] + " X " + volSize[2] );
  						for (int u = 0 ; u <detSize[1]; u++ ){
  							System.out.println("u=" + u );
  							for (int v = 0; v < detSize[2]; v++){
  								System.out.print(projImage.getAtIndex(0, u, v) + "\t");
  							}
  							System.out.println();
  						}
  					}
  
  					double tempVal = volImage.getAtIndex(i, j, k) * amplitude ;
  
  					for ( int iu = iumin, iiu = 0; iu <= iumax; iu++, iiu++ ){
  						double temp = tempVal * weightU[iiu];
  						for (int iv = ivmin, iiv = 0; iv <= ivmax; iv++, iiv++){
  							projImage.addAtIndex(ip, iu, iv, (float)(temp*weightV[iiv]));
  						} //iv
  					} //iu
  
  					coordBottom = coordTop;
  					coordTop = coordTop + coordStep;
  
  				} //k
  			}//j		
  		}//i
  	}
  
  	protected void distanceDrivenBackView( Grid3D projImage, Grid3D volImage, Projection proj, final int ip ){
  
  		SimpleMatrix mat = proj.computeP();
  		SimpleVector cameraCenter = proj.computeCameraCenter();
  
  		// buffer have move of voxel
  		SimpleVector halfVoxelMoveX = mat.getCol(0);
  		halfVoxelMoveX.multiplyBy(dx/2);
  		SimpleVector halfVoxelMoveY = mat.getCol(1);
  		halfVoxelMoveY.multiplyBy(dy/2);
  		SimpleVector fullVoxelMoveZ = mat.getCol(2);
  		fullVoxelMoveZ.multiplyBy(dz);
  
  		SimpleVector point3d = new SimpleVector(4);
  		point3d.setElementValue(3, 1);
  
  		SimpleVector point2d, point2dMinus, point2dPlus;
  
  		double cx = cameraCenter.getElement(0);
  		double cy = cameraCenter.getElement(1);
  		double cz = cameraCenter.getElement(2);
  		double dsx0, dsy0, dsxy0_sqr;
  		double coordLeft, coordRight, coordBottom, coordTop, coordStep;
  		double ds0, amplitude;
  		int iumin, iumax, ivmin, ivmax;
  
  		if (Debug) {
  			System.out.println( "Camera Center: " + cx + ", " + cy);
  			System.out.println( "Voxel Spaceing: " + dx + ", " + dy);
  			System.out.println( "Projection matrix: ");
  			printSimpleMatrix(mat);
  			System.out.println( "Half voxel moves: ");
  			printSimpleVector(halfVoxelMoveX);
  			printSimpleVector(halfVoxelMoveY);
  		} 
  
  		for (int i=1; i < maxI; i++){
  
  			point3d.setElementValue(0, i*dx - offsetX  );
  
  			for (int j = 0; j < maxJ; j++){
  				point3d.setElementValue(1, j*dy - offsetY );
  				point3d.setElementValue(2, -dz/2-offsetZ);
  				point2d = SimpleOperators.multiply(mat, point3d);
  
  				dsx0 = Math.abs(point3d.getElement(0) - cx);
  				dsy0 = Math.abs(point3d.getElement(1) - cy);
  				dsxy0_sqr = dsx0*dsx0 + dsy0*dsy0;
  
  				if (Debug1){
  					System.out.println("dsx0="+ dsx0 + ", dsy0=" + dsy0 + ", dsxy0_sqr= " + dsxy0_sqr);
  				}
  
  				if ( dsy0 > dsx0){
  					ds0 = dsy0;
  					point2dMinus = SimpleOperators.subtract(point2d, halfVoxelMoveX);
  					point2dPlus = SimpleOperators.add(point2d, halfVoxelMoveX);					
  				}else{
  					ds0 = dsx0;
  					point2dMinus = SimpleOperators.subtract(point2d, halfVoxelMoveY);
  					point2dPlus = SimpleOperators.add(point2d, halfVoxelMoveY);
  				}
  
  				coordLeft =  point2dMinus.getElement(0) / point2dMinus.getElement(2) + 0.5;
  				coordRight =  point2dPlus.getElement(0) / point2dPlus.getElement(2) + 0.5;
  
  				if ( coordLeft > coordRight ){
  					double temp = coordRight;
  					coordRight = coordLeft;
  					coordLeft = temp;	
  				}
  
  				coordLeft = Math.max(coordLeft, 0);
  				coordRight = Math.min(coordRight, maxU);
  
  				if (coordLeft >= coordRight){
  					continue;
  				}
  
  				iumin = (int)coordLeft;
  				iumax = (int)coordRight;
  				iumax = Math.min(iumax, maxU-1);
  
  				iumax = Math.min(iumax, iumin + MAX_WEIGHT_LENGTH - 1);
  
  				double[] weightU = new double[]{1,1,1,1,1,1};
  
  				if ( iumax == iumin){
  					weightU[0] = coordRight - coordLeft;
  				}else{
  					weightU[0] = iumin + 1 - coordLeft;
  					weightU[iumax - iumin] = coordRight - iumax;
  				};
  
  				coordBottom = point2d.getElement(1) / point2d.getElement(2) + 0.5;
  				point2d.add(fullVoxelMoveZ);
  				coordTop = point2d.getElement(1) / point2d.getElement(2) + 0.5;
  				coordStep = coordTop - coordBottom;
  
  
  				for ( int k = 0; k < maxK; k++ ){
  
  					double dsz0 = Math.abs(k*dz-offsetZ-cz); 
  
  					if ( coordBottom <= 0 ){
  						coordBottom = coordTop;
  						coordTop = coordTop + coordStep;
  						continue;
  					}
  
  					if ( coordTop >= maxV ){
  						break;
  					}
  
  					ivmin = (int)coordBottom;
  					ivmax = (int)coordTop;
  
  					double[] weightV = new double[]{1,1,1,1,1,1};
  
  					//ivmax = Math.min(ivmax, ivmin + MAX_LENGTH_V - 1);			
  					if ( ivmin == ivmax ){
  						weightV[0] = coordTop - coordBottom;
  					}else{
  						weightV[0] = ivmin + 1 - coordBottom;
  						weightV[ivmax - ivmin] = coordTop - ivmax;
  					};
  
  					if (Debug1){
  						printSimpleVector(point2d);
  						System.out.println( "Coords: coordLeft= " + coordLeft + ", \t coordRight=" + coordRight);
  						System.out.println( "Coords: coordBottom= " + coordBottom + ", \t coordTop=" + coordTop);
  						//System.out.println( "iu: iumin= " + iumin + ", \t iumax=" + iumax);
  						System.out.println("weights along u : " + weightU[0] + "\t" + weightU[1] + "\t" + weightU[2] + "\t" + weightU[3] + "\t" + weightU[4] + "\t" + weightU[5] );
  						System.out.println("weights along v : " + weightV[0] + "\t" + weightV[1] + "\t" + weightV[2] + "\t" + weightV[3] + "\t" + weightV[4] + "\t" + weightV[5] );
  					}
  
  					amplitude = Math.sqrt(dsz0*dsz0 + dsxy0_sqr) / ds0;
  
  					double tempVal = 0;
  
  					for ( int iu = iumin, iiu = 0; iu <= iumax; iu++, iiu++ ){
  						double sum = 0;
  						for (int iv = ivmin, iiv = 0; iv <= ivmax; iv++, iiv++){
  							sum += weightV[iiv] * projImage.getAtIndex(ip, iu, iv);
  						} //iv
  						tempVal += sum * weightU[iiu];
  					} //iu
  
  					tempVal = tempVal * amplitude;
  					volImage.addAtIndex(i, j, k, (float) tempVal);
  
  					//update for next voxel
  					coordBottom = coordTop;
  					coordTop = coordTop + coordStep;
  
  				} //k
  			} //j
  		} //i
  
  	}
  
  	@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;
  	}
  
  
  }
  /*
   * Copyright (C) 2010-2014 Meng Wu
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
  */