package edu.stanford.rsl.tutorial.cone;
import java.io.IOException;
import java.nio.FloatBuffer;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;
import java.util.concurrent.TimeUnit;
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.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;
import edu.stanford.rsl.conrad.utils.RegKeys;
public class ConeBeamBackprojector {
final boolean debug = false;
final boolean verbose = false;
private Trajectory geometry;
//image variables
private int imgSizeX;
private int imgSizeY;
private int imgSizeZ;
private Projection[] projMats;
private int maxProjs;
private double spacingX;
private double spacingY;
private double spacingZ;
private double originX;
private double originY;
private double originZ;
//cl variables
private CLContext context;
private CLDevice device;
private CLBuffer<FloatBuffer> projMatrices;
private CLCommandQueue queue;
private CLKernel kernel;
// Length of arrays to process
private int localWorkSize;
private int globalWorkSizeX;
private int globalWorkSizeY;
private CLImageFormat format;
private CLProgram program;
//normalization parameter
float normalizer;
public ConeBeamBackprojector() {
configure();
initCL();
}
public void configure(){
geometry = Configuration.getGlobalConfiguration().getGeometry();
imgSizeX = geometry.getReconDimensionX();
imgSizeY = geometry.getReconDimensionY();
imgSizeZ = geometry.getReconDimensionZ();
projMats = geometry.getProjectionMatrices();
maxProjs = geometry.getProjectionStackSize();
spacingX = geometry.getVoxelSpacingX();
spacingY = geometry.getVoxelSpacingY();
spacingZ = geometry.getVoxelSpacingZ();
originX = -geometry.getOriginX();
originY = -geometry.getOriginY();
originZ = -geometry.getOriginZ();
normalizer = (float) (geometry.getSourceToDetectorDistance()*geometry.getSourceToAxisDistance()*Math.PI / maxProjs);
}
private void initCL(){
context = OpenCLUtil.getStaticContext();
device = context.getMaxFlopsDevice();
queue = device.createCommandQueue();
// load sources, create and build program
program = null;
try {
program = context.createProgram(this.getClass().getResourceAsStream("ConeBeamBackProjector.cl")).build();
} catch (IOException e) {
// TODO Auto-generated catch block
e.printStackTrace();
System.exit(-1);
}
kernel = program.createCLKernel("backProjectPixelDrivenCL");
// create image from input grid
format = new CLImageFormat(ChannelOrder.INTENSITY, ChannelType.FLOAT);
localWorkSize = Math.min(device.getMaxWorkGroupSize(), 16);
globalWorkSizeX = OpenCLUtil.roundUp(localWorkSize, imgSizeX);
globalWorkSizeY = OpenCLUtil.roundUp(localWorkSize, imgSizeY);
projMatrices = context.createFloatBuffer(maxProjs*3*4, Mem.READ_ONLY);
for(int p = 0; p < maxProjs; p++) {
for(int row = 0; row < 3; row++) {
for(int col = 0; col < 4; col++) {
projMatrices.getBuffer().put((float)projMats[p].computeP().getElement(row, col));
}
}
}
projMatrices.getBuffer().rewind();
queue.putWriteBuffer(projMatrices, true).finish();
}
public void unload(){
if(program != null && !program.isReleased())
program.release();
if(projMatrices != null && !projMatrices.isReleased())
projMatrices.release();
}
public Grid3D backprojectPixelDriven(Grid2D sino, int projIdx) {
configure();
if(projIdx >= maxProjs || 0 > projIdx){
System.err.println("ConeBeamBackprojector: Invalid projection index");
return null;
}
Grid3D grid = new Grid3D(imgSizeX,imgSizeY,imgSizeZ);
grid.setOrigin(-originX, -originY, -originZ);
grid.setSpacing(spacingX, spacingY, spacingZ);
for(int x = 0; x < imgSizeX; x++) {
double xTrans = x*spacingX-originX;
for(int y = 0; y < imgSizeY ; y++) {
double yTrans = y*spacingY-originY;
for(int z = 0; z < imgSizeZ; 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);
float val = (float) (InterpolationOperators.interpolateLinear(sino, coordU, coordV)/(point2d.getElement(2)*point2d.getElement(2))); //
//if(Float.isInfinite(val) || Float.isNaN(val))
// val = 0;
grid.setAtIndex(x, y, z, val); // "set" instead of "add", because #proj==1
//}
}
}
}
return grid;
}
public Grid3D backprojectPixelDriven(final Grid3D sino) {
configure();
final Grid3D grid = new Grid3D(imgSizeX,imgSizeY,imgSizeZ);
int nThreads = Integer.valueOf(Configuration.getGlobalConfiguration().getRegistryEntry(RegKeys.MAX_THREADS));
// TODO Error-Checking for thread number
ExecutorService executorService = Executors.newFixedThreadPool(nThreads);
for(int i = 0; i < imgSizeX; i++) {
final int x = i;
executorService.execute(new Runnable(){
@Override
public void run(){
System.out.println("Working on slice "+String.valueOf(x+1)+" of "+imgSizeX);
double xTrans = x*spacingX-originX;
for(int y = 0; y < imgSizeY ; y++) {
double yTrans = y*spacingY-originY;
for(int z = 0; z < imgSizeZ; 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);
/*
// TEST // TODO
if(0==p)
System.out.println("Sino - Min: " + sino.getGridOperator().min(sino) + " Max: " + sino.getGridOperator().max(sino));
int uplusOne = Math.min((int)coordU+1, geo.getDetectorWidth()-1);
int vplusOne = Math.min((int)coordV+1, geo.getDetectorHeight()-1);
float[] values = new float[]{sino.getAtIndex((int)coordU, (int)coordV, p), sino.getAtIndex(uplusOne, (int)coordV, p),
sino.getAtIndex((int)coordU, vplusOne, p), sino.getAtIndex(uplusOne, vplusOne, p)};
float sumVal = values[0] + values[1] + values[2] + values[3];
if (sumVal != 0)
System.out.println("NOT NULL");
// /TEST
*/
float val = (float) (InterpolationOperators.interpolateLinear(sino, p, coordU, coordV)/(point2d.getElement(2)*point2d.getElement(2)));
/*
// TEST // TODO
if(sumVal/4 != val || 0 != val)
System.out.println("[" + x + ", " + y + ", " +z + "] = " + val + " = interp(" + p + ", " + coordU + ", " + coordV + "), " + sumVal/4);
// /TEST
*/
grid.addAtIndex(x, y, z, val);
}
}
}
}});
}
executorService.shutdown();
try {
executorService.awaitTermination(Long.MAX_VALUE, TimeUnit.NANOSECONDS);
} catch (InterruptedException e) {
System.out.println("Exception waiting for thread termination: "+e);
}
return grid;
}
public void backprojectPixelDrivenCL(OpenCLGrid3D volume, OpenCLGrid2D[] sino) {
for(int p = 0; p < maxProjs; p++) {
CLImage2d<FloatBuffer> sinoGrid = context.createImage2d(sino[p].getDelegate().getCLBuffer().getBuffer(), sino[p].getSize()[0], sino[p].getSize()[1],format,Mem.READ_ONLY);
kernel.putArg(sinoGrid)
.putArg(volume.getDelegate().getCLBuffer())
.putArg(projMatrices)
.putArg(p)
.putArg(imgSizeX).putArg(imgSizeY).putArg(imgSizeZ)
.putArg((float)originX).putArg((float)originY).putArg((float)originZ)
.putArg((float)spacingX).putArg((float)spacingY).putArg((float)spacingZ)
.putArg(normalizer);
queue
.putCopyBufferToImage(sino[p].getDelegate().getCLBuffer(), sinoGrid).finish()
.put2DRangeKernel(kernel, 0, 0, globalWorkSizeX, globalWorkSizeY,localWorkSize, localWorkSize)
.finish();
kernel.rewind();
sinoGrid.release();
}
volume.getDelegate().notifyDeviceChange();
}
public Grid3D backprojectPixelDrivenCL(Grid3D sino) {
configure();
OpenCLGrid2D [] sinoCL = new OpenCLGrid2D[sino.getSize()[2]];
for (int i=0; i < sinoCL.length; i++){ sinoCL[i] = new OpenCLGrid2D(sino.getSubGrid(i)); sinoCL[i].getDelegate().prepareForDeviceOperation();}
Grid3D grid = new Grid3D(imgSizeX,imgSizeY,imgSizeZ);
OpenCLGrid3D gridCL = new OpenCLGrid3D(grid);
gridCL.getDelegate().prepareForDeviceOperation();
backprojectPixelDrivenCL(gridCL, sinoCL);
gridCL.setOrigin(-originX, -originY, -originZ);
gridCL.setSpacing(spacingX, spacingY, spacingZ);
for (int i=0; i < sinoCL.length; i++) sinoCL[i].release();
grid = new Grid3D(gridCL);
gridCL.release();
unload();
return grid;
}
public Grid3D backprojectPixelDrivenCL(Grid2D sino , int projIdx) {
configure();
OpenCLGrid2D sinoCL = new OpenCLGrid2D(sino);
Grid3D grid = new Grid3D(imgSizeX,imgSizeY,imgSizeZ);
OpenCLGrid3D gridCL = new OpenCLGrid3D(grid);
gridCL.getDelegate().prepareForDeviceOperation();
backprojectPixelDrivenCL(gridCL, sinoCL, projIdx);
gridCL.setOrigin(-originX, -originY, -originZ);
gridCL.setSpacing(spacingX, spacingY, spacingZ);
sinoCL.release();
grid = new Grid3D(gridCL);
gridCL.release();
unload();
return grid;
}
public void backprojectPixelDrivenCL(OpenCLGrid3D volume, OpenCLGrid2D sino, int projIdx) {
//TODO MOEGLICHE FEHLERQUELLE
CLImage2d<FloatBuffer> sinoGrid = context.createImage2d(sino.getDelegate().getCLBuffer().getBuffer(), sino.getSize()[0], sino.getSize()[1],format,Mem.READ_ONLY);
kernel.putArg(sinoGrid)
.putArg(volume.getDelegate().getCLBuffer())
.putArg(projMatrices)
.putArg(projIdx)
.putArg(imgSizeX).putArg(imgSizeY).putArg(imgSizeZ)
.putArg((float)originX).putArg((float)originY).putArg((float)originZ)
.putArg((float)spacingX).putArg((float)spacingY).putArg((float)spacingZ)
.putArg(normalizer);
queue
.putCopyBufferToImage(sino.getDelegate().getCLBuffer(), sinoGrid).finish()
.put2DRangeKernel(kernel, 0, 0, globalWorkSizeX, globalWorkSizeY,localWorkSize, localWorkSize)
.finish();
kernel.rewind();
volume.getDelegate().notifyDeviceChange();
sinoGrid.release();
}
public void fastBackprojectPixelDrivenCL(OpenCLGrid2D sinoCL, OpenCLGrid3D gridCL, int projIdx) {
configure();
gridCL.getDelegate().prepareForDeviceOperation();
sinoCL.getDelegate().prepareForDeviceOperation();
backprojectPixelDrivenCL(gridCL, sinoCL, projIdx);
}
public void fastBackprojectPixelDrivenCL(OpenCLGrid3D sinoCL, OpenCLGrid3D gridCL) {
configure();
gridCL.getDelegate().prepareForDeviceOperation();
sinoCL.getDelegate().prepareForDeviceOperation();
OpenCLGrid2D[] sinoBuf = new OpenCLGrid2D[sinoCL.getSize()[2]];
for(int i = 0; i< sinoCL.getSize()[2];i++){
sinoBuf[i] = new OpenCLGrid2D(sinoCL.getSubGrid(i));
sinoBuf[i].getDelegate().prepareForDeviceOperation();
}
backprojectPixelDrivenCL(gridCL,sinoBuf);
for(int i = 0; i< sinoCL.getSize()[2];i++)
sinoBuf[i].release();
unload();
}
}
/*
* Copyright (C) 2010-2014 Andreas Maier
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/