/*
* Copyright (C) 2014 Michael Manhart
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/
package edu.stanford.rsl.conrad.opencl;
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.apps.gui.Citeable;
import edu.stanford.rsl.conrad.geometry.trajectories.Trajectory;
import edu.stanford.rsl.conrad.numerics.SimpleMatrix;
import edu.stanford.rsl.conrad.utils.CONRAD;
import edu.stanford.rsl.conrad.utils.Configuration;
public class OpenCLEdgeForwardProjector implements Citeable {
/**
* forward projector configuration
*/
public static enum EdgeDirection {EDGES_UV, EDGES_U, EDGES_V};
private boolean isInitialized = false;
private boolean isConfigured = false;
private float [] voxelSize;
private int [] volumeSize;
private float [] volumeEdgeMinPoint;
private float [] volumeEdgeMaxPoint;
private Trajectory g_fwd;
private Trajectory g_bwd;
/**
* OpenCL
*/
static int bpBlockSize[] = {16, 16};
private CLContext cl_context;
private CLProgram cl_program;
private CLDevice cl_device;
private CLKernel cl_kernelFunctionEdgeProjector;
private CLKernel cl_kernelFunctionEdgeProjectorU;
private CLKernel cl_kernelFunctionEdgeProjectorV;
private CLCommandQueue cl_commandQueue;
private CLImage3d<FloatBuffer> cl_VolumeTexture3D = null;
private CLBuffer<FloatBuffer> cl_VolumeEdgeMaxPoint = null;
private CLBuffer<FloatBuffer> cl_VolumeEdgeMinPoint = null;
private CLBuffer<FloatBuffer> cl_VoxelElementSize = null;
private CLBuffer<FloatBuffer> cl_InvARmatrixFwd = null;
private CLBuffer<FloatBuffer> cl_SrcPointFwd = null;
private CLBuffer<FloatBuffer> cl_InvARmatrixBwd = null;
private CLBuffer<FloatBuffer> cl_SrcPointBwd = null;
private CLBuffer<FloatBuffer> cl_Projection = null;
/**
* configuration of the projection and volume geometries
*/
public void configure(Trajectory g_fwd, Trajectory g_bwd, int[] volumeSize, float[] voxelSize) {
this.voxelSize = new float [3];
this.volumeSize = new int [4];
this.voxelSize[0] = voxelSize[0];
this.voxelSize[1] = voxelSize[1];
this.voxelSize[2] = voxelSize[2];
this.volumeSize[0] = volumeSize[0];
this.volumeSize[1] = volumeSize[1];
this.volumeSize[2] = volumeSize[2];
this.volumeSize[3] = volumeSize[0]*volumeSize[1]*volumeSize[2];
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);
}
this.g_fwd = g_fwd;
this.g_bwd = g_bwd;
isConfigured = true;
}
/**
* Initiates communication with the graphics card.
*/
private void initialize() {
if(!isConfigured) {
System.err.println("OpenCLForwardProjectorDynamicVolume.initialize(): Need to configure before initalization!");
return;
}
if (!isInitialized) {
// Initialize JOCL.
cl_context = OpenCLUtil.createContext();
try {
// get the fastest device
cl_device = cl_context.getMaxFlopsDevice();
// create the command queue
cl_commandQueue = cl_device.createCommandQueue();
// initialize the program
if (cl_program == null || !cl_program.getContext().equals(this.cl_context)){
cl_program = cl_context.createProgram(TestOpenCL.class.getResourceAsStream("OpenCLEdgeForwardProjector.cl")).build();
}
// create the computing kernel
cl_kernelFunctionEdgeProjector = cl_program.createCLKernel("projectKernel");
cl_kernelFunctionEdgeProjectorU = cl_program.createCLKernel("projectKernel_UEdges");
cl_kernelFunctionEdgeProjectorV = cl_program.createCLKernel("projectKernel_VEdges");
// volume properties
cl_VoxelElementSize = cl_context.createFloatBuffer(voxelSize.length, Mem.READ_ONLY);
cl_VoxelElementSize.getBuffer().put(voxelSize);
cl_VoxelElementSize.getBuffer().rewind();
cl_VolumeEdgeMinPoint = cl_context.createFloatBuffer(volumeEdgeMinPoint.length, Mem.READ_ONLY);
cl_VolumeEdgeMinPoint.getBuffer().put(volumeEdgeMinPoint);
cl_VolumeEdgeMinPoint.getBuffer().rewind();
cl_VolumeEdgeMaxPoint = cl_context.createFloatBuffer(volumeEdgeMaxPoint.length, Mem.READ_ONLY);
cl_VolumeEdgeMaxPoint.getBuffer().put(volumeEdgeMaxPoint);
cl_VolumeEdgeMaxPoint.getBuffer().rewind();
// projection memory
cl_Projection = cl_context.createFloatBuffer(g_fwd.getDetectorWidth()*g_fwd.getDetectorHeight(), Mem.WRITE_ONLY);
cl_commandQueue
.putWriteBuffer(cl_VoxelElementSize,true)
.putWriteBuffer(cl_VolumeEdgeMinPoint, true)
.putWriteBuffer(cl_VolumeEdgeMaxPoint, true)
.finish();
} catch (IOException e) {
// TODO: handle exception
e.printStackTrace();
}
cl_InvARmatrixFwd = cl_context.createFloatBuffer(9*g_fwd.getNumProjectionMatrices(), Mem.READ_ONLY);
cl_SrcPointFwd = cl_context.createFloatBuffer(3*g_fwd.getNumProjectionMatrices(), Mem.READ_ONLY);
prepareProjectionMatrices(g_fwd, cl_InvARmatrixFwd, cl_SrcPointFwd);
cl_InvARmatrixBwd = cl_context.createFloatBuffer(9*g_bwd.getNumProjectionMatrices(), Mem.READ_ONLY);
cl_SrcPointBwd = cl_context.createFloatBuffer(3*g_bwd.getNumProjectionMatrices(), Mem.READ_ONLY);
prepareProjectionMatrices(g_bwd, cl_InvARmatrixBwd, cl_SrcPointBwd);
isInitialized = true;
}
}
/**
* Compute inverted projection matrices and source positions for all projections and write in GPU memory
* @param projectionNumber
*/
private void prepareProjectionMatrices(Trajectory g, CLBuffer<FloatBuffer> cl_InvARmatrix, CLBuffer<FloatBuffer> cl_SrcPoint){
float [] cann = new float[3*4];
float [] invAR = new float[3*3];
float [] srcP = new float[3];
for (int i=0; i < g.getNumProjectionMatrices(); ++i){
SimpleMatrix projMat = g.getProjectionMatrix(i).computeP();
double [][] mat = new double [3][4];
projMat.copyTo(mat);
computeCanonicalProjectionMatrix(cann, invAR, srcP, new Jama.Matrix(mat));
cl_InvARmatrix.getBuffer().put(invAR);
cl_SrcPoint.getBuffer().put(srcP);
}
cl_InvARmatrix.getBuffer().rewind();
cl_SrcPoint.getBuffer().rewind();
cl_commandQueue
.putWriteBuffer(cl_SrcPoint, true)
.putWriteBuffer(cl_InvARmatrix, true)
.finish();
}
public void setVolume(float[] volumeBuffer) {
if(!isConfigured) {
System.err.println("OpenCLForwardProjectorDynamicVolume.setVolume(): Need to configure before forward projection!");
return;
}
if(volumeBuffer.length != volumeSize[3])
{
System.err.println("OpenCLForwardProjectorDynamicVolume.setVolume(): Invalid volume buffer size!");
return;
}
initialize();
if(cl_VolumeTexture3D != null) {
cl_VolumeTexture3D.release();
cl_VolumeTexture3D = null;
}
// create 3d texture for projected volumes
CLBuffer<FloatBuffer> hvolumeBuffer = cl_context.createFloatBuffer(volumeBuffer.length, Mem.READ_ONLY);
hvolumeBuffer.getBuffer().put(volumeBuffer);
hvolumeBuffer.getBuffer().rewind();
CLImageFormat format = new CLImageFormat(ChannelOrder.INTENSITY, ChannelType.FLOAT);
cl_VolumeTexture3D = cl_context.createImage3d(hvolumeBuffer.getBuffer(), volumeSize[0], volumeSize[1], volumeSize[2], format, Mem.READ_ONLY);
cl_commandQueue
.putWriteImage(cl_VolumeTexture3D, true)
.finish();
hvolumeBuffer.release();
}
public void applyForwardProjection(int projectionID, boolean fwd, EdgeDirection edgeDirection, float[] projection_out) {
if(!isConfigured) {
System.err.println("OpenCLForwardProjectorDynamicVolume.applyForwardProjection(): Need to configure before forward projection!");
return;
}
if(cl_VolumeTexture3D == null) {
System.err.println("OpenCLForwardProjectorDynamicVolume.applyForwardProjection(): Need to set volume before forward projection!");
return;
}
Trajectory g = fwd?g_fwd:g_bwd;
int projSize = g.getDetectorWidth()*g.getDetectorHeight();
if(projection_out.length != projSize) {
System.err.println("OpenCLForwardProjectorDynamicVolume.applyForwardProjection(): Wrong projection buffer size!");
return;
}
CLKernel cl_kernel = null;
switch(edgeDirection) {
case EDGES_UV:
cl_kernel = cl_kernelFunctionEdgeProjector;
break;
case EDGES_U:
cl_kernel = cl_kernelFunctionEdgeProjectorU;
break;
case EDGES_V:
cl_kernel = cl_kernelFunctionEdgeProjectorV;
break;
}
// write kernel parameters
cl_kernel.rewind();
cl_kernel
.putArg(cl_Projection)
.putArg(g.getDetectorWidth())
.putArg(g.getDetectorHeight())
.putArg(0.5f)
.putArg(cl_VolumeTexture3D)
.putArg(cl_VoxelElementSize)
.putArg(cl_VolumeEdgeMinPoint)
.putArg(cl_VolumeEdgeMaxPoint)
.putArg(fwd?cl_SrcPointFwd:cl_SrcPointBwd)
.putArg(fwd?cl_InvARmatrixFwd:cl_InvARmatrixBwd)
.putArg(projectionID);
int[] realLocalSize = new int[2];
realLocalSize[0] = Math.min(cl_device.getMaxWorkGroupSize(),bpBlockSize[0]);
realLocalSize[1] = Math.max(1, Math.min(cl_device.getMaxWorkGroupSize()/realLocalSize[0], bpBlockSize[1]));
// rounded up to the nearest multiple of localWorkSize
int[] globalWorkSize = {g.getDetectorWidth(), g.getDetectorHeight()};
if ((globalWorkSize[0] % realLocalSize[0] ) != 0){
globalWorkSize[0] = ((globalWorkSize[0] / realLocalSize[0]) + 1) * realLocalSize[0];
}
if ((globalWorkSize[1] % realLocalSize[1] ) != 0){
globalWorkSize[1] = ((globalWorkSize[1] / realLocalSize[1]) + 1) * realLocalSize[1];
}
// add kernel function to the queue
cl_commandQueue
.put2DRangeKernel(cl_kernel, 0, 0, globalWorkSize[0], globalWorkSize[1], realLocalSize[0], realLocalSize[1])
.finish()
.putReadBuffer(cl_Projection, true)
.finish();
// copy result from device to host
cl_Projection.getBuffer().rewind();
cl_Projection.getBuffer().get(projection_out);
cl_Projection.getBuffer().rewind();
}
/**
*
*
* 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()};
du[0] = 1;
dv[1] = 1;
// 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);
// 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);
}
/**
* release all CL related objects and free memory
*/
public void release(){
if (cl_commandQueue != null)
cl_commandQueue.release();
if (cl_VolumeTexture3D != null) {
cl_VolumeTexture3D.release();
cl_VolumeTexture3D = null;
}
if (cl_VolumeEdgeMaxPoint != null)
cl_VolumeEdgeMaxPoint.release();
if (cl_VolumeEdgeMinPoint != null)
cl_VolumeEdgeMinPoint.release();
if (cl_VoxelElementSize != null)
cl_VoxelElementSize.release();
if (cl_InvARmatrixFwd != null)
cl_InvARmatrixFwd.release();
if (cl_InvARmatrixBwd != null)
cl_InvARmatrixBwd.release();
if (cl_SrcPointFwd != null)
cl_SrcPointFwd.release();
if (cl_SrcPointBwd != null)
cl_SrcPointBwd.release();
if (cl_Projection != null)
cl_Projection.release();
if (cl_kernelFunctionEdgeProjector != null)
cl_kernelFunctionEdgeProjector.release();
if (cl_kernelFunctionEdgeProjectorU != null)
cl_kernelFunctionEdgeProjectorU.release();
if (cl_kernelFunctionEdgeProjectorV != null)
cl_kernelFunctionEdgeProjectorV.release();
if (cl_program != null)
cl_program.release();
if (cl_context != null)
cl_context.release();
isConfigured = false;
isInitialized = false;
}
@Override
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;
}
@Override
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.";
}
}