package edu.stanford.rsl.tutorial.motion.compensation; import ij.IJ; import java.io.IOException; import java.nio.FloatBuffer; import java.nio.IntBuffer; import java.util.ArrayList; import java.util.Arrays; import com.jogamp.opencl.CLBuffer; import com.jogamp.opencl.CLCommandQueue; import com.jogamp.opencl.CLContext; import com.jogamp.opencl.CLDevice; import com.jogamp.opencl.CLImage2d; 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.apps.gui.Citeable; import edu.stanford.rsl.conrad.data.numeric.Grid2D; import edu.stanford.rsl.conrad.data.numeric.Grid3D; import edu.stanford.rsl.conrad.numerics.SimpleMatrix; import edu.stanford.rsl.conrad.opencl.OpenCLUtil; 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; /** * DO NOT USE IN LARGEVOLUMEMODE * @author Marco * */ public class OpenCLCompensatedBackProjector extends VOIBasedReconstructionFilter implements Runnable, Citeable{ /** * */ private static final long serialVersionUID = -8615490043940236889L; private ImageGridBuffer inputQueue = new ImageGridBuffer(); static int bpBlockSize[] = {32, 16}; private static boolean debug = false; /** * The OpenCL context */ private CLContext context; /** * The OpenCL program */ private CLProgram program; /** * The OpenCL device */ private CLDevice device; /** * The OpenCL kernel function binding */ private CLKernel kernelFunction; /** * The OpenCL command queue */ private CLCommandQueue commandQueue; /** * The 2D projection texture reference */ private CLImage2d<FloatBuffer> projectionTex = null; /** * The volume data that is to be reconstructed */ protected float h_volume[]; /** * The global variable of the module which stores the * view matrix. */ private CLBuffer<FloatBuffer> projectionMatrix = null; private CLBuffer<IntBuffer> volStride = null; private CLBuffer<FloatBuffer> volumePointer = null; private CLBuffer<FloatBuffer> projectionArray = null; protected ImageGridBuffer projections = new ImageGridBuffer(); 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 OpenCLCompensatedBackProjector () { super(); } @Override public void prepareForSerialization(){ super.prepareForSerialization(); projectionMatrix = null; volStride = null; volumePointer = null; projectionArray = null; projections = null; projectionsAvailable =null; projectionsDone = null; h_volume = null; initialized = false; // JOCL members program = null; device = null; kernelFunction = null; commandQueue = null; projectionTex = null; context = null; } @Override public void configure() throws Exception{ configured = true; } protected void init(){ if (!initialized) { largeVolumeMode = false; int reconDimensionX = getGeometry().getReconDimensionX(); int reconDimensionY = getGeometry().getReconDimensionY(); int reconDimensionZ = getGeometry().getReconDimensionZ(); projectionsAvailable = new ArrayList<Integer>(); projectionsDone = new ArrayList<Integer>(); // Initialize JOCL. context = OpenCLUtil.createContext(); try { // get the fastest device device = context.getMaxFlopsDevice(); // create the command queue commandQueue = device.createCommandQueue(); // initialize the program if (program==null || !program.getContext().equals(this.context)){ program = context.createProgram(OpenCLCompensatedBackProjector.class.getResourceAsStream("compensatedBackprojectCL.cl")).build(); } } catch (Exception e) { if (commandQueue != null) commandQueue.release(); if (kernelFunction != null) kernelFunction.release(); if (program != null) program.release(); // destory context if (context != null) context.release(); // TODO: handle exception e.printStackTrace(); } // check space on device: long memory = device.getMaxMemAllocSize(); long availableMemory = (memory); long requiredMemory = (long)((( ((double) reconDimensionX) * reconDimensionY * ((double) reconDimensionZ) * 4) + (((double)Configuration.getGlobalConfiguration().getGeometry().getDetectorHeight()) * Configuration.getGlobalConfiguration().getGeometry().getDetectorWidth() * 4))); if (debug) { System.out.println("Total available Memory on OpenCL card:" + availableMemory); System.out.println("Required Memory on OpenCL card:" + requiredMemory); } if (requiredMemory > availableMemory){ nSteps = (int)OpenCLUtil.iDivUp (requiredMemory, availableMemory); if (debug) System.out.println("Switching to large volume mode with nSteps = " + nSteps); largeVolumeMode = true; } if (debug) { //TODO replace /* CUdevprop prop = new CUdevprop(); JCudaDriver.cuDeviceGetProperties(prop, dev); System.out.println(prop.toFormattedString()); */ } // create the computing kernel kernelFunction = program.createCLKernel("backprojectKernel"); // create the reconstruction volume; int memorysize = reconDimensionX * reconDimensionY * reconDimensionZ * 4; if (largeVolumeMode){ subVolumeZ = OpenCLUtil.iDivUp(reconDimensionZ, nSteps); if(debug) System.out.println("SubVolumeZ: " + subVolumeZ); h_volume = new float[reconDimensionX * reconDimensionY * subVolumeZ]; memorysize = reconDimensionX * reconDimensionY * subVolumeZ * 4; if(debug)System.out.println("Memory: " + memorysize); } else { h_volume = new float[reconDimensionX * reconDimensionY * reconDimensionZ]; } // 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]; // copy volume to device volumePointer = context.createFloatBuffer(h_volume.length, Mem.WRITE_ONLY); volumePointer.getBuffer().put(h_volume); volumePointer.getBuffer().rewind(); // copy volume stride to device volStride = context.createIntBuffer(volStrideHost.length, Mem.READ_ONLY); volStride.getBuffer().put(volStrideHost); volStride.getBuffer().rewind(); commandQueue. putWriteBuffer(volumePointer, true). putWriteBuffer(volStride, true). finish(); initialized = true; } } private synchronized void unload(){ if (initialized) { if ((projectionVolume != null) && (!largeVolumeMode)) { commandQueue.putReadBuffer(volumePointer, true).finish(); volumePointer.getBuffer().rewind(); volumePointer.getBuffer().get(h_volume); volumePointer.getBuffer().rewind(); 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 commandQueue.release(); if (projectionTex != null) projectionTex.release(); if (projectionMatrix != null) projectionMatrix.release(); if (volStride != null) volStride.release(); if (projectionArray != null) projectionArray.release(); if (volumePointer != null) volumePointer.release(); kernelFunction.release(); program.release(); // destory context context.release(); commandQueue = null; projectionArray = null; projectionMatrix = null; projectionTex = null; volStride = null; volumePointer = null; kernelFunction = null; program = null; context = null; 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); } } // Obtain the global pointer to the view matrix from // the module if (projectionMatrix == null) projectionMatrix = context.createFloatBuffer(pMatFloat.length, Mem.READ_ONLY); projectionMatrix.getBuffer().put(pMatFloat); projectionMatrix.getBuffer().rewind(); commandQueue.putWriteBuffer(projectionMatrix, true).finish(); } 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 array that will contain the // projection data. projectionArray = context.createFloatBuffer(projection.getWidth()*projection.getHeight(), Mem.READ_ONLY); } // Copy the projection data to the array projectionArray.getBuffer().put(proj); projectionArray.getBuffer().rewind(); // set the texture CLImageFormat format = new CLImageFormat(ChannelOrder.INTENSITY, ChannelType.FLOAT); projectionTex = context.createImage2d(projectionArray.getBuffer(), projection.getWidth(), projection.getHeight(), format, Mem.READ_ONLY); //projectionArray.release(); } else { System.out.println("Projection was null!!"); } } @Override public String getBibtexCitation() { String bibtex = "@inproceedings{Rohkohl08-CCR,\n" + " author = {{Scherl}, H. and {Keck}, B. and {Kowarschik}, M. and {Hornegger}, J.},\n" + " title = {{Fast GPU-Based CT Reconstruction using the Common Unified Device Architecture (CUDA)}},\n" + " booktitle = {{Nuclear Science Symposium, Medical Imaging Conference 2007}},\n" + " publisher = {IEEE},\n" + " volume={6},\n" + " address = {Honolulu, HI, United States},\n" + " year = {2007}\n" + " pages= {4464--4466},\n" + "}"; return bibtex; } @Override public String getMedlineCitation() { return "Scherl H, Keck B, Kowarschik M, Hornegger J. Fast GPU-Based CT Reconstruction using the Common Unified Device Architecture (CUDA). In Nuclear Science Symposium, Medical Imaging Conference Record, IEEE, Honolulu, HI, United States, 2008 6:4464-6."; } public void waitForResult(double[] motionfield) { OpenCLRun(motionfield); } @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, float respoffset){ // load projection matrix initProjectionMatrix(projectionNumber); // load projection Grid2D projection = projections.get(projectionNumber); initProjectionData(projection); if (!largeVolumeMode) { //projections.remove(projectionNumber); } // backproject for each slice // OpenCL Grids are only two dimensional! int reconDimensionZ = dimz; double voxelSpacingX = getGeometry().getVoxelSpacingX(); double voxelSpacingY = getGeometry().getVoxelSpacingY(); double voxelSpacingZ = getGeometry().getVoxelSpacingZ(); // write kernel parameters kernelFunction.rewind(); kernelFunction .putArg(volumePointer) .putArg(respoffset) .putArg((int) lineOffset) .putArg(reconDimensionZ) .putArg((float) voxelSpacingX) .putArg((float) voxelSpacingY) .putArg((float) voxelSpacingZ) .putArg((float) offsetX) .putArg((float) offsetY) .putArg((float) offsetZ) .putArg(projectionTex) .putArg(volStride) .putArg(projectionMatrix); int[] realLocalSize = {Math.min(device.getMaxWorkGroupSize(),bpBlockSize[0]), Math.min(device.getMaxWorkGroupSize(),bpBlockSize[1])}; // rounded up to the nearest multiple of localWorkSize int[] globalWorkSize = {getGeometry().getReconDimensionX(), getGeometry().getReconDimensionY()}; // Call the OpenCL kernel, writing the results into the volume which is pointed at commandQueue .putWriteImage(projectionTex, false) .finish() .put2DRangeKernel(kernelFunction, 0, 0, globalWorkSize[0], globalWorkSize[1], realLocalSize[0], realLocalSize[1]) //.finish() //.putReadBuffer(dOut, true) .finish(); } public void OpenCLRun(double[] motionfield) { 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 OpenCL Buffer"); } else { IJ.showStatus("Backprojecting with OpenCL"); } IJ.showProgress(status); } if (!largeVolumeMode) { workOnProjectionData(motionfield); } 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 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 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); volumePointer.getBuffer().rewind(); volumePointer.getBuffer().put(h_volume); volumePointer.getBuffer().rewind(); commandQueue.putWriteBuffer(volumePointer, true).finish(); 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 OpenCL"); IJ.showProgress(currentStep/all); } //System.out.println("Current: " + p); float respoffset = (float) Math.round(motionfield[p]/voxelSpacingZ); try { projectSingleProjection(p, reconDimensionZ, respoffset); } catch (Exception e){ System.out.println("Backprojection of projection " + p + " was not successful."); e.printStackTrace(); } } // Gather volume commandQueue.putReadBuffer(volumePointer, true).finish(); volumePointer.getBuffer().rewind(); volumePointer.getBuffer().get(h_volume); volumePointer.getBuffer().rewind(); // 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 OpenCL"); 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; // exception for the case "interestedInVolume == null" and largeVolume is enabled if (interestedInVolume == null) { projectionVolume.setAtIndex(i, j, index, value); } else { 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(double[] motionfield){ if (projectionsAvailable.size() > 0){ Integer current = projectionsAvailable.get(0); projectionsAvailable.remove(0); int p = current.intValue(); projectSingleProjection(p, getGeometry().getReconDimensionZ(), (float)(motionfield[p]/getGeometry().getVoxelSpacingZ())); projectionsDone.add(current); } } private synchronized void checkProjectionData(){ if (projectionsAvailable.size() > 0){ Integer current = projectionsAvailable.get(0); projectionsAvailable.remove(current); projectionsDone.add(current); } } protected Grid3D reconstructCL(double[] motionfield){ init(); int n = inputQueue.size(); for (int i = 0; i < n; i++){ backproject(inputQueue.get(i), i); } waitForResult(motionfield); if (Configuration.getGlobalConfiguration().getUseHounsfieldScaling()) applyHounsfieldScaling(); //projectionVolume.show(); return projectionVolume; } public void loadInputQueue(ImageGridBuffer inp) throws IOException { inputQueue = inp; projections = inp; } } /* * Copyright (C) 2010-2014 Marco B�gel * CONRAD is developed as an Open Source project under the GNU General Public License (GPL). */