package edu.stanford.rsl.tutorial.fan;
import java.io.IOException;
import java.nio.FloatBuffer;
import java.util.ArrayList;
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.CLMemory.Mem;
import com.jogamp.opencl.CLProgram;
import com.jogamp.opencl.CLImageFormat.ChannelOrder;
import com.jogamp.opencl.CLImageFormat.ChannelType;
import edu.stanford.rsl.conrad.data.numeric.Grid1D;
import edu.stanford.rsl.conrad.data.numeric.Grid2D;
import edu.stanford.rsl.conrad.data.numeric.NumericGrid;
import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators;
import edu.stanford.rsl.conrad.data.numeric.InterpolationOperators;
import edu.stanford.rsl.conrad.data.numeric.opencl.OpenCLGrid1D;
import edu.stanford.rsl.conrad.data.numeric.opencl.OpenCLGrid2D;
import edu.stanford.rsl.conrad.geometry.shapes.simple.Box;
import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND;
import edu.stanford.rsl.conrad.geometry.shapes.simple.StraightLine;
import edu.stanford.rsl.conrad.geometry.transforms.Transform;
import edu.stanford.rsl.conrad.geometry.transforms.Translation;
import edu.stanford.rsl.conrad.numerics.SimpleVector;
import edu.stanford.rsl.conrad.opencl.OpenCLUtil;
public class FanBeamProjector2D{
final double samplingRate = 3.d;
private double focalLength, maxBeta, deltaBeta, maxT, deltaT;
int maxTIndex, maxBetaIndex;
/**
* Creates a new instance of a fan-beam projector
*
* @param focalLength the focal length
* @param maxBeta the maximal rotation angle
* @param deltaBeta the step size between source positions
* @param maxT the length of the detector array
* @param deltaT the size of one detector element
*/
public FanBeamProjector2D(double focalLength,
double maxBeta, double deltaBeta,
double maxT, double deltaT) {
this.focalLength = focalLength;
this.maxBeta = maxBeta;
this.maxT = maxT;
this.deltaBeta = deltaBeta;
this.deltaT = deltaT;
//this.maxBetaIndex = (int) (maxBeta / deltaBeta + 1);
this.maxBetaIndex = (int) (maxBeta / deltaBeta );
this.maxTIndex = (int) (maxT / deltaT);
}
public Grid2D projectRayDriven(Grid2D grid) {
Grid2D sino = new Grid2D(maxTIndex, maxBetaIndex);
sino.setSpacing(deltaT, deltaBeta);
// create translation to the grid origin
Translation trans = new Translation(-grid.getSize()[0] / 2.0,
-grid.getSize()[1] / 2.0, -1);
// build the inverse translation
Transform inverse = trans.inverse();
// set up image bounding box and translate to origin
Box b = new Box(grid.getSize()[0], grid.getSize()[1], 2);
b.applyTransform(trans);
// iterate over the rotation angle
for (int i = 0; i < maxBetaIndex; i++) {
// compute the current rotation angle and its sine and cosine
double beta = deltaBeta * i;
double cosBeta = Math.cos(beta);
double sinBeta = Math.sin(beta);
// System.out.println(beta / Math.PI * 180);
// compute source position
PointND a = new PointND(focalLength * cosBeta, focalLength
* sinBeta, 0.d);
// compute end point of detector
PointND p0 = new PointND(-maxT / 2.f * sinBeta, maxT / 2.f
* cosBeta, 0.d);
// create an unit vector that points along the detector
SimpleVector dirDetector = p0.getAbstractVector().multipliedBy(-1);
dirDetector.normalizeL2();
// iterate over the detector elements
for (int t = 0; t < maxTIndex; t++) {
// calculate current bin position
// the detector elements' position are centered
double stepsDirection = 0.5f * deltaT + t * deltaT;
PointND p = new PointND(p0);
p.getAbstractVector().add(dirDetector.multipliedBy(stepsDirection));
// create a straight line between detector bin and source
StraightLine line = new StraightLine(a, p);
// find the line's intersection with the box
ArrayList<PointND> points = b.intersect(line);
// if we have two intersections build the integral
// otherwise continue with the next bin
if (2 != points.size()) {
if (points.size() == 0) {
line.getDirection().multiplyBy(-1.d);
points = b.intersect(line);
if (points.size() == 0)
continue;
} else {
continue; // last possibility:
// a) it is only one intersection point (exactly one of the boundary vertices) or
// b) it are infinitely many intersection points (along one of the box boundaries).
// c) our code is wrong
}
}
// Extract intersections
PointND start = points.get(0);
PointND end = points.get(1);
// get the normalized increment
SimpleVector increment = new SimpleVector(end.getAbstractVector());
increment.subtract(start.getAbstractVector());
double distance = increment.normL2();
increment.divideBy(distance * samplingRate);
double sum = .0;
start = inverse.transform(start);
double incrementLength = increment.normL2();
// compute the integral along the line
for (double tLine = 0.0; tLine < distance * samplingRate; ++tLine) {
PointND current = new PointND(start);
current.getAbstractVector().add(increment.multipliedBy(tLine));
if (grid.getSize()[0] <= current.get(0) + 1
|| grid.getSize()[1] <= current.get(1) + 1
|| current.get(0) < 0 || current.get(1) < 0)
continue;
sum += InterpolationOperators.interpolateLinear(grid,
current.get(0), current.get(1));
}
// normalize by the number of interpolation points
sum /= samplingRate;
// write integral value into the sinogram.
sino.setAtIndex(t, i, (float) sum);
}
}
return sino;
}
public Grid1D projectRayDriven1D(Grid2D grid, int Betaindex) {
Grid1D sino=new Grid1D(this.maxTIndex);
sino.setSpacing(this.deltaT);
// create translation to the grid origin
Translation trans = new Translation(-grid.getSize()[0] / 2.0,-grid.getSize()[1] / 2.0, -1);
// build the inverse translation
Transform inverse = trans.inverse();
// set up image bounding box and translate to origin
Box b = new Box(grid.getSize()[0], grid.getSize()[1], 2);
b.applyTransform(trans);
// compute the current rotation angle and its sine and cosine
double beta = deltaBeta * Betaindex;
double cosBeta = Math.cos(beta);
double sinBeta = Math.sin(beta);
// compute source position
PointND a = new PointND(focalLength * cosBeta, focalLength
* sinBeta, 0.d);
// compute end point of detector
PointND p0 = new PointND(-maxT / 2.f * sinBeta, maxT / 2.f
* cosBeta, 0.d);
// create an unit vector that points along the detector
SimpleVector dirDetector = p0.getAbstractVector().multipliedBy(-1);
dirDetector.normalizeL2();
// iterate over the detector elements
for (int t = 0; t < maxTIndex; t++) {
// calculate current bin position
// the detector elements' position are centered
double stepsDirection = 0.5f * deltaT + t * deltaT;
PointND p = new PointND(p0);
p.getAbstractVector().add(dirDetector.multipliedBy(stepsDirection));
// create a straight line between detector bin and source
StraightLine line = new StraightLine(a, p);
// find the line's intersection with the box
ArrayList<PointND> points = b.intersect(line);
// if we have two intersections build the integral
// otherwise continue with the next bin
if (2 != points.size()) {
if (points.size() == 0) {
line.getDirection().multiplyBy(-1.d);
points = b.intersect(line);
if (points.size() == 0)
continue;
} else {
continue; // last possibility:
// a) it is only one intersection point (exactly one of the boundary vertices) or
// b) it are infinitely many intersection points (along one of the box boundaries).
// c) our code is wrong
}
}
// Extract intersections
PointND start = points.get(0);
PointND end = points.get(1);
// get the normalized increment
SimpleVector increment = new SimpleVector(end.getAbstractVector());
increment.subtract(start.getAbstractVector());
double distance = increment.normL2();
increment.divideBy(distance * samplingRate);
double sum = .0;
start = inverse.transform(start);
// compute the integral along the line
for (double tLine = 0.0; tLine < distance * samplingRate; ++tLine) {
PointND current = new PointND(start);
current.getAbstractVector().add(increment.multipliedBy(tLine));
if (grid.getSize()[0] <= current.get(0) + 1
|| grid.getSize()[1] <= current.get(1) + 1
|| current.get(0) < 0 || current.get(1) < 0)
continue;
sum += InterpolationOperators.interpolateLinear(grid,
current.get(0), current.get(1));
}
// normalize by the number of interpolation points
sum /= samplingRate;
// write integral value into the sinogram.
sino.setAtIndex(t, (float) sum);
}
return sino;
}
public Grid2D projectRayDrivenCL(Grid2D grid) {
boolean debug = false;
// create context
CLContext context = OpenCLUtil.createContext();
if (debug)
System.out.println("Context: " + context);
//show OpenCL devices in System
CLDevice[] devices = context.getDevices();
if (debug){
for (CLDevice dev: devices)
System.out.println(dev);
}
// select device
CLDevice device = context.getMaxFlopsDevice();
if (debug)
System.out.println("Device: " + device);
int imageSize = grid.getSize()[0] * grid.getSize()[1];
// Length of arrays to process
int localWorkSize = Math.min(device.getMaxWorkGroupSize(), 8); // Local work size dimensions
int globalWorkSizeT = OpenCLUtil.roundUp(localWorkSize, maxTIndex); // rounded up to the nearest multiple of localWorkSize
int globalWorkSizeBeta = OpenCLUtil.roundUp(localWorkSize, maxBetaIndex); // rounded up to the nearest multiple of localWorkSize
// load sources, create and build program
CLProgram program = null;
try {
program = context.createProgram(this.getClass().getResourceAsStream("FanBeamProjector.cl"))
.build();
} catch (IOException e) {
// TODO Auto-generated catch block
e.printStackTrace();
System.exit(-1);
}
// create image from input grid
CLImageFormat format = new CLImageFormat(ChannelOrder.INTENSITY, ChannelType.FLOAT);
CLBuffer<FloatBuffer> imageBuffer = context.createFloatBuffer(imageSize, Mem.READ_ONLY);
for (int i=0;i<grid.getBuffer().length;++i){
imageBuffer.getBuffer().put(grid.getBuffer()[i]);
}
imageBuffer.getBuffer().rewind();
CLImage2d<FloatBuffer> imageGrid = context.createImage2d(
imageBuffer.getBuffer(), grid.getSize()[0], grid.getSize()[1],
format,Mem.READ_ONLY);
imageBuffer.release();
// create memory for sinogram
CLBuffer<FloatBuffer> sinogram = context.createFloatBuffer(maxTIndex * maxBetaIndex, Mem.WRITE_ONLY);
// copy params
CLKernel kernel = program.createCLKernel("projectRayDriven2DCL");
kernel.putArg(imageGrid).putArg(sinogram)
.putArg((float)maxT).putArg((float)deltaT)
.putArg((float)maxBeta).putArg((float)deltaBeta)
.putArg((float)focalLength).putArg(maxTIndex).putArg(maxBetaIndex); // TODO: Spacing :)
// createCommandQueue
CLCommandQueue queue = device.createCommandQueue();
queue
.putWriteImage(imageGrid, true)
.finish()
.put2DRangeKernel(kernel, 0, 0, globalWorkSizeBeta, globalWorkSizeT,
localWorkSize, localWorkSize).putBarrier()
.putReadBuffer(sinogram, true)
.finish();
// write sinogram back to grid2D
Grid2D sino = new Grid2D(maxTIndex,maxBetaIndex);
sino.setSpacing(deltaT, deltaBeta);
sinogram.getBuffer().rewind();
for (int i = 0; i < sino.getBuffer().length; ++i) {
sino.getBuffer()[i] = sinogram.getBuffer().get();
}
// clean up
queue.release();
imageGrid.release();
sinogram.release();
kernel.release();
program.release();
context.release();
return sino;
}
public void fastProjectRayDrivenCL(OpenCLGrid2D sinoCL, OpenCLGrid2D gridCL) {
sinoCL.getDelegate().prepareForDeviceOperation();
gridCL.getDelegate().prepareForDeviceOperation();
CLContext context = OpenCLUtil.getStaticContext();
CLDevice device = context.getMaxFlopsDevice();
int localWorkSize = Math.min(device.getMaxWorkGroupSize(), 8); // Local work size dimensions
int globalWorkSizeT = OpenCLUtil.roundUp(localWorkSize, maxTIndex); // rounded up to the nearest multiple of localWorkSize
int globalWorkSizeBeta = OpenCLUtil.roundUp(localWorkSize, maxBetaIndex); // rounded up to the nearest multiple of localWorkSize
// load sources, create and build program
CLProgram program = null;
try {
program = context.createProgram(this.getClass().getResourceAsStream("FanBeamProjector.cl"))
.build();
} catch (IOException e) {
// TODO Auto-generated catch block
e.printStackTrace();
System.exit(-1);
}
// create image from input grid
CLImageFormat format = new CLImageFormat(ChannelOrder.INTENSITY, ChannelType.FLOAT);
CLImage2d<FloatBuffer> imageGrid = context.createImage2d(gridCL.getDelegate().getCLBuffer().getBuffer(), gridCL.getSize()[0], gridCL.getSize()[1],format,Mem.READ_ONLY);
// create memory for sinogram
CLBuffer<FloatBuffer> sinogram = sinoCL.getDelegate().getCLBuffer();
// copy params
CLKernel kernel = program.createCLKernel("projectRayDriven2DCL");
kernel.putArg(imageGrid).putArg(sinogram)
.putArg((float)maxT).putArg((float)deltaT)
.putArg((float)maxBeta).putArg((float)deltaBeta)
.putArg((float)focalLength).putArg(maxTIndex).putArg(maxBetaIndex); // TODO: Spacing :)
// createCommandQueue
CLCommandQueue queue = device.createCommandQueue();
queue
.putCopyBufferToImage(gridCL.getDelegate().getCLBuffer(),imageGrid).finish()
.put2DRangeKernel(kernel, 0, 0, globalWorkSizeBeta, globalWorkSizeT,localWorkSize, localWorkSize).putBarrier()
.finish();
kernel.rewind();
sinoCL.getDelegate().notifyDeviceChange();
sinoCL.setSpacing(deltaT, deltaBeta);
// clean up
queue.release();
kernel.release();
program.release();
}
public Grid1D projectRayDriven1DCL(Grid2D grid,int index) {
boolean debug = false;
// create context
CLContext context = OpenCLUtil.createContext();
if (debug)
System.out.println("Context: " + context);
//show OpenCL devices in System
CLDevice[] devices = context.getDevices();
if (debug){
for (CLDevice dev: devices)
System.out.println(dev);
}
// select device
CLDevice device = context.getMaxFlopsDevice();
if (debug)
System.out.println("Device: " + device);
int imageSize = grid.getSize()[0] * grid.getSize()[1];
// Length of arrays to process
int localWorkSize = Math.min(device.getMaxWorkGroupSize(), 8); // Local work size dimensions
int globalWorkSizeT = OpenCLUtil.roundUp(localWorkSize, maxTIndex); // rounded up to the nearest multiple of localWorkSize
//int globalWorkSizeBeta = OpenCLUtil.roundUp(localWorkSize, maxBetaIndex); // rounded up to the nearest multiple of localWorkSize
// load sources, create and build program
CLProgram program = null;
try {
program = context.createProgram(this.getClass().getResourceAsStream("FanBeamProjector.cl"))
.build();
} catch (IOException e) {
// TODO Auto-generated catch block
e.printStackTrace();
System.exit(-1);
}
// create image from input grid
CLImageFormat format = new CLImageFormat(ChannelOrder.INTENSITY, ChannelType.FLOAT);
CLBuffer<FloatBuffer> imageBuffer = context.createFloatBuffer(imageSize, Mem.READ_ONLY);
// for (int i = 0; i < grid.getSize()[0]; ++i) {
// imageBuffer.getBuffer().put(grid.getSubGrid(i).getBuffer());
// }
for (int i=0;i<grid.getSize()[1];++i){
for (int j=0;j<grid.getSize()[0];++j)
imageBuffer.getBuffer().put(grid.getAtIndex(j, i));
}
imageBuffer.getBuffer().rewind();
CLImage2d<FloatBuffer> imageGrid = context.createImage2d(
imageBuffer.getBuffer(), grid.getSize()[0], grid.getSize()[1],
format);
imageBuffer.release();
// create memory for sinogram
CLBuffer<FloatBuffer> sinogram = context.createFloatBuffer(maxTIndex , Mem.WRITE_ONLY);//different from 2D openCL
// copy params
CLKernel kernel = program.createCLKernel("projectRayDriven1DCL");
kernel.putArg(imageGrid).putArg(sinogram)
.putArg((float)maxT).putArg((float)deltaT)
.putArg((float)maxBeta).putArg((float)deltaBeta)
.putArg((float)focalLength).putArg(maxTIndex).putArg(maxBetaIndex).putArg(index); // TODO: Spacing :)
// createCommandQueue
CLCommandQueue queue = device.createCommandQueue();
queue
.putWriteImage(imageGrid, true)
.finish()
.put1DRangeKernel(kernel, 0, globalWorkSizeT,//***********************different from 2D opencl
localWorkSize).putBarrier()
.putReadBuffer(sinogram, true)
.finish();
// write sinogram back to grid2D
Grid1D sino = new Grid1D(maxTIndex);
sino.setSpacing(deltaT);
sinogram.getBuffer().rewind();
for (int i = 0; i < sino.getBuffer().length; ++i) {
//sino.getBuffer()[i] = sinogram.getBuffer().get();
sino.setAtIndex(i, sinogram.getBuffer().get());
}
// clean up
queue.release();
imageGrid.release();
sinogram.release();
kernel.release();
program.release();
context.release();
return sino;
}
public void fastProjectRayDrivenCL(OpenCLGrid1D sinoCL, OpenCLGrid2D gridCL,int index) {
CLContext context = OpenCLUtil.getStaticContext();
CLDevice device = context.getMaxFlopsDevice();
int localWorkSize = Math.min(device.getMaxWorkGroupSize(), 8); // Local work size dimensions
int globalWorkSizeT = OpenCLUtil.roundUp(localWorkSize, maxTIndex); // rounded up to the nearest multiple of localWorkSize
// load sources, create and build program
CLProgram program = null;
try {
program = context.createProgram(this.getClass().getResourceAsStream("FanBeamProjector.cl"))
.build();
} catch (IOException e) {
// TODO Auto-generated catch block
e.printStackTrace();
System.exit(-1);
}
// create image from input grid
CLImageFormat format = new CLImageFormat(ChannelOrder.INTENSITY, ChannelType.FLOAT);
CLImage2d<FloatBuffer> imageGrid = context.createImage2d(gridCL.getDelegate().getCLBuffer().getBuffer(), gridCL.getSize()[0], gridCL.getSize()[1],format);
// copy params
CLKernel kernel = program.createCLKernel("projectRayDriven1DCL");
kernel.putArg(imageGrid).putArg(sinoCL.getDelegate().getCLBuffer())
.putArg((float)maxT).putArg((float)deltaT)
.putArg((float)maxBeta).putArg((float)deltaBeta)
.putArg((float)focalLength).putArg(maxTIndex).putArg(maxBetaIndex).putArg(index);
// createCommandQueue
CLCommandQueue queue = device.createCommandQueue();
queue
.putCopyBufferToImage(gridCL.getDelegate().getCLBuffer(),imageGrid).finish()
.put1DRangeKernel(kernel, 0, globalWorkSizeT,localWorkSize).putBarrier()
.finish();
kernel.rewind();
sinoCL.getDelegate().notifyDeviceChange();
// clean up
queue.release();
imageGrid.release();
kernel.release();
program.release();
}
}
/*
* Copyright (C) 2010-2014 Andreas Maier
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/