content: package edu.stanford.rsl.tutorial.cone;
  
  
  import java.io.IOException;
  import java.nio.FloatBuffer;
  
  import com.jogamp.opencl.CLBuffer;
  import com.jogamp.opencl.CLCommandQueue;
  import com.jogamp.opencl.CLContext;
  import com.jogamp.opencl.CLDevice;
  import com.jogamp.opencl.CLImage3d;
  import com.jogamp.opencl.CLImageFormat;
  import com.jogamp.opencl.CLKernel;
  import com.jogamp.opencl.CLProgram;
  import com.jogamp.opencl.CLImageFormat.ChannelOrder;
  import com.jogamp.opencl.CLImageFormat.ChannelType;
  import com.jogamp.opencl.CLMemory.Mem;
  
  import edu.stanford.rsl.conrad.data.numeric.Grid2D;
  import edu.stanford.rsl.conrad.data.numeric.Grid3D;
  import edu.stanford.rsl.conrad.data.numeric.InterpolationOperators;
  import edu.stanford.rsl.conrad.data.numeric.opencl.OpenCLGrid2D;
  import edu.stanford.rsl.conrad.data.numeric.opencl.OpenCLGrid3D;
  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;
  import edu.stanford.rsl.conrad.opencl.OpenCLUtil;
  import edu.stanford.rsl.conrad.opencl.TestOpenCL;
  import edu.stanford.rsl.conrad.utils.CONRAD;
  import edu.stanford.rsl.conrad.utils.Configuration;
  
  public class ConeBeamProjector {
  
  	final boolean debug = false;
  	final boolean verbose = false;
  
  	static int bpBlockSize[] = {16, 16};
  
  	//opencl variables
  	protected CLContext context;
  	protected CLDevice device;
  	private CLProgram program;
  	private CLImageFormat format;
  	private CLBuffer<FloatBuffer> gVolumeEdgeMaxPoint = null;
  	private CLBuffer<FloatBuffer> gVolumeEdgeMinPoint = null;
  	private CLBuffer<FloatBuffer> gVoxelElementSize = null;
  	protected CLBuffer<FloatBuffer> gInvARmatrix = null;
  	protected CLBuffer<FloatBuffer> gSrcPoint = null;
  	private CLBuffer<FloatBuffer> sinogram = null;
  	private CLImage3d<FloatBuffer> imageGrid = null;
  	protected CLCommandQueue queue = null;
  	private CLKernel kernelFunction;
  	// Length of arrays to process
  	int localWorkSize;
  	int globalWorkSizeU;
  	int globalWorkSizeV; 
  	
  	//imaging variables
  	private int width;
  	private int height;
  	protected Trajectory geometry;
  	private int currentStep = 0;
  	private float [] voxelSize = null;
  	private float [] volumeSize = null;
  	private float [] volumeEdgeMinPoint = null;
  	private float [] volumeEdgeMaxPoint = null;
  	// buffer for 3D volume:
  	private int subVolumeZ;
  
  	public ConeBeamProjector() {
  		configure();
  		initCL();
  	}
  	
  	private void initCL(){
  		context = OpenCLUtil.getStaticContext();
  		device = context.getMaxFlopsDevice();
  		program = null;
  		// initialize the program
  		if (program==null || !program.getContext().equals(context)){
  			try {
  				program = context.createProgram(TestOpenCL.class.getResourceAsStream("projectCL.cl")).build();
  			} catch (IOException e) {
  				// TODO Auto-generated catch block
  				e.printStackTrace();
  			}
  		}
  
  		// create image from input grid
  		format = new CLImageFormat(ChannelOrder.INTENSITY, ChannelType.FLOAT);
  		
  		// Length of arrays to process
  		localWorkSize = Math.min(device.getMaxWorkGroupSize(), 16); // Local work size dimensions
  		globalWorkSizeU = OpenCLUtil.roundUp(localWorkSize, width); // rounded up to the nearest multiple of localWorkSize
  		globalWorkSizeV = OpenCLUtil.roundUp(localWorkSize, height); // rounded up to the nearest multiple of localWorkSize
  		
  		queue = device.createCommandQueue();
  		kernelFunction = program.createCLKernel("projectKernel");
  		setEdgeMaxima();
  		prepareAllProjections();
  
  		gVoxelElementSize = context.createFloatBuffer(voxelSize.length, Mem.READ_ONLY);
  		gVoxelElementSize.getBuffer().put(voxelSize);
  		gVoxelElementSize.getBuffer().rewind();
  		
  		queue
  		.putWriteBuffer(gVoxelElementSize,true)
  		.putWriteBuffer(gVolumeEdgeMinPoint, true)
  		.putWriteBuffer(gVolumeEdgeMaxPoint, true)
  		.finish();
  		
  	}
  
  	public Grid2D projectPixelDriven(Grid3D grid, int projIdx) {
  		geometry = Configuration.getGlobalConfiguration().getGeometry();
  		int maxV = geometry.getDetectorHeight();
  		int maxU = geometry.getDetectorWidth();
  		int imgSizeX = geometry.getReconDimensionX();
  		int imgSizeY = geometry.getReconDimensionY();
  		int imgSizeZ = geometry.getReconDimensionZ();
  		Projection[] projMats = geometry.getProjectionMatrices();
  		int maxProjs = geometry.getProjectionStackSize();
  		if(projIdx+1 > maxProjs || 0 > projIdx){
  			System.err.println("ConeBeamProjector: Invalid projection index");
  			return null;
  		}
  		Grid2D sino = new Grid2D(maxU,maxV); //
  		double spacingX = geometry.getVoxelSpacingX();
  		double spacingY = geometry.getVoxelSpacingY();
  		double spacingZ = geometry.getVoxelSpacingZ();
  		double originX = -geometry.getOriginX();
  		double originY = -geometry.getOriginY();
  		double originZ = -geometry.getOriginZ();
  		for (int x = 0; x < imgSizeX - 1; x++) {
  			double xTrans = x * spacingX - originX;
  			for (int y = 0; y < imgSizeY - 1; y++) {
  				double yTrans = y * spacingY - originY;
  				for (int z = 0; z < imgSizeZ - 1; z++) {
  					SimpleVector point3d = new SimpleVector(xTrans, yTrans, z * spacingZ - originZ, 1);
  					// for(int p = 0; p < maxProjs; p++) {
  					int p = projIdx; //
  					SimpleVector point2d = SimpleOperators.multiply(projMats[p].computeP(), point3d);
  					double coordU = point2d.getElement(0) / point2d.getElement(2);
  					double coordV = point2d.getElement(1) / point2d.getElement(2);
  					if (coordU >= maxU - 1 || coordV >= maxV - 1 || coordU <= 0 || coordV <= 0)
  						continue;
  					float val = grid.getAtIndex(x, y, z);
  					InterpolationOperators.addInterpolateLinear(sino, coordU, coordV, val); //
  					// }
  				}
  			}
  		}
  
  		return sino;
  	}
  
  
  	public Grid3D projectPixelDriven(Grid3D grid) {
  		geometry = Configuration.getGlobalConfiguration().getGeometry();
  		int maxV = geometry.getDetectorHeight();
  		int maxU = geometry.getDetectorWidth();
  		int imgSizeX = geometry.getReconDimensionX();
  		int imgSizeY = geometry.getReconDimensionY();
  		int imgSizeZ = geometry.getReconDimensionZ();
  		Projection[] projMats = geometry.getProjectionMatrices();
  		int maxProjs = geometry.getProjectionStackSize();
  		Grid3D sino = new Grid3D(maxU,maxV,maxProjs);
  		double spacingX = geometry.getVoxelSpacingX();
  		double spacingY = geometry.getVoxelSpacingY();
  		double spacingZ = geometry.getVoxelSpacingZ();
  		double originX = -geometry.getOriginX();
  		double originY = -geometry.getOriginY();
  		double originZ = -geometry.getOriginZ();
  		for(int x = 0; x < imgSizeX-1; x++) {
  			double xTrans = x*spacingX-originX;
  			for(int y = 0; y < imgSizeY-1 ; y++) {
  				double yTrans = y*spacingY-originY;
  				for(int z = 0; z < imgSizeZ-1; z++) {
  					SimpleVector point3d = new SimpleVector(xTrans, yTrans, z*spacingZ-originZ, 1);
  					for(int p = 0; p < maxProjs; p++) {
  						SimpleVector point2d = SimpleOperators.multiply(projMats[p].computeP(), point3d);
  						double coordU = point2d.getElement(0) / point2d.getElement(2);
  						double coordV = point2d.getElement(1) / point2d.getElement(2);
  						if (coordU >= maxU-1 || coordV>= maxV -1 || coordU <=0 || coordV <= 0)
  							continue;
  						float val = grid.getAtIndex(x, y, z);
  						InterpolationOperators.addInterpolateLinear(sino, coordU, coordV, p, val);
  					}
  				}
  			}
  		}
  
  		return sino;
  	}
  
  
  	private void setEdgeMaxima(){
  
  		int test = currentStep;
  		volumeEdgeMaxPoint[2] = (float) ((test * subVolumeZ) + subVolumeZ -0.5 - CONRAD.SMALL_VALUE);
  		volumeEdgeMinPoint[2] = (float) ((test * subVolumeZ) -0.5 - CONRAD.SMALL_VALUE);
  		
  		gVolumeEdgeMaxPoint = context.createFloatBuffer(volumeEdgeMaxPoint.length, Mem.READ_ONLY);
  		gVolumeEdgeMinPoint = context.createFloatBuffer(volumeEdgeMinPoint.length, Mem.READ_ONLY);
  
  		gVolumeEdgeMaxPoint.getBuffer().put(volumeEdgeMaxPoint);
  		gVolumeEdgeMinPoint.getBuffer().put(volumeEdgeMinPoint);
  
  		gVolumeEdgeMaxPoint.getBuffer().rewind();
  		gVolumeEdgeMinPoint.getBuffer().rewind();
  	}
  
  	protected void prepareAllProjections(){
  
  		float [] cann = new float[3*4];
  		float [] invAR = new float[3*3];
  		float [] srcP = new float[3];
  
  		//if (gInvARmatrix == null)
  		gInvARmatrix = context.createFloatBuffer(invAR.length*geometry.getNumProjectionMatrices(), Mem.READ_ONLY);
  		//if (gSrcPoint == null)
  		gSrcPoint = context.createFloatBuffer(srcP.length*geometry.getNumProjectionMatrices(), Mem.READ_ONLY);
  
  		for (int i=0; i < geometry.getNumProjectionMatrices(); ++i){
  			SimpleMatrix projMat = geometry.getProjectionMatrix(i).computeP();
  			double [][] mat = new double [3][4];
  			projMat.copyTo(mat);
  			computeCanonicalProjectionMatrix(cann, invAR, srcP, new Jama.Matrix(mat));
  
  			gInvARmatrix.getBuffer().put(invAR);
  			gSrcPoint.getBuffer().put(srcP);
  		}
  
  		gInvARmatrix.getBuffer().rewind();
  		gSrcPoint.getBuffer().rewind();
  
  		queue
  		.putWriteBuffer(gSrcPoint, true)
  		.putWriteBuffer(gInvARmatrix, true)
  		.finish();
  	}
  
  	public void computeCanonicalProjectionMatrix(float [] canonicalProjMatrix, float [] invARmatrix, float [] srcPoint, Jama.Matrix projectionMatrix){
  
  		double [] du = {geometry.getPixelDimensionX(), 0};
  		double [] dv = {0, geometry.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)); 
  
  	}
  
  	public static Jama.Matrix computeSrcPt(Jama.Matrix projectionMatrix, Jama.Matrix invertedProjMatrix) {
  		Jama.Matrix at = projectionMatrix.getMatrix(0, 2, 3, 3);
  		//at = at.times(-1.0);
  		return invertedProjMatrix.times(at);
  	}
  
  	public void configure(){
  		Configuration.loadConfiguration();
  		geometry = Configuration.getGlobalConfiguration().getGeometry();
  		voxelSize = new float [3];
  		volumeSize = new float [3];
  
  		voxelSize[0] = (float) geometry.getVoxelSpacingX();
  		voxelSize[1] = (float) geometry.getVoxelSpacingY();
  		voxelSize[2] = (float) geometry.getVoxelSpacingZ();
  
  		volumeSize[0] = geometry.getReconDimensionX();
  		volumeSize[1] = geometry.getReconDimensionY();
  		volumeSize[2] = geometry.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 = geometry.getDetectorWidth();
  		height = geometry.getDetectorHeight();
  
  		subVolumeZ = (int) volumeSize[2];
  		
  		if (debug) System.out.println("Projection Matrices: " + geometry.getNumProjectionMatrices());
  	}
  
  	public void unload(){
  
  		if(sinogram != null && !sinogram.isReleased())
  			sinogram.release();
  		if(imageGrid != null && !imageGrid.isReleased())
  			imageGrid.release();
  
  	}
  
  	public int getMaxProjections() {
  		return geometry.getProjectionStackSize();
  	}
  
  	public void projectRayDrivenCL(OpenCLGrid2D[] sinoCL, OpenCLGrid3D gridCL){
  
  		imageGrid = context.createImage3d(gridCL.getDelegate().getCLBuffer().getBuffer(), (int)volumeSize[0], (int)volumeSize[1], (int)volumeSize[2],format, Mem.READ_ONLY);
  
  		queue
  		.putCopyBufferToImage(gridCL.getDelegate().getCLBuffer(), imageGrid)
  		.finish();
  
  		for(int p = 0; p < geometry.getProjectionStackSize(); p++) {
  	
  			kernelFunction
  			.putArg(sinoCL[p].getDelegate().getCLBuffer())
  			.putArg(width)
  			.putArg(height)
  			.putArg(1.f)
  			.putArg(imageGrid)
  			.putArg(gVoxelElementSize)
  			.putArg(gVolumeEdgeMinPoint)
  			.putArg(gVolumeEdgeMaxPoint)
  			.putArg(gSrcPoint)
  			.putArg(gInvARmatrix)
  			.putArg(p);
  
  			queue
  			.put2DRangeKernel(kernelFunction, 0, 0, globalWorkSizeU, globalWorkSizeV,localWorkSize, localWorkSize)
  			.finish();
  
  			kernelFunction.rewind();
  			sinoCL[p].getDelegate().notifyDeviceChange();
  			sinoCL[p].getDelegate().prepareForHostOperation();
  		}
  		
  		kernelFunction.release();
  		queue.release();
  	}
  	//2d method
  	/**
  	 * loads the actual OpenCL kernel and performs the projection
  	 * @param projectionNumber the projection number.
  	 * @return the image as image processor
  	 */
  
  	public void projectRayDrivenCL(OpenCLGrid2D sinoCL, OpenCLGrid3D gridCL, int projIdx){
  		
  		imageGrid = context.createImage3d(gridCL.getDelegate().getCLBuffer().getBuffer(), (int)gridCL.getSize()[0], (int)gridCL.getSize()[1], (int)gridCL.getSize()[2],format, Mem.READ_ONLY);
  
  		queue
  		.putCopyBufferToImage(gridCL.getDelegate().getCLBuffer(), imageGrid)
  		.finish();
  
  		kernelFunction
  		.putArg(sinoCL.getDelegate().getCLBuffer())
  		.putArg(width)
  		.putArg(height)
  		.putArg(1.f)
  		.putArg(imageGrid)
  		.putArg(gVoxelElementSize)
  		.putArg(gVolumeEdgeMinPoint)
  		.putArg(gVolumeEdgeMaxPoint)
  		.putArg(gSrcPoint)
  		.putArg(gInvARmatrix)
  		.putArg(projIdx);
  	
  		queue
  		.put2DRangeKernel(kernelFunction, 0, 0, globalWorkSizeU, globalWorkSizeV,localWorkSize, localWorkSize).finish();
  
  		kernelFunction.rewind();
  		sinoCL.getDelegate().notifyDeviceChange();
  	}
  
  	public Grid2D projectRayDrivenCL(Grid3D grid, int projIdx) throws Exception {
  		configure();
  
  		OpenCLGrid3D gridCL = new OpenCLGrid3D(grid);
  		OpenCLGrid2D sinoCL = new OpenCLGrid2D(new Grid2D(width,height));
  		sinoCL.getDelegate().prepareForDeviceOperation();
  
  		projectRayDrivenCL(sinoCL, gridCL, projIdx);
  
  		gridCL.release();
  
  		Grid2D sino = new Grid2D(sinoCL);
  		unload();
  		return sino;
  	}
  	
  	public Grid3D projectRayDrivenCL(Grid3D grid) {
  
  		configure();
  
  		OpenCLGrid3D gridCL = new OpenCLGrid3D(grid);
  		OpenCLGrid2D [] sinoCL = new OpenCLGrid2D[geometry.getProjectionStackSize()];
  		for (int i=0; i < geometry.getProjectionStackSize(); i++) {
  			sinoCL[i] = new OpenCLGrid2D(new Grid2D(width,height));
  			sinoCL[i].getDelegate().prepareForDeviceOperation();
  		}
  		
  		projectRayDrivenCL(sinoCL, gridCL);
  		gridCL.release();
  
  		Grid3D sino = new Grid3D(width,height,geometry.getProjectionStackSize());
  		for(int i = 0; i<sinoCL.length; i++){
  			sino.setSubGrid(i, new Grid2D(sinoCL[i]));
  			sinoCL[i].release();
  		}
  
  		unload();
  		return sino;
  
  	}
  
  	// calculates the sinogram out of a volume grid at the projectionindex projIdx
  	public void fastProjectRayDrivenCL(OpenCLGrid2D sinoCL, OpenCLGrid3D grid, int projIdx) {
  		sinoCL.getDelegate().prepareForDeviceOperation();
  		grid.getDelegate().prepareForDeviceOperation();
  		
  		projectRayDrivenCL(sinoCL, grid, projIdx);
  
  		unload();
  	}
  	
  	// calculates the sinogram out of a volume grid at the projectionindex projIdx
  	public void fastProjectRayDrivenCL(OpenCLGrid3D sinoCL, OpenCLGrid3D grid) throws Exception {
  
  		OpenCLGrid2D sinoCLBuffer = new OpenCLGrid2D(new Grid2D(width,height));
  		sinoCLBuffer.setOrigin(0,0);
  		sinoCLBuffer.setSpacing(geometry.getPixelDimensionX(),geometry.getPixelDimensionY());
  		sinoCL.getDelegate().prepareForDeviceOperation();
  
  		for(int pIdx = 0; pIdx < geometry.getProjectionStackSize(); pIdx++){
  			sinoCLBuffer.getDelegate().prepareForDeviceOperation();
  			fastProjectRayDrivenCL(sinoCLBuffer,grid,pIdx);
  			sinoCLBuffer.getDelegate().notifyDeviceChange();
  			queue.putCopyBuffer(sinoCLBuffer.getDelegate().getCLBuffer(), sinoCL.getDelegate().getCLBuffer(),0,pIdx*(int)sinoCLBuffer.getDelegate().getCLBuffer().getCLSize(),sinoCLBuffer.getDelegate().getCLBuffer().getCLSize(),null).finish();
  		}
  
  		sinoCLBuffer.release();
  		sinoCL.getDelegate().notifyDeviceChange();
  
  		unload();
  	}
  
  }
  /*
   * Copyright (C) 2010-2014 Andreas Maier
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
   */
  