package edu.stanford.rsl.conrad.cuda; //TODO: Use our own matrices instead of Jama.Matrix import java.util.ArrayList; 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.SimpleOperators; 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.FileUtil; import edu.stanford.rsl.conrad.utils.ImageUtil; import edu.stanford.rsl.conrad.utils.XmlUtils; /** * 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, Choi, Martin * */ public class CUDAForwardProjectorWithMotion 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; } /** * The XML filename where the rigid motion parameters are stored */ private String rotationTranslationFilename = null; 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]; // load the motion transforms // load data from XML file SimpleMatrix[] motion = readInMotionMatrices(); SimpleMatrix projMat = geometry.getProjectionMatrix(projectionNumber).computeP(); double [][] mat = new double [3][4]; SimpleOperators.multiplyMatrixProd(projMat, motion[projectionNumber]).copyTo(mat); // apply motion // projMat.copyTo(mat);// No motion applied 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); } public SimpleMatrix[] readInMotionMatrices() { SimpleMatrix[] motion = new SimpleMatrix[geometry.getNumProjectionMatrices()]; // load data from XML file ArrayList<double[][][]> RotTrans = null; try { if (rotationTranslationFilename == null) rotationTranslationFilename = FileUtil.myFileChoose(".xml", false); RotTrans = (ArrayList<double[][][]>)XmlUtils.importFromXML(rotationTranslationFilename); } catch (Exception e) { // TODO Auto-generated catch block e.printStackTrace(); } for (int i = 0; i < geometry.getNumProjectionMatrices(); i++) { motion[i] = new SimpleMatrix(4,4); SimpleMatrix rotation = new SimpleMatrix(RotTrans.get(0)[i]); SimpleMatrix translation = new SimpleMatrix(RotTrans.get(1)[i]); // motion[i].setSubMatrixValue(0, 0, rotation.transposed()); // motion[i].setSubMatrixValue(0, 3, translation.multipliedBy(-1)); motion[i].setSubMatrixValue(0, 0, rotation); //Correct motion[i].setSubMatrixValue(0, 3, translation); //Correct motion[i].setRowValue(3, new SimpleVector(0,0,0,1)); } return motion; } /** * 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); CUDAForwardProjectorWithMotion cudaForwardProjector = new CUDAForwardProjectorWithMotion(); 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, Martin Berger, Jang Hwang Choi * CONRAD is developed as an Open Source project under the GNU General Public License (GPL). */