content: package edu.stanford.rsl.conrad.cuda;
  
  import java.util.ArrayList;
  import java.util.Arrays;
  
  import edu.stanford.rsl.conrad.data.numeric.Grid2D;
  import edu.stanford.rsl.conrad.io.ImagePlusDataSink;
  import edu.stanford.rsl.conrad.numerics.SimpleMatrix;
  import edu.stanford.rsl.conrad.reconstruction.VOIBasedReconstructionFilter;
  import edu.stanford.rsl.conrad.utils.CONRAD;
  import edu.stanford.rsl.conrad.utils.Configuration;
  import edu.stanford.rsl.conrad.utils.ImageGridBuffer;
  import edu.stanford.rsl.conrad.utils.ImageUtil;
  import edu.stanford.rsl.conrad.utils.VisualizationUtil;
  import ij.IJ;
  import ij.ImagePlus;
  import jcuda.Pointer;
  import jcuda.Sizeof;
  import jcuda.driver.CUDA_ARRAY_DESCRIPTOR;
  import jcuda.driver.CUDA_MEMCPY2D;
  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.JCuda;
  import jcuda.runtime.dim3;
  /**
  This is the existing CUDABackProjector with added motion compensation. The motion compensated kernel is called and the respiration signal is added, other than that the class is the same.
  
  @author Marco Bögel
  
   */
  
  public class CUDACompensatedBackProjector extends VOIBasedReconstructionFilter implements Runnable {
  
  	double mot[] = new double[200];
  	/**
  	 * 
  	 */
  	private static final long serialVersionUID = 7732291211252379464L;
  
  	/**
  	 * The CUDA module containing the kernel
  	 */
  	private CUmodule module = null;
  	//private static Object lock = new Object();
  	private static boolean debug = true;
  
  	// Pre-determined kernel block size
  	static int bpBlockSize[] = {32, 16};
  
  	/**
  	 * The handle for the CUDA function of the kernel that is to be called
  	 */
  	private CUfunction function = null;
  
  	/**
  	 * The volume data that is to be reconstructed
  	 */
  	protected float h_volume[];
  
  	/**
  	 * The 2D projection texture reference
  	 */
  	private CUtexref projectionTex = null;
  
  	/**
  	 * The grid size of the kernel execution
  	 */
  	private dim3 gridSize = null;
  
  	/**
  	 * the context
  	 */
  	private CUcontext cuCtx = null;
  
  
  	/**
  	 * The global variable of the module which stores the
  	 * view matrix.
  	 */
  	private CUdeviceptr projectionMatrix = null;
  	private CUdeviceptr volStride = null;
  	private CUdeviceptr volumePointer = null;
  	private CUarray projectionArray = null;
  
  
  	protected ImageGridBuffer projections;
  	protected ArrayList<Integer> projectionsAvailable;
  	protected ArrayList<Integer> projectionsDone;
  	private boolean largeVolumeMode = false;
  	private int nSteps = 1;
  	private int subVolumeZ = 0;
  
  	private boolean initialized = false;
  
  	public CUDACompensatedBackProjector () {
  		super();
  	}
  
  	@Override
  	public void prepareForSerialization(){
  		super.prepareForSerialization();
  		projectionMatrix = null;
  		volStride = null;
  		volumePointer = null;
  		projectionArray = null;
  		projections = null;
  		projectionsAvailable =null;
  		projectionsDone = null;
  		cuCtx = null;
  		gridSize = null;
  		projectionTex = null;
  		h_volume = null;
  		initialized = false;
  		function = null;
  		module = null;
  	}
  
  	@Override
  	public void configure() throws Exception{
  		boolean success = true;
  		configured = success;
  	}
  
  	public void reset(){
  		projectionArray = null;
  		JCuda.cudaThreadExit();
  	}
  
  	protected void init(){
  		if (!initialized) {
  			largeVolumeMode = false;
  
  			int reconDimensionX = getGeometry().getReconDimensionX();
  			int reconDimensionY = getGeometry().getReconDimensionY();
  			int reconDimensionZ = getGeometry().getReconDimensionZ();
  			projections = new ImageGridBuffer();
  			projectionsAvailable = new ArrayList<Integer>();
  			projectionsDone = new ArrayList<Integer>();
  			// 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 = CUDAUtil.getBestDevice();
  			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) reconDimensionX) * reconDimensionY * ((double) reconDimensionZ) * Sizeof.FLOAT) 
  					+ (((double)Configuration.getGlobalConfiguration().getGeometry().getDetectorHeight()) * Configuration.getGlobalConfiguration().getGeometry().getDetectorWidth() * 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, "respirationCompensatedBackprojectWithCuda.ptx");
  
  			// Obtain a function pointer to the kernel function. This function
  			// will later be called. 
  			// 
  			function = new CUfunction();
  			JCudaDriver.cuModuleGetFunction(function, module,
  			"_Z17backprojectKernelPfiiiffffff");
  			// create the reconstruction volume;
  			int memorysize = reconDimensionX * reconDimensionY * reconDimensionZ * Sizeof.FLOAT;
  			if (largeVolumeMode){
  				subVolumeZ = CUDAUtil.iDivUp(reconDimensionZ, nSteps);
  				if(debug) System.out.println("SubVolumeZ: " + subVolumeZ);
  				h_volume = new float[reconDimensionX * reconDimensionY * subVolumeZ];
  				memorysize = reconDimensionX * reconDimensionY * subVolumeZ * Sizeof.FLOAT;
  				if(debug)System.out.println("Memory: " + memorysize);
  			} else {
  				h_volume = new float[reconDimensionX * reconDimensionY * reconDimensionZ];	
  			}
  			// copy volume to device
  			volumePointer = new CUdeviceptr();
  			JCudaDriver.cuMemAlloc(volumePointer, memorysize);
  			JCudaDriver.cuMemcpyHtoD(volumePointer, Pointer.to(h_volume), memorysize);
  
  			// compute adapted volume size 
  			//    volume size in x = multiple of bpBlockSize[0]
  			//    volume size in y = multiple of bpBlockSize[1]
  
  			int adaptedVolSize[] = new int[3];
  			if ((reconDimensionX % bpBlockSize[0] ) == 0){
  				adaptedVolSize[0] = reconDimensionX;
  			} else {
  				adaptedVolSize[0] = ((reconDimensionX / bpBlockSize[0]) + 1) * bpBlockSize[0];
  			}
  			if ((reconDimensionY % bpBlockSize[1] ) == 0){
  				adaptedVolSize[1] = reconDimensionY;
  			} else {
  				adaptedVolSize[1] = ((reconDimensionY / bpBlockSize[1]) + 1) * bpBlockSize[1];
  			}
  			adaptedVolSize[2] = reconDimensionZ;
  			int volStrideHost [] = new int[2];
  			// compute volstride and copy it to constant memory
  			volStrideHost[0] = adaptedVolSize[0];
  			volStrideHost[1] = adaptedVolSize[0] * adaptedVolSize[1];
  
  
  			volStride = new CUdeviceptr();
  			JCudaDriver.cuModuleGetGlobal(volStride, new int[1], module, "gVolStride");
  			JCudaDriver.cuMemcpyHtoD(volStride, Pointer.to(volStrideHost), Sizeof.INT * 2);
  
  			// Calculate new grid size
  			gridSize = new dim3(
  					CUDAUtil.iDivUp(adaptedVolSize[0], bpBlockSize[0]), 
  					CUDAUtil.iDivUp(adaptedVolSize[1], bpBlockSize[1]), 
  					adaptedVolSize[2]);
  
  
  			// Obtain the global pointer to the view matrix from
  			// the module
  			projectionMatrix = new CUdeviceptr();
  			JCudaDriver.cuModuleGetGlobal(projectionMatrix, new int[1], module, "gProjMatrix");
  
  			initialized = true;
  		}
  
  	}
  
  	private synchronized void unload(){
  		if (initialized) {
  
  			if (projectionArray != null) {
  				JCudaDriver.cuArrayDestroy(projectionArray);
  			}
  
  			int reconDimensionX = getGeometry().getReconDimensionX();
  			int reconDimensionY = getGeometry().getReconDimensionY();
  			int reconDimensionZ = getGeometry().getReconDimensionZ();
  
  
  			if ((projectionVolume != null) && (!largeVolumeMode)) {
  				// fetch data
  				int memorysize = reconDimensionX * reconDimensionY * reconDimensionZ * 4;
  				JCudaDriver.cuMemcpyDtoH(Pointer.to(h_volume), volumePointer, memorysize);
  				int width = projectionVolume.getSize()[0];
  				int height = projectionVolume.getSize()[1];
  				if (this.useVOImap) {
  					for (int k = 0; k < projectionVolume.getSize()[2]; k++){
  						for (int j = 0; j < height; j++){
  							for (int i = 0; i < width; i++){			
  								float value = h_volume[(((height * k) + j) * width) + i];
  								if (voiMap[i][j][k]) {
  									projectionVolume.setAtIndex(i, j, k, value);
  								} else {
  									projectionVolume.setAtIndex(i, j, k, 0);
  								}
  							}
  						}
  					}
  				} else {
  					for (int k = 0; k < projectionVolume.getSize()[2]; k++){
  						for (int j = 0; j < height; j++){
  							for (int i = 0; i < width; i++){			
  								float value = h_volume[(((height * k) + j) * width) + i];
  								projectionVolume.setAtIndex(i, j, k, value);
  							}
  						}
  					}
  				}
  			} else {
  				System.out.println("Check ProjectionVolume. It seems null.");
  			}
  
  
  			h_volume = null;
  			// free memory on device
  
  			JCudaDriver.cuMemFree(volumePointer);
  			// destory context
  			JCudaDriver.cuCtxDestroy(cuCtx);
  
  			reset();
  
  			initialized = false;
  		}
  	}
  
  	private synchronized void initProjectionMatrix(int projectionNumber){
  		// load projection Matrix for current Projection.
  		SimpleMatrix pMat = getGeometry().getProjectionMatrix(projectionNumber).computeP();
  		float [] pMatFloat = new float[pMat.getCols() * pMat.getRows()];
  		for (int j = 0; j< pMat.getRows(); j++) {
  			for (int i = 0; i< pMat.getCols(); i++) {
  
  				pMatFloat[(j * pMat.getCols()) + i] = (float) pMat.getElement(j, i);
  			}
  		}
  		JCudaDriver.cuMemcpyHtoD(projectionMatrix, Pointer.to(pMatFloat), Sizeof.FLOAT * pMatFloat.length);
  	}
  
  	private synchronized void initProjectionData(Grid2D projection){
  		initialize(projection);
  		if (projection != null){ 
  			float [] proj= new float[projection.getWidth() * projection.getHeight()];
  
  			for(int i = 0; i< projection.getWidth(); i++){
  				for (int j =0; j < projection.getHeight(); j++){
  					proj[(j*projection.getWidth()) + i] = projection.getPixelValue(i, j);
  				}
  			}
  
  			if (projectionArray == null) {
  				// Create the 2D array that will contain the
  				// projection data. 
  				projectionArray = new CUarray();
  				CUDA_ARRAY_DESCRIPTOR ad = new CUDA_ARRAY_DESCRIPTOR();
  				ad.Format = CUarray_format.CU_AD_FORMAT_FLOAT;
  				ad.Width = projection.getWidth();
  				ad.Height = projection.getHeight();
  				ad.NumChannels = 1;//projection.getNChannels();
  				JCudaDriver.cuArrayCreate(projectionArray, ad);
  			}
  
  			// Copy the projection data to the array  
  			CUDA_MEMCPY2D copy2 = new CUDA_MEMCPY2D();
  			copy2.srcMemoryType = CUmemorytype.CU_MEMORYTYPE_HOST;
  			copy2.srcHost = Pointer.to(proj);
  			copy2.srcPitch = projection.getWidth() * Sizeof.FLOAT;
  			copy2.dstMemoryType = CUmemorytype.CU_MEMORYTYPE_ARRAY;
  			copy2.dstArray = projectionArray;
  			copy2.WidthInBytes = projection.getWidth() * Sizeof.FLOAT;
  			copy2.Height = projection.getHeight();
  			JCudaDriver.cuMemcpy2D(copy2);
  
  			// Obtain the texture reference from the module, 
  			// set its parameters and assign the projection  
  			// array as its reference.
  			projectionTex = new CUtexref();
  			JCudaDriver.cuModuleGetTexRef(projectionTex, module, "gTex2D");
  			JCudaDriver.cuTexRefSetFilterMode(projectionTex,
  					CUfilter_mode.CU_TR_FILTER_MODE_LINEAR);
  			JCudaDriver.cuTexRefSetAddressMode(projectionTex, 0,
  					CUaddress_mode.CU_TR_ADDRESS_MODE_CLAMP);
  			JCudaDriver.cuTexRefSetFlags(projectionTex,
  					JCudaDriver.CU_TRSF_READ_AS_INTEGER);
  			JCudaDriver.cuTexRefSetFormat(projectionTex,
  					CUarray_format.CU_AD_FORMAT_FLOAT, 4);
  			JCudaDriver.cuTexRefSetArray(projectionTex, projectionArray,
  					JCudaDriver.CU_TRSA_OVERRIDE_FORMAT);
  
  			// Set the texture references as parameters for the function call
  			JCudaDriver.cuParamSetTexRef(function, JCudaDriver.CU_PARAM_TR_DEFAULT,
  					projectionTex);
  		} else {
  			System.out.println("Projection was null!!");
  		}
  	}
  
  	@Override
  	public String getName() {
  		return "CUDA Backprojector";
  	}
  
  
  	@Override
  	public String getBibtexCitation() {
  		return "@article{Boegel13-RMC,\n"+
  		"  number={1},\n"+
  		"  author={Marco B{\"o}gel and Hannes Hofmann and Joachim Hornegger and Rebecca Fahrig and Stefan Britzen and Andreas Maier},\n"+
  		"  keywords={cardiac reconstruction; c-arm ct; motion compensation; diaphragm tracking},\n"+
  		"  url={http://www5.informatik.uni-erlangen.de/Forschung/Publikationen/2013/Boegel13-RMC.pdf},\n"+
  		"  doi={10.1155/2013/520540},\n"+
  		"  journal={International Journal of Biomedical Imaging},\n"+
  		"  volume={2013},\n"+
  		"  title={{Respiratory Motion Compensation Using Diaphragm Tracking for Cone-Beam C-Arm CT: A Simulation and a Phantom Study}},\n"+
  		"  year={2013},\n"+
  		"  pages={1--10}\n" +
  		"}";
  	}
  
  
  
  
  	@Override
  	public String getMedlineCitation() {
  		return "Bögel M, Hofmann H, Hornegger J, Fahrig R, Britzen S, Maier A. Respiratory Motion Compensation Using Diaphragm Tracking for Cone-Beam C-Arm CT: A Simulation and a Phantom Study. International Journal of Biomedical Imaging, vol. 2013, no. 1, pp. 1-10, 2013 ";
  	}
  
  	public void waitForResult() {
  		cudaRun();
  	}
  
  	@Override
  	public void backproject(Grid2D projection, int projectionNumber)
  	{
  		appendProjection(projection, projectionNumber);
  	}
  
  	private void appendProjection(Grid2D projection, int projectionNumber){	
  		projections.add(projection, projectionNumber);
  		projectionsAvailable.add(new Integer(projectionNumber));
  	}
  
  	private synchronized void projectSingleProjection(int projectionNumber, int dimz){
  		// load projection matrix
  		initProjectionMatrix(projectionNumber);
  		// load projection
  		Grid2D projection = projections.get(projectionNumber); 
  		initProjectionData(projection);
  		if (!largeVolumeMode) {
  			projections.remove(projectionNumber);
  		}
  		// backproject for each slice
  		// CUDA Grids are only two dimensional!
  		int [] zed = new int[1];
  		int reconDimensionZ = dimz;
  		double voxelSpacingX = getGeometry().getVoxelSpacingX();
  		double voxelSpacingY = getGeometry().getVoxelSpacingY();
  		double voxelSpacingZ = getGeometry().getVoxelSpacingZ();
  
  		zed[0] = reconDimensionZ;
  		double motionfieldShifted[] = new double[200];
  
  		double val = Configuration.getGlobalConfiguration().getRespiratoryMotionFieldEntry(0);
  		for(int i = 0; i < 200; i++) {
  			if(i < 20){
  				motionfieldShifted[i] = 1.3*0.05*i*val;}
  			else
  				motionfieldShifted[i] =1.3*Configuration.getGlobalConfiguration().getRespiratoryMotionFieldEntry(i-20);
  		}
  
  		mot[projectionNumber] = -1.0* Configuration.getGlobalConfiguration().getRespiratoryMotionFieldEntry(projectionNumber);
  		//		Pointer diaMotion = Pointer.to(new int[]{(int) Math.round(-1.0*Configuration.getGlobalConfiguration().getRespiratoryMotionFieldEntry(projectionNumber)) });
  		//		Pointer diaMotion = Pointer.to(new int[]{(int) Math.round(-1.0*motionfieldShifted[projectionNumber]) });
  		System.out.println(Configuration.getGlobalConfiguration().getRespiratoryMotionFieldEntry(projectionNumber));
  		double motionfield[] = new double[200];
  		motionfield[0]=0.0;
  		motionfield[1]=-0.18582188544104383;
  		motionfield[2]=-0.3716437708821445;
  		motionfield[3]=-0.5574656563231883;
  		motionfield[4]=-0.7432875417642322;
  		motionfield[5]=-0.9291094272053044;
  		motionfield[6]=-1.1149313126463767;
  		motionfield[7]=-1.300753198087449;
  		motionfield[8]=-1.4865750835284928;
  		motionfield[9]=-1.672396968969565;
  		motionfield[10]=-1.8582188544106089;
  		motionfield[11]=-2.044040739851681;
  		motionfield[12]=-2.2298626252927534;
  		motionfield[13]=-2.415684510733797;
  		motionfield[14]=-2.6015063961748695;
  		motionfield[15]=-2.7873282816159417;
  		motionfield[16]=-2.9731501670569855;
  		motionfield[17]=-3.158972052498058;
  		motionfield[18]=-3.34479393793913;
  		motionfield[19]=-3.530615823380174;
  		motionfield[20]=-3.7164377088212177;
  		motionfield[21]=-3.9022595942623184;
  		motionfield[22]=-4.088081479703362;
  		motionfield[23]=-4.459705151586803;
  		motionfield[24]=-4.88441504816808;
  		motionfield[25]=-5.309124944749357;
  		motionfield[26]=-5.7338348413306335;
  		motionfield[27]=-6.15854473791191;
  		motionfield[28]=-6.583254634493187;
  		motionfield[29]=-7.007964531074435;
  		motionfield[30]=-7.432674427655741;
  		motionfield[31]=-7.857384324237017;
  		motionfield[32]=-8.282094220818294;
  		motionfield[33]=-8.706804117399543;
  		motionfield[34]=-9.131514013980848;
  		motionfield[35]=-9.556223910562096;
  		motionfield[36]=-9.980933807143401;
  		motionfield[37]=-10.40564370372465;
  		motionfield[38]=-10.830353600305955;
  		motionfield[39]=-11.255063496887203;
  		motionfield[40]=-11.679773393468508;
  		motionfield[41]=-12.104483290049757;
  		motionfield[42]=-12.529193186631034;
  		motionfield[43]=-12.953903083212339;
  		motionfield[44]=-13.378612979793587;
  		motionfield[45]=-13.767469952177692;
  		motionfield[46]=-14.127644585203996;
  		motionfield[47]=-14.487819218230328;
  		motionfield[48]=-14.84799385125666;
  		motionfield[49]=-15.208168484282993;
  		motionfield[50]=-15.568343117309325;
  		motionfield[51]=-15.928517750335658;
  		motionfield[52]=-16.288692383361962;
  		motionfield[53]=-16.648867016388294;
  		motionfield[54]=-17.009041649414655;
  		motionfield[55]=-17.36921628244096;
  		motionfield[56]=-17.72939091546729;
  		motionfield[57]=-18.089565548493624;
  		motionfield[58]=-18.449740181519957;
  		motionfield[59]=-18.80991481454629;
  		motionfield[60]=-19.170089447572593;
  		motionfield[61]=-19.530264080598926;
  		motionfield[62]=-19.890438713625258;
  		motionfield[63]=-20.25061334665159;
  		motionfield[64]=-20.610787979677923;
  		motionfield[65]=-20.970962612704227;
  		motionfield[66]=-21.33113724573056;
  		motionfield[67]=-21.58981084329045;
  		motionfield[68]=-21.645482369917403;
  		motionfield[69]=-21.701153896544355;
  		motionfield[70]=-21.756825423171307;
  		motionfield[71]=-21.81249694979826;
  		motionfield[72]=-21.86816847642521;
  		motionfield[73]=-21.92384000305219;
  		motionfield[74]=-21.979511529679144;
  		motionfield[75]=-22.035183056306096;
  		motionfield[76]=-22.090854582933048;
  		motionfield[77]=-22.14652610956003;
  		motionfield[78]=-22.20219763618698;
  		motionfield[79]=-22.25786916281396;
  		motionfield[80]=-22.313540689440913;
  		motionfield[81]=-22.369212216067865;
  		motionfield[82]=-22.424883742694817;
  		motionfield[83]=-22.48055526932177;
  		motionfield[84]=-22.53622679594872;
  		motionfield[85]=-22.591898322575673;
  		motionfield[86]=-22.647569849202654;
  		motionfield[87]=-22.703241375829606;
  		motionfield[88]=-22.758912902456558;
  		motionfield[89]=-22.791821477052054;
  		motionfield[90]=-22.642626435395727;
  		motionfield[91]=-22.493431393739428;
  		motionfield[92]=-22.34423635208313;
  		motionfield[93]=-22.19504131042683;
  		motionfield[94]=-22.045846268770504;
  		motionfield[95]=-21.896651227114205;
  		motionfield[96]=-21.747456185457906;
  		motionfield[97]=-21.598261143801608;
  		motionfield[98]=-21.449066102145252;
  		motionfield[99]=-21.299871060488954;
  		motionfield[100]=-21.150676018832655;
  		motionfield[101]=-21.001480977176357;
  		motionfield[102]=-20.852285935520058;
  		motionfield[103]=-20.70309089386376;
  		motionfield[104]=-20.553895852207432;
  		motionfield[105]=-20.404700810551105;
  		motionfield[106]=-20.255505768894807;
  		motionfield[107]=-20.106310727238508;
  		motionfield[108]=-19.95711568558218;
  		motionfield[109]=-19.807920643925883;
  		motionfield[110]=-19.658725602269584;
  		motionfield[111]=-19.509530560613257;
  		motionfield[112]=-19.255930117765246;
  		motionfield[113]=-18.98927899976826;
  		motionfield[114]=-18.722627881771274;
  		motionfield[115]=-18.455976763774288;
  		motionfield[116]=-18.18932564577733;
  		motionfield[117]=-17.922674527780345;
  		motionfield[118]=-17.65602340978336;
  		motionfield[119]=-17.3893722917864;
  		motionfield[120]=-17.122721173789444;
  		motionfield[121]=-16.85607005579243;
  		motionfield[122]=-16.58941893779547;
  		motionfield[123]=-16.322767819798514;
  		motionfield[124]=-16.0561167018015;
  		motionfield[125]=-15.789465583804542;
  		motionfield[126]=-15.522814465807528;
  		motionfield[127]=-15.25616334781057;
  		motionfield[128]=-14.989512229813613;
  		motionfield[129]=-14.722861111816627;
  		motionfield[130]=-14.456209993819641;
  		motionfield[131]=-14.189558875822684;
  		motionfield[132]=-13.922907757825698;
  		motionfield[133]=-13.656256639828712;
  		motionfield[134]=-13.369791282952463;
  		motionfield[135]=-13.073418806636539;
  		motionfield[136]=-12.777046330320616;
  		motionfield[137]=-12.480673854004749;
  		motionfield[138]=-12.184301377688797;
  		motionfield[139]=-11.887928901372902;
  		motionfield[140]=-11.591556425056979;
  		motionfield[141]=-11.295183948741084;
  		motionfield[142]=-10.998811472425189;
  		motionfield[143]=-10.702438996109265;
  		motionfield[144]=-10.406066519793399;
  		motionfield[145]=-10.109694043477475;
  		motionfield[146]=-9.813321567161552;
  		motionfield[147]=-9.516949090845657;
  		motionfield[148]=-9.220576614529762;
  		motionfield[149]=-8.924204138213838;
  		motionfield[150]=-8.627831661897943;
  		motionfield[151]=-8.33145918558202;
  		motionfield[152]=-8.035086709266125;
  		motionfield[153]=-7.73871423295023;
  		motionfield[154]=-7.442341756634335;
  		motionfield[155]=-7.145969280318411;
  		motionfield[156]=-6.879726764496922;
  		motionfield[157]=-6.651146699293491;
  		motionfield[158]=-6.42256663409006;
  		motionfield[159]=-6.193986568886601;
  		motionfield[160]=-5.965406503683141;
  		motionfield[161]=-5.7368264384797385;
  		motionfield[162]=-5.508246373276279;
  		motionfield[163]=-5.279666308072848;
  		motionfield[164]=-5.051086242869417;
  		motionfield[165]=-4.822506177665957;
  		motionfield[166]=-4.593926112462498;
  		motionfield[167]=-4.365346047259095;
  		motionfield[168]=-4.1367659820556355;
  		motionfield[169]=-3.908185916852233;
  		motionfield[170]=-3.6796058516487733;
  		motionfield[171]=-3.4510257864453138;
  		motionfield[172]=-3.2224457212418542;
  		motionfield[173]=-2.9938656560384516;
  		motionfield[174]=-2.765285590834992;
  		motionfield[175]=-2.536705525631561;
  		motionfield[176]=-2.30812546042813;
  		motionfield[177]=-2.0795453952246987;
  		motionfield[178]=-1.8827432910657933;
  		motionfield[179]=-1.7971640505627988;
  		motionfield[180]=-1.7115848100598043;
  		motionfield[181]=-1.6260055695568099;
  		motionfield[182]=-1.5404263290538438;
  		motionfield[183]=-1.454847088550821;
  		motionfield[184]=-1.3692678480478548;
  		motionfield[185]=-1.283688607544832;
  		motionfield[186]=-1.1981093670418659;
  		motionfield[187]=-1.1125301265388714;
  		motionfield[188]=-1.026950886035877;
  		motionfield[189]=-0.9413716455328824;
  		motionfield[190]=-0.855792405029888;
  		motionfield[191]=-0.7702131645268935;
  		motionfield[192]=-0.6846339240239274;
  		motionfield[193]=-0.5990546835209329;
  		motionfield[194]=-0.5134754430179385;
  		motionfield[195]=-0.427896202514944;
  		motionfield[196]=-0.3423169620119779;
  		motionfield[197]=-0.256737721508955;
  		motionfield[198]=-0.17115848100598896;
  		motionfield[199]=-0.08557924050296606;
  
  		Pointer diaMotion = Pointer.to(new int[]{(int) (Math.round(motionfield[projectionNumber])) });
  
  		if(projectionNumber == 199){
  			VisualizationUtil.createPlot(motionfield, "motionfield real", "frame", "z").show();		
  			VisualizationUtil.createPlot(mot, "motionfield scaled", "frame", "z").show();
  			VisualizationUtil.createPlot(motionfieldShifted, "shifted", "frame", "z").show();
  		}
  		//Heart respiration compensation
  
  
  		/*
  		 * Lung Compensation		
  		 */   /*  MotionFieldReader motion = null;
  		try {
  			motion = new MotionFieldReader();
  		} catch (Exception e) {
  			// TODO Auto-generated catch block
  			e.printStackTrace();
  		}*/
  		//double diaMot = Configuration.getGlobalConfiguration().getRespiratoryMotionFieldEntry(projectionNumber);
  
  		//double diaPos = Configuration.getGlobalConfiguration().getDiaphragmPositionFieldEntry(projectionNumber);
  		//float fMot = (float) (diaMot/Configuration.getGlobalConfiguration().getMaxMotion());
  		//float amplitude = 0.5f;
  		//Pointer m = Pointer.to(new float[]{(float) motion.getGlobalCompensationLinearScaling()});
  		//Pointer m = Pointer.to(new float[]{0.f});
  		//Pointer m = Pointer.to(new float[]{(float) motion.getGlobalCompensationLinearMinMax(fMot, diaPos, reconDimensionZ)});
  		//System.out.println(motion.getGlobalCompensationLinearMinMax(fMot, diaPos, reconDimensionZ*voxelSpacingZ-offsetZ));
  		//Pointer m = Pointer.to(new float[]{(float) motion.getInterpolatedGlobalCompensationLinearScaling(fMot)});
  		//System.out.println(motion.getInterpolatedGlobalCompensationLinearScaling(fMot));
  		//Pointer diaPosition = Pointer.to(new float[]{(float) diaPos});
  		//Pointer diaMotion = Pointer.to(new float[]{(float)(-amplitude*fMot) });	
  		//for old constant version set m to zero		
  
  		//	Pointer interX = Pointer.to(motion.getBilinearInterpolationDataPosition(diaMot));
  		//	Pointer interY = Pointer.to(motion.getBilinearInterpolationDataMotion(diaMot));
  
  		//		System.out.println((int) Math.round(Configuration.getGlobalConfiguration().getRespiratoryMotionFieldEntry(projectionNumber)));
  		Pointer dOut = Pointer.to(volumePointer);
  		Pointer pWidth = Pointer.to(new int[]{(int) lineOffset});
  		Pointer pZOffset = Pointer.to(zed);
  		float [] vsx = new float[]{(float) voxelSpacingX};
  		Pointer pvsx = Pointer.to(vsx);
  		Pointer pvsy = Pointer.to(new float[]{(float) voxelSpacingY});
  		Pointer pvsz = Pointer.to(new float[]{(float) voxelSpacingZ});
  		Pointer pox = Pointer.to(new float[]{(float) offsetX});
  		Pointer poy = Pointer.to(new float[]{(float) offsetY});
  		Pointer poz = Pointer.to(new float[]{(float) offsetZ});
  
  		int offset = 0;
  		//System.out.println(dimz + " " + zed[0] + " " + offsetZ + " " + voxelSpacingZ);
  
  		offset = CUDAUtil.align(offset, Sizeof.POINTER);
  		JCudaDriver.cuParamSetv(function, offset, dOut, Sizeof.POINTER);
  		offset += Sizeof.POINTER;
  		/*// new for lung motion
  		offset = CUDAUtil.align(offset,Sizeof.FLOAT);
  		JCudaDriver.cuParamSetv(function, offset, m, Sizeof.FLOAT);
  		offset += Sizeof.FLOAT;
  		//new for lung motion
  		offset = CUDAUtil.align(offset,Sizeof.FLOAT);
  		JCudaDriver.cuParamSetv(function, offset, diaPosition, Sizeof.FLOAT);
  		offset += Sizeof.FLOAT;
  		 */
  		offset = CUDAUtil.align(offset,Sizeof.INT);
  		JCudaDriver.cuParamSetv(function, offset, diaMotion, Sizeof.INT);
  		offset += Sizeof.INT;
  
  		offset = CUDAUtil.align(offset, Sizeof.INT);
  		JCudaDriver.cuParamSetv(function, offset, pWidth, Sizeof.INT);
  		offset += Sizeof.INT;
  
  		offset = CUDAUtil.align(offset, Sizeof.INT);
  		JCudaDriver.cuParamSetv(function, offset, pZOffset, Sizeof.INT);
  		offset += Sizeof.INT;
  
  		offset = CUDAUtil.align(offset, Sizeof.FLOAT);
  		JCudaDriver.cuParamSetv(function, offset, pvsx, Sizeof.FLOAT);
  		offset += Sizeof.FLOAT;
  
  		offset = CUDAUtil.align(offset, Sizeof.FLOAT);
  		JCudaDriver.cuParamSetv(function, offset, pvsy, Sizeof.FLOAT);
  		offset += Sizeof.FLOAT;
  
  		offset = CUDAUtil.align(offset, Sizeof.FLOAT);
  		JCudaDriver.cuParamSetv(function, offset, pvsz, Sizeof.FLOAT);
  		offset += Sizeof.FLOAT;
  
  
  		offset = CUDAUtil.align(offset, Sizeof.FLOAT);
  		JCudaDriver.cuParamSetv(function, offset, pox, Sizeof.FLOAT);
  		offset += Sizeof.FLOAT;
  
  		offset = CUDAUtil.align(offset, Sizeof.FLOAT);
  		JCudaDriver.cuParamSetv(function, offset, poy, Sizeof.FLOAT);
  		offset += Sizeof.FLOAT;
  
  		offset = CUDAUtil.align(offset, Sizeof.FLOAT);
  		JCudaDriver.cuParamSetv(function, offset, poz, Sizeof.FLOAT);
  		offset += Sizeof.FLOAT;
  
  		JCudaDriver.cuParamSetSize(function, offset);
  
  		// Call the CUDA kernel, writing the results into the volume which is pointed at
  		JCudaDriver.cuFuncSetBlockShape(function, bpBlockSize[0], bpBlockSize[1], 1);
  		JCudaDriver.cuLaunchGrid(function, gridSize.x, gridSize.y);
  		JCudaDriver.cuCtxSynchronize();
  
  	}
  	public void cudaRun() {
  		try {
  			while (projectionsAvailable.size() > 0) {
  				Thread.sleep(CONRAD.INVERSE_SPEEDUP);
  				if (showStatus) {
  					float status = (float) (1.0 / projections.size());
  					if (largeVolumeMode) {
  						IJ.showStatus("Streaming Projections to CUDA Buffer");
  					} else {
  						IJ.showStatus("Backprojecting with CUDA");
  					}
  					IJ.showProgress(status);
  				}
  				if (!largeVolumeMode) {			
  					workOnProjectionData();
  				} else {
  					checkProjectionData();
  				}
  			}
  			System.out.println("large Volume " + largeVolumeMode);
  			if (largeVolumeMode){
  				// we have collected all projections.
  				// now we can reconstruct subvolumes and stich them together.
  				int reconDimensionX = getGeometry().getReconDimensionX();
  				int reconDimensionY = getGeometry().getReconDimensionY();
  				int reconDimensionZ = getGeometry().getReconDimensionZ();
  				double voxelSpacingX = getGeometry().getVoxelSpacingX();
  				double voxelSpacingY = getGeometry().getVoxelSpacingY();
  				double voxelSpacingZ = getGeometry().getVoxelSpacingZ();
  				useVOImap = false;
  				initialize(projections.get(0));
  				double originalOffsetZ = offsetZ;
  				double originalReconDimZ = reconDimensionZ;
  				reconDimensionZ = subVolumeZ;
  				int memorysize = reconDimensionX * reconDimensionY * subVolumeZ * Sizeof.FLOAT;
  				int maxProjectionNumber = projections.size();
  				float all = nSteps * maxProjectionNumber*2;
  				for (int n =0; n < nSteps; n++){ // For each subvolume
  					// set all to 0;
  					Arrays.fill(h_volume, 0);
  					JCudaDriver.cuMemcpyHtoD(volumePointer, Pointer.to(h_volume), memorysize);
  					offsetZ = originalOffsetZ - (reconDimensionZ*voxelSpacingZ*n);
  					for (int p = 0; p < maxProjectionNumber; p ++){ // For all projections
  						float currentStep = (n*maxProjectionNumber*2) + p;
  						if (showStatus) {
  							IJ.showStatus("Backprojecting with CUDA");
  							IJ.showProgress(currentStep/all);
  						}
  						//System.out.println("Current: " + p);
  						try {
  							projectSingleProjection(p, reconDimensionZ);
  						} catch (Exception e){
  							System.out.println("Backprojection of projection " + p + " was not successful.");
  							e.printStackTrace();
  						}
  					}
  					// Gather volume
  					JCudaDriver.cuMemcpyDtoH(Pointer.to(h_volume), volumePointer, memorysize);
  
  					// move data to ImagePlus;
  					if (projectionVolume != null) {
  						for (int k = 0; k < reconDimensionZ; k++){
  							int index = (n*subVolumeZ) + k;
  							if (showStatus) {
  								float currentStep = (n*maxProjectionNumber*2) + maxProjectionNumber + k;
  								IJ.showStatus("Fetching Volume from CUDA");
  								IJ.showProgress(currentStep/all);
  							}
  							if (index < originalReconDimZ) {
  								for (int j = 0; j < projectionVolume.getSize()[1]; j++){
  									for (int i = 0; i < projectionVolume.getSize()[0]; i++){
  										float value = h_volume[(((projectionVolume.getSize()[1] * k) + j) * projectionVolume.getSize()[0]) + i];
  										double[][] voxel = new double [4][1];
  										voxel[0][0] = (voxelSpacingX * i) - offsetX;
  										voxel[1][0] = (voxelSpacingY * j) - offsetY;
  										voxel[2][0] = (voxelSpacingZ * index) - originalOffsetZ;
  										if (interestedInVolume.contains(voxel[0][0], voxel[1][0], voxel[2][0])) {
  											projectionVolume.setAtIndex(i, j, index, value);
  										} else {
  											projectionVolume.setAtIndex(i, j, index, 0);
  										}
  									}
  								}
  							}
  						}
  					}
  				}
  			}
  
  
  		} catch (InterruptedException e) {
  
  			e.printStackTrace();
  		}
  		if (showStatus) IJ.showProgress(1.0);
  		unload();
  		if (debug) System.out.println("Unloaded");
  	}
  
  	private synchronized void workOnProjectionData(){
  		if (projectionsAvailable.size() > 0){
  			Integer current = projectionsAvailable.get(0);
  			projectionsAvailable.remove(0);
  			projectSingleProjection(current.intValue(),  
  					getGeometry().getReconDimensionZ());
  			projectionsDone.add(current);
  		}	
  	}
  
  	private synchronized void checkProjectionData(){
  		if (projectionsAvailable.size() > 0){
  			Integer current = projectionsAvailable.get(0);
  			projectionsAvailable.remove(current);
  			projectionsDone.add(current);
  		}	
  	}
  
  	public void reconstructOffline(ImagePlus imagePlus) throws Exception {
  		ImagePlusDataSink sink = new ImagePlusDataSink();
  		configure();
  		init();
  		for (int i = 0; i < imagePlus.getStackSize(); i++){
  			backproject(ImageUtil.wrapImageProcessor(imagePlus.getStack().getProcessor(i+1)), i);
  		}
  		waitForResult();
  		if (Configuration.getGlobalConfiguration().getUseHounsfieldScaling()) applyHounsfieldScaling();
  		int [] size = projectionVolume.getSize();
  		System.out.println(size [0] + " " + size [1] + " " + size[2]);
  		for (int k = 0; k < projectionVolume.getSize()[2]; k++){
  			sink.process(projectionVolume.getSubGrid(k), k);
  		}
  		sink.close();
  		ImagePlus revan = ImageUtil.wrapGrid3D(sink.getResult(),"Compensated Reconstruction of" + imagePlus.getTitle());
  		revan.setTitle(imagePlus.getTitle() + " reconstructed");
  		revan.show();
  		reset();
  	}
  
  	@Override
  	protected void reconstruct() throws Exception {
  		init();
  		for (int i = 0; i < nImages; i++){
  			backproject(inputQueue.get(i), i);
  		}
  		waitForResult();
  		if (Configuration.getGlobalConfiguration().getUseHounsfieldScaling()) applyHounsfieldScaling();
  		int [] size = projectionVolume.getSize();
  
  		for (int k = 0; k < size[2]; k++){
  			sink.process(projectionVolume.getSubGrid(k), k);
  		}
  		sink.close();
  	}
  
  	@Override
  	public String getToolName(){
  		return "CUDA Compensated Backprojector";
  	}
  
  }
  
  /*
   * Copyright (C) 2010-2014 - Marco Bgel 
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
  */
  