content: package edu.stanford.rsl.conrad.cuda;
  
  //TODO: Use our own matrices instead of Jama.Matrix
  
  import ij.ImagePlus;
  import ij.ImageStack;
  import ij.process.FloatProcessor;
  import ij.process.ImageProcessor;
  
  import jcuda.Pointer;
  import jcuda.Sizeof;
  import jcuda.driver.CUDA_ARRAY3D_DESCRIPTOR;
  import jcuda.driver.CUDA_MEMCPY3D;
  import jcuda.driver.CUaddress_mode;
  import jcuda.driver.CUarray;
  import jcuda.driver.CUarray_format;
  import jcuda.driver.CUcontext;
  import jcuda.driver.CUdevice;
  import jcuda.driver.CUdeviceptr;
  import jcuda.driver.CUdevprop;
  import jcuda.driver.CUfilter_mode;
  import jcuda.driver.CUfunction;
  import jcuda.driver.CUmemorytype;
  import jcuda.driver.CUmodule;
  import jcuda.driver.CUtexref;
  import jcuda.driver.JCudaDriver;
  import jcuda.runtime.dim3;
  
  import edu.stanford.rsl.apps.gui.Citeable;
  import edu.stanford.rsl.apps.gui.GUIConfigurable;
  import edu.stanford.rsl.conrad.geometry.Projection;
  import edu.stanford.rsl.conrad.geometry.trajectories.ProjectionTableFileTrajectory;
  import edu.stanford.rsl.conrad.geometry.trajectories.Trajectory;
  import edu.stanford.rsl.conrad.numerics.SimpleMatrix;
  import edu.stanford.rsl.conrad.numerics.SimpleVector;
  import edu.stanford.rsl.conrad.utils.CONRAD;
  import edu.stanford.rsl.conrad.utils.Configuration;
  import edu.stanford.rsl.conrad.utils.ImageUtil;
  
  
  
  /**
   * Forward projection expects input of a volumetric phantom scaled to mass density. Projection result {@latex.inline $p(\\vec{x})$} is then the accumulated mass along the ray {@latex.inline $\\vec{x}$} which consists of the line segments {@latex.inline $x_i$} in {@latex.inline $[\\textnormal{cm}]$} with the mass densities {@latex.inline $\\mu_i$} in {@latex.inline $[\\frac{\\textnormal{g}}{\\textnormal{cm}^3}]$}.
   * The actual projection is then computed as:<br>
   * {@latex.inline $$p(\\vec{x}) = \\sum_{i} x_i \\cdot \\mu_i$$}<BR>
   * The projection values are then returned in {@latex.inline $[\\frac{\\textnormal{g}}{\\textnormal{cm}^2}]$}
   * @author akmaier
   *
   */
  public class CUDAForwardProjector implements GUIConfigurable, Citeable {
  
  	static int bpBlockSize[] = {16, 16};
  
  	private static boolean debug = true;
  	/**
  	 * The CUDA module containing the kernel
  	 */
  	private CUmodule module = null;
  	/**
  	 * The handle for the CUDA function of the kernel that is to be called
  	 */
  	private CUfunction function = null;
  	/**
  	 * The 3D volume texture reference
  	 */
  	private CUtexref gTex3D = null;
  	private CUarray gVolume = null;
  
  	private CUdeviceptr gVolumeEdgeMaxPoint = null;
  	private CUdeviceptr gVolumeEdgeMinPoint = null;
  	private CUdeviceptr gInvARmatrix = null;
  	private CUdeviceptr gSrcPoint = null;
  	private CUdeviceptr gProjection = null;
  	/**
  	 * the context
  	 */
  	private CUcontext cuCtx = null;
  
  	private boolean initialized = false;
  	private boolean largeVolumeMode = false;
  	private int nSteps;
  	private int currentStep = 0;
  	private float [] voxelSize = null;
  	private float [] volumeSize = null;
  	private boolean configured = false;
  	private float [] volumeEdgeMinPoint = null;
  	private float [] volumeEdgeMaxPoint = null;
  	// buffer for 3D volume:
  	private float [] h_volume;
  	private int subVolumeZ;
  
  	private ImagePlus tex3D = null;
  
  	/**
  	 * Gets the volume to project
  	 * @return the volume as image plus
  	 */
  	public ImagePlus getTex3D() {
  		return tex3D;
  	}
  
  	/**
  	 * Sets the volume to project
  	 * @param tex3d
  	 */
  	public void setTex3D(ImagePlus tex3d) {
  		tex3D = tex3d;
  	}
  
  	private float [] projection;
  	private int width;
  	private int height;
  	private Trajectory geometry;
  
  	/**
  	 * 
  	 * 	
  	 *  Method: computeCanonicalProjectionMatrix<br>
  	 *  Author: Sungwon Yoon<br>
  	 *  Description:<br>
  	 *  <pre> 
  	 *         W -> W projection matrix = [ AR   t ]
  	 *         C -> C projection matrix = T0 * [AR  t] * T4
  	 *                 
  	 *                                 [ [ du(0)  dv(0) ]^-1   -0.5 ]
  	 *                where      T0 =  [ [ du(1)  dv(1) ]      -0.5 ]  
  	 *                                 [     0      0            1  ]      ,
  	 *
  	 *                                 [  dx   0    0   -(L-1)/2*dx  ]
  	 *                           T4 =  [  0   dy    0   -(M-1)/2*dy  ]
  	 *                                 [  0    0   dz   -(N-1)/2*dz  ]
  	 *                                 [  0    0    0         1      ]
  	 *  </pre>
  	 *
  	 *  C -> C projection matrix can be written as
  	 *  <pre>
  	 *        C -> C projection matrix = [ T0 * AR * T4(1:3,1:3)    T0 * (AR * T4(1:3,4) + t) ] 
  	 *  
  	 *  Therefore, the new invARmatrix = T4(1:3,1:3)^-1 * (AR)^-1 * T0^-1
  	 *                                   [ 1/dx    0      0   ]             [du(0)  dv(0)  0]
  	 *                                 = [  0     1/dy    0   ] * (AR)^-1 * [du(1)  dv(1)  0]
  	 *                                   [  0      0     1/dz ]             [  0      0    1]
  	 *
  	 *            and the new srcPoint = -T4(1:3,1:3)^-1 * T4(1:3,4) - T4(1:3,1:3)^-1 * (AR)^-1 * t
  	 *                                     [[ -0.5 * (L-1) ]     [ 1/dx    0      0   ]                 ]
  	 *                                 = - [[ -0.5 * (M-1) ]  +  [  0     1/dy    0   ] * srcPoint^{W}  ]
  	 *                                     [[ -0.5 * (N-1) ]     [  0      0     1/dz ]                 ]
  	 *  </pre>  
  	 *  
  	 *  <BR><BR>
  	 *  This implementation is consistent with Andreas Keil's Projection class and Benni's inversion.
  	 *                                   
  	 * @param canonicalProjMatrix is filled with a 3x4 projection matrix in this canonical format
  	 * @param invARmatrix is filled with the inverse of AR in canonical format
  	 * @param srcPoint is filled with the 3x1 source point in canonical format 
  	 * @param projectionMatrix the Matrix on which the conversion is based.
  	 */
  
  	public void computeCanonicalProjectionMatrix(float [] canonicalProjMatrix, float [] invARmatrix, float [] srcPoint, Jama.Matrix projectionMatrix){
  
  		double [] du = {Configuration.getGlobalConfiguration().getGeometry().getPixelDimensionX(), 0};
  		double [] dv = {0, Configuration.getGlobalConfiguration().getGeometry().getPixelDimensionY()};
  
  
  		// Assumed the vectors have only one non-zero element
  		// use 
  		du[0] = 1;
  		dv[1] = 1;
  
  
  
  		/*  ---------------------------------------------------
  		 *
  		 *  Canonical projection matrix used for BP computation
  		 *
  		 *  ---------------------------------------------------
  		 */
  		// T0 matrix
  		Jama.Matrix T0 = new Jama.Matrix (3,3);
  		double denom = 1.0f / (du[0]*dv[1] - dv[0]*du[1]);
  		T0.set(0,0, denom * dv[1]);
  		T0.set(0,1,-denom * dv[0]);
  		T0.set(1,0, -denom * du[1]);
  		T0.set(1,1, denom * du[0]);
  		T0.set(0,2, -0.5);
  		T0.set(1,2, -0.5);
  		T0.set(2,2, 1.0);
  		//T0.print(NumberFormat.getInstance(), 8);
  
  		// T4 matrix 
  		Jama.Matrix T4 = new Jama.Matrix(4,4);
  		T4.set(0,0, voxelSize[0]);
  		T4.set(1,1, voxelSize[1]);
  		T4.set(2,2, voxelSize[2]);
  		for (int k=0; k<3; ++k){
  			T4.set(k,3, -0.5 * (volumeSize[k] - 1.0) * voxelSize[k]);
  		}
  		T4.set(3,3, 1.0);
  
  		// New projection matrix in Canonical coord sys
  		Jama.Matrix tmpMatrix;
  		Jama.Matrix newProjMat;
  		tmpMatrix = projectionMatrix.times(T4);
  		newProjMat = T0.times(tmpMatrix);
  		for (int r=0; r<3; ++r) {
  			for (int c=0; c<4; ++c) {
  				// 3x3 matrix 1st row (indices 0-3), 2nd row (indices 4-7),
  				//   and 3rd row (indices 8-11)
  				canonicalProjMatrix[4*r + c] = (float) newProjMat.get(r,c);
  			}
  		}
  
  		/*  -----------------------------------------------------------
  		 *
  		 *  Canonical inverse projection matrix used for FP computation
  		 *
  		 *  -----------------------------------------------------------
  		 */
  		// Inverse of T0 matrix
  		//T0.inverse().print(NumberFormat.getInstance(), 8);
  		T0.set(0,0, du[0]);
  		T0.set(0,1, dv[0]);
  		T0.set(1,0, du[1]);
  		T0.set(1,1, dv[1]);
  		T0.set(0,2, 0);//0.5 * (du[0]+dv[0]));
  		T0.set(1,2, 0);//0.5 * (du[1]+dv[1]));
  		//T0.print(NumberFormat.getInstance(), 8);
  
  		// Inverse scaling by dx, dy, dz
  		Jama.Matrix invVoxelScale = new Jama.Matrix(3,3);
  		invVoxelScale.set(0,0, 1.0/voxelSize[0]);
  		invVoxelScale.set(1,1, 1.0/voxelSize[1]);
  		invVoxelScale.set(2,2, 1.0/voxelSize[2]);
  
  		// New invARmatrix in the Canonical coord sys
  		Jama.Matrix invARmatrixMat = projectionMatrix.getMatrix(0, 2, 0, 2).inverse();
  		tmpMatrix = invARmatrixMat.times(T0);    // invARmatrix_ * T0^{-1}
  		Jama.Matrix invAR = invVoxelScale.times(tmpMatrix);     // invVoxelScale * (invARmatrix_ * T0)
  		for (int r=0; r<3; ++r) {
  			for (int c=0; c<3; ++c) {
  				// 3x3 matrix 1st row (indices 0-2), 2nd row (indices 3-5),
  				//   and 3rd row (indices 6-8)
  				invARmatrix[3*r + c] = (float) invAR.get(r,c);
  
  			}
  		}
  		//invAR.print(NumberFormat.getInstance(), 6);
  		// New srcPoint in the Canonical coord sys
  
  		Jama.Matrix srcPtW = computeSrcPt(projectionMatrix, invARmatrixMat);
  		srcPoint[0] = (float) -(-0.5 * (volumeSize[0] -1.0) + invVoxelScale.get(0,0) * srcPtW.get(0, 0)); 
  		srcPoint[1] = (float) -(-0.5 * (volumeSize[1] -1.0) + invVoxelScale.get(1,1) * srcPtW.get(1, 0)); 
  		srcPoint[2] = (float) -(-0.5 * (volumeSize[2] -1.0) + invVoxelScale.get(2,2) * srcPtW.get(2, 0)); 
  
  
  	}
  
  	/**
  	 * computes the location of the Source Point given a 3x4 projection matrix and and inverted 3x3 AR Projection matrix.
  	 * Used in computeCanonicalProjectionMatrix
  	 * 
  	 * @param projectionMatrix the original projection matrix
  	 * @param invertedProjMatrix the inverted AR projection matrix
  	 * @return the source point
  	 */
  	private Jama.Matrix computeSrcPt(Jama.Matrix projectionMatrix, Jama.Matrix invertedProjMatrix) {
  		Jama.Matrix at = projectionMatrix.getMatrix(0, 2, 3, 3);
  		at.times(-1.0);
  		return invertedProjMatrix.times(at);
  	}
  
  	/**
  	 * Initiates communication with the CUDA card.
  	 */
  	private void init(){
  		if (!initialized) {
  			largeVolumeMode = false;
  			// Initialize the JCudaDriver. Note that this has to be done from 
  			// the same thread that will later use the JCudaDriver API.
  			JCudaDriver.setExceptionsEnabled(true);
  			JCudaDriver.cuInit(0);
  			CUdevice dev = new CUdevice();
  			JCudaDriver.cuDeviceGet(dev, 0);
  			cuCtx = new CUcontext();
  			JCudaDriver.cuCtxCreate(cuCtx, 0, dev);
  			// check space on device:
  			int [] memory = new int [1];
  			int [] total = new int [1]; 
  			JCudaDriver.cuDeviceTotalMem(memory, dev);
  			JCudaDriver.cuMemGetInfo(memory, total);
  			int availableMemory = (int) (CUDAUtil.correctMemoryValue(memory[0]) / ((long) 1024 * 1024));
  			int requiredMemory = (int)(((
  					((double) volumeSize[0]) * volumeSize[1] * ((double) volumeSize[2]) * Sizeof.FLOAT) 
  					+ (((double) height) * width * Sizeof.FLOAT)) 
  					/ (1024.0 * 1024));
  			if (debug) {
  				System.out.println("Total available Memory on CUDA card:" + availableMemory);
  				System.out.println("Required Memory on CUDA card:" + requiredMemory);
  			}
  			if (requiredMemory > availableMemory){
  				nSteps = CUDAUtil.iDivUp (requiredMemory, (int)(availableMemory));
  				if (debug) System.out.println("Switching to large volume mode with nSteps = " + nSteps);
  				largeVolumeMode = true;
  			}
  			if (debug) {
  				CUdevprop prop = new CUdevprop();
  				JCudaDriver.cuDeviceGetProperties(prop, dev);
  				System.out.println(prop.toFormattedString());
  			}
  			// Load the CUBIN file containing the kernel
  			module = new CUmodule();
  			JCudaDriver.cuModuleLoad(module, "projectWithCuda.sm_10.cubin");
  
  			// Obtain a function pointer to the kernel function. This function
  			// will later be called. 
  			// 
  			function = new CUfunction();
  			JCudaDriver.cuModuleGetFunction(function, module, "_Z13projectKernelPfjf");
  
  
  			int memorysize = (int) (volumeSize[0] * volumeSize[1] * volumeSize[2] * Sizeof.FLOAT);
  			if (largeVolumeMode){
  				subVolumeZ = CUDAUtil.iDivUp((int) volumeSize[2], nSteps);
  				if(debug) System.out.println("SubVolumeZ: " + subVolumeZ);
  				h_volume = new float[(int) (volumeSize[0] * volumeSize[1] * subVolumeZ)];
  				memorysize = (int) (volumeSize[0] * volumeSize[1] * subVolumeZ * Sizeof.FLOAT);
  				if(debug)System.out.println("Memory: " + memorysize);
  			} else {
  				h_volume = new float[(int) (volumeSize[0] * volumeSize[1] * volumeSize[2])];
  				subVolumeZ = (int) volumeSize[2];
  				nSteps = 1;
  			}
  
  			copyVolumeToCard();
  
  			gTex3D = new CUtexref();
  			// Obtain the 3D texture reference for the volume data from 
  			// the module, set its parameters and assign the 3D volume 
  			// data array as its reference.
  			JCudaDriver.cuModuleGetTexRef(gTex3D, module, "gTex3D");
  			JCudaDriver.cuTexRefSetFilterMode(gTex3D,
  					CUfilter_mode.CU_TR_FILTER_MODE_LINEAR);
  			JCudaDriver.cuTexRefSetAddressMode(gTex3D, 0,
  					CUaddress_mode.CU_TR_ADDRESS_MODE_CLAMP);
  			JCudaDriver.cuTexRefSetAddressMode(gTex3D, 1,
  					CUaddress_mode.CU_TR_ADDRESS_MODE_CLAMP);
  			JCudaDriver.cuTexRefSetFormat(gTex3D,
  					CUarray_format.CU_AD_FORMAT_FLOAT, 1);
  			JCudaDriver.cuTexRefSetFlags(gTex3D,
  					JCudaDriver.CU_TRSF_READ_AS_INTEGER);
  			JCudaDriver.cuTexRefSetArray(gTex3D, gVolume,
  					JCudaDriver.CU_TRSA_OVERRIDE_FORMAT);
  
  
  			CUDAUtil.copyFloatArrayToDevice(voxelSize, module, "gVoxelElementSize");
  
  			// copy volume to device
  			gProjection = new CUdeviceptr();
  			JCudaDriver.cuMemAlloc(gProjection, width * height * Sizeof.FLOAT);
  
  
  			initialized = true;
  		}
  
  	}
  
  
  	/**
  	 * Load the inverted projection matrix for the current projection and resets the projection data.
  	 * @param projectionNumber
  	 */
  	private void prepareProjection(int projectionNumber){
  		float [] cann = new float[3*4];
  		float [] invAR = new float[3*3];
  		float [] srcP = new float[3];
  		SimpleMatrix projMat = geometry.getProjectionMatrix(projectionNumber).computeP();
  		double [][] mat = new double [3][4];
  		projMat.copyTo(mat);
  		computeCanonicalProjectionMatrix(cann, invAR, srcP, new Jama.Matrix(mat));
  		if (gInvARmatrix == null){
  			gInvARmatrix = CUDAUtil.copyFloatArrayToDevice(invAR, module, "gInvARmatrix");
  			gSrcPoint = CUDAUtil.copyFloatArrayToDevice(srcP, module, "gSrcPoint");
  		} else {
  			CUDAUtil.updateFloatArrayOnDevice(gInvARmatrix, invAR, module);
  			CUDAUtil.updateFloatArrayOnDevice(gSrcPoint, srcP, module);
  		}
  		// reset Projection Data
  		JCudaDriver.cuMemsetD32(gProjection, 0, width * height);
  	}
  
  
  
  	/**
  	 * Main method checks different computation methods and compares to a reference implementation.
  	 * The P matrix
  	 * [[   -1.957647643397161     -1.5159990314154403    2.6654943287669416E-5    284.3151786700395 ]; 
  	 *  [ -0.22689199744970734     0.21980776263540414      -2.4371669460715766   219.19724835256451 ];
  	 *  [-9.088284790690905E-4    8.644729632667306E-4    1.7220251080057994E-5                  1.0 ]]
  	 * 
  	 * should yield the following source point in texture coordinates
  	 * 1402.969301
  	 * -851.170185
  	 * 228.742784
  	 * 
  	 * and this inverse (up to a scalar factor):
  	 * -0.563161 	-0.006972 	-985.825785
  	 * -0.592038 	 0.008988 	1273.025240
  	 * -0.000967 	-0.819165 	206.591022
  	 * 
  	 * This part should be moved to a test-case some time later.
  	 * 
  	 * @param args
  	 */
  	public static void main(String[] args){
  		Configuration.loadConfiguration();
  		/*
  		double [] matrix = {-1.957647643397161,    -1.5159990314154403,    2.6654943287669416E-5,   284.3151786700395,
  				-0.22689199744970734,    0.21980776263540414,    -2.4371669460715766,   219.19724835256451,
  				-9.088284790690905E-4,    8.644729632667306E-4,    1.7220251080057994E-5,    1.0 };
  		*/
  		Projection proj = new Projection();
  		proj.setPMatrixSerialization("[[-1.957647643397161    -1.5159990314154403    2.6654943287669416E-5   284.3151786700395]; " +
  				"[-0.22689199744970734    0.21980776263540414    -2.4371669460715766   219.19724835256451 ];" +
  		"[ -9.088284790690905E-4    8.644729632667306E-4    1.7220251080057994E-5    1.0 ]]");
  		SimpleVector center = proj.computeCameraCenter();
  		System.out.println("center = "+ center);
  		System.out.println(proj.computeP());
  		System.out.println(proj.getRTKinv());
  		double scale = 206.591022 / proj.getRTKinv().getElement(2, 2);
  		System.out.println(proj.getRTKinv().multipliedBy(scale) +  "\n" + scale + " " + proj.computeSourceToDetectorDistance(new SimpleVector(0.320, 0.320))[0]);
  
  		SimpleMatrix m = proj.computeP().getSubMatrix(0, 0, 3, 3);
  		System.out.println("M="+ m.inverse(SimpleMatrix.InversionType.INVERT_QR).multipliedBy(2));
  
  		float [] cann = new float[3*4];
  		float [] invAR = new float[3*3];
  		float [] srcP = new float[3];
  		SimpleMatrix projMat = proj.computeP();
  		double [][] mat = new double [3][4];
  		projMat.copyTo(mat);
  		CUDAForwardProjector cudaForwardProjector = new CUDAForwardProjector();
  		cudaForwardProjector.voxelSize = new float [] {0.5f, 0.5f, 0.5f};
  		cudaForwardProjector.volumeSize = new float [] {512f, 512f, 512f};
  		cudaForwardProjector.computeCanonicalProjectionMatrix(cann, invAR, srcP, new Jama.Matrix(mat));
  		for (int i =0; i < 9; i++){
  			System.out.println(invAR[i]);
  		}
  		for (int i =0; i < 3; i++){
  			System.out.println(srcP[i]);
  		}
  	}
  
  	/**
  	 * destroys the CUDA context and frees the allocated memory
  	 */
  	private void unload(){
  		JCudaDriver.cuMemFree(gProjection);
  		// destory context
  		JCudaDriver.cuCtxDestroy(cuCtx);
  	}
  
  
  	/**
  	 * loads the actual CUDA kernel and performs the projection
  	 * @param projectionNumber the projection number.
  	 * @return
  	 */
  	private ImageProcessor project(int projectionNumber){
  		init();
  		prepareProjection(projectionNumber);
  
  		Pointer dOut = Pointer.to(gProjection);
  		Pointer pStride = Pointer.to(new int[]{width});
  		Pointer pstepsize = Pointer.to(new float[]{(float) 1});
  
  		int offset = 0;
  
  		offset = CUDAUtil.align(offset, Sizeof.POINTER);
  		JCudaDriver.cuParamSetv(function, offset, dOut, Sizeof.POINTER);
  		offset += Sizeof.POINTER;
  
  		offset = CUDAUtil.align(offset, Sizeof.INT);
  		JCudaDriver.cuParamSetv(function, offset, pStride, Sizeof.INT);
  		offset += Sizeof.INT;
  
  		offset = CUDAUtil.align(offset, Sizeof.FLOAT);
  		JCudaDriver.cuParamSetv(function, offset, pstepsize, Sizeof.FLOAT);
  		offset += Sizeof.FLOAT;
  
  		JCudaDriver.cuParamSetSize(function, offset);
  
  		dim3 gridSize = new dim3(
  				CUDAUtil.iDivUp(width, bpBlockSize[0]), 
  				CUDAUtil.iDivUp(height, bpBlockSize[0]), 
  				1);
  
  
  		//System.out.println("Grid: " + gridSize);
  
  		JCudaDriver.cuFuncSetBlockShape(function, bpBlockSize[0], bpBlockSize[1], 1);
  		JCudaDriver.cuLaunchGrid(function, gridSize.x, gridSize.y);
  		JCudaDriver.cuCtxSynchronize();
  
  		JCudaDriver.cuMemcpyDtoH(Pointer.to(projection), gProjection, width * height * Sizeof.FLOAT);
  		FloatProcessor fl = new FloatProcessor(width, height, projection, null);
  		
  
  		// TODO: Normalization is never considered in the backprojectors, 
  		// 		 thus, iteratively applying forward and backward projections
  		//		 would yield to a scaling issue!
  		// conversion from [g*mm/cm^3] = [g*0.1cm/cm^3] to [g/cm^2]
  		// fl.multiply(1.0 / 10);
  		
  		if (geometry instanceof ProjectionTableFileTrajectory){
  			fl.flipVertical();
  		}
  		return fl;
  	}
  
  	/**
  	 * Starts projection and returns Projection Data, as ImagePlus
  	 * @return the projection stack
  	 */
  	public ImagePlus project(){
  		ImageStack stack = new ImageStack(width, height);
  		for (int i = 0; i < geometry.getNumProjectionMatrices(); i++){
  			stack.addSlice("Projection " + i, project(i).duplicate());
  		}
  		if (largeVolumeMode){
  			// play it again, Sam!
  			while (currentStep < nSteps) {
  				currentStep++;
  
  				if (currentStep == nSteps) break;
  				if (debug) System.out.println("Processing step " + currentStep + " of " + nSteps);
  				copyVolumeToCard();
  				for (int i = 0; i < geometry.getNumProjectionMatrices(); i++){
  					ImageUtil.addProcessors(stack.getProcessor(i+1), project(i));
  					//stack.addSlice("slice " + i, project(i).duplicate());
  				}
  			}
  		}
  		ImagePlus image = new ImagePlus();
  		image.setStack("Forward Projection of " + tex3D.getTitle(), stack);
  		unload();
  		//image.show();
  		return image;
  	}
  
  	private void copyVolumeToCard(){
  
  		// write the volume from ImagePlus to a float array.
  		for (int k = 0; k < subVolumeZ; k++){
  			int index = ((nSteps - currentStep -1)*subVolumeZ) + k + 1;
  			if (index <= tex3D.getStackSize()){
  				ImageProcessor currentSlice = tex3D.getStack().getProcessor(index);
  				for (int i = 0; i< volumeSize[0]; i++){
  					for(int j =0; j< volumeSize[1]; j++){
  						boolean flip =  false;
  
  						int index2 = (int) (
  								(
  										(
  												(((int)volumeSize[1]) * (k))
  												+ ((int)j)
  										) * ((int)volumeSize[0])
  								)
  								+ ((int)i));
  						if (flip) {
  							index2 = (int) (
  									// Note that we must flip the coordinates of the volume slices in order to be compatible with the CUDA Backprojector...
  									(
  											(
  													(((int)volumeSize[1]) * (subVolumeZ - k -1))
  													+ ((int)volumeSize[1]-j-1)
  											) * ((int)volumeSize[0])
  									)
  									+ ((int)volumeSize[0] - i -1));
  						}
  						if (index2 >= (((int)volumeSize[0]) * ((int)volumeSize[1])* ((int)subVolumeZ))){
  							System.out.println(k + " " + i + " " + j);
  							break;
  						}
  						h_volume[index2] = currentSlice.getPixelValue(i, j);
  					}
  				}
  			} else {
  				System.out.println("Not in Volume: " + index + " " + nSteps + " " + currentStep);
  				for (int i = 0; i< volumeSize[0]; i++){
  					for(int j =0; j< volumeSize[1]; j++){
  						h_volume[(int)((((((int)volumeSize[1]) * k) + j) * ((int)volumeSize[0]))+i)] = 0;
  					}
  				}
  			}
  		}
  		int test = currentStep;
  		volumeEdgeMaxPoint[2] = (float) ((test * subVolumeZ) + subVolumeZ -0.5 - CONRAD.SMALL_VALUE);
  		volumeEdgeMinPoint[2] = (float) ((test * subVolumeZ) -0.5 - CONRAD.SMALL_VALUE);
  		if (debug) System.out.println("New volume z min: " + volumeEdgeMinPoint[2] + " new volume z max: " + volumeEdgeMaxPoint[2]);
  		if (gVolumeEdgeMaxPoint == null){
  			gVolumeEdgeMaxPoint = CUDAUtil.copyFloatArrayToDevice(volumeEdgeMaxPoint, module, "gVolumeEdgeMaxPoint");
  			gVolumeEdgeMinPoint = CUDAUtil.copyFloatArrayToDevice(volumeEdgeMinPoint, module, "gVolumeEdgeMinPoint");
  		} else {
  			CUDAUtil.updateFloatArrayOnDevice(gVolumeEdgeMaxPoint, volumeEdgeMaxPoint, module);
  			CUDAUtil.updateFloatArrayOnDevice(gVolumeEdgeMinPoint, volumeEdgeMinPoint, module);
  		}  
  
  		if (gVolume == null) {
  			gVolume = new CUarray();
  			// Create the 3D array that will contain the volume data
  			// and will be accessed via the 3D texture
  			CUDA_ARRAY3D_DESCRIPTOR allocateArray = new CUDA_ARRAY3D_DESCRIPTOR();
  			allocateArray.Width = (int) volumeSize[0];
  			allocateArray.Height = (int) volumeSize[1];
  			allocateArray.Depth = subVolumeZ;
  			allocateArray.Format = CUarray_format.CU_AD_FORMAT_FLOAT;
  			allocateArray.NumChannels = 1;
  			JCudaDriver.cuArray3DCreate(gVolume, allocateArray);
  		}
  
  		// Copy the volume data data to the 3D array
  		CUDA_MEMCPY3D copy = new CUDA_MEMCPY3D();
  		copy.srcMemoryType = CUmemorytype.CU_MEMORYTYPE_HOST;
  		copy.srcHost = Pointer.to(h_volume);
  		copy.srcPitch = (int) volumeSize[0] * Sizeof.FLOAT;
  		copy.srcHeight = (int) volumeSize[1];
  
  		copy.dstMemoryType = CUmemorytype.CU_MEMORYTYPE_ARRAY;
  		copy.dstArray = gVolume;
  		copy.dstPitch = (int) volumeSize[0] * Sizeof.FLOAT;
  		copy.dstHeight = (int) volumeSize[1];
  		copy.WidthInBytes = (int) volumeSize[0] * Sizeof.FLOAT;
  		copy.Height = (int) volumeSize[1];
  		copy.Depth = subVolumeZ;
  		JCudaDriver.cuMemcpy3D(copy);
  	}
  
  	/**
  	 * Start GUI configuration. Reads from global Configuration.
  	 */
  	public void configure() throws Exception {
  		// TODO Auto-generated method stub
  		voxelSize = new float [3];
  		volumeSize = new float [3];
  		Configuration config = Configuration.getGlobalConfiguration();
  		voxelSize[0] = (float) config.getGeometry().getVoxelSpacingX();
  		voxelSize[1] = (float) config.getGeometry().getVoxelSpacingY();
  		voxelSize[2] = (float) config.getGeometry().getVoxelSpacingZ();
  		volumeSize[0] = config.getGeometry().getReconDimensionX();
  		volumeSize[1] = config.getGeometry().getReconDimensionY();
  		volumeSize[2] = config.getGeometry().getReconDimensionZ();
  		volumeEdgeMinPoint = new float[3];
  		for (int i=0; i < 3; i ++){
  			volumeEdgeMinPoint[i] = (float) (-0.5 + CONRAD.SMALL_VALUE);
  		}
  		volumeEdgeMaxPoint = new float[3];
  		for (int i=0; i < 3; i ++){
  			volumeEdgeMaxPoint[i] = (float) (volumeSize[i] -0.5 - CONRAD.SMALL_VALUE);
  		}
  		width = config.getGeometry().getDetectorWidth();
  		height = config.getGeometry().getDetectorHeight();
  		geometry = config.getGeometry();
  
  		if (debug) System.out.println("Projection Matrices: " + geometry.getNumProjectionMatrices());
  		projection = new float[width * height];
  		configured = true;
  	}
  
  	/**
  	 * returns whether the projector was already configured or not.
  	 */
  	public boolean isConfigured() {
  		return configured;
  	}
  
  	/**
  	 * Returns a reference to literature describing this algorithm in Bibtex format
  	 */
  	public String getBibtexCitation() {
  		String bibtex = "@ARTICLE{Galigekere03-CBR,\n" +
  		"  author = {{Galigekere}, R. R. and {Wiesent}, K. and {Holdsworth}, D. W.},\n" +
  		"  title = \"{{Cone-Beam Reprojection Using Projection-Matrices}}\",\n" +
  		"  journal = {{IEEE Transactions on Medical Imaging}},\n" +
  		"  year = 2003,\n" +
  		"  volume = 22,\n"+
  		"  number = 10,\n" +
  		"  pages = {1202-1214}\n" +
  		"}";
  		return bibtex;
  	}
  
  	/**
  	 * Returns a reference to literature describing this algorithm in Medline
  	 */
  	public String getMedlineCitation() {
  		return "Galigekere RR, Wiesent K, and Holdsworth DW. Cone-Beam Reprojection Using Projection-Matrices. IEEE Transactions on Medical Imaging 22(10):1202-14 2003.";
  	}
  
  
  
  }
  
  /*
   * Copyright (C) 2010-2014 - Andreas Maier 
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
  */
  