package edu.stanford.rsl.tutorial.fan;
import java.io.IOException;
import java.nio.FloatBuffer;
import java.util.ArrayList;
import java.util.Vector;
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.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.SimpleOperators;
import edu.stanford.rsl.conrad.numerics.SimpleVector;
import edu.stanford.rsl.conrad.opencl.OpenCLUtil;
/**This backprojector removed the distance weighting from the backprojector in the file tutorial.fan,
* since the ART algorithm only needs a backprojection without any filtering
**/
public class FanBeamBackprojector2D {
final double samplingRate = 3.d;
private float focalLength,
deltaT, deltaBeta,
maxT, maxBeta;
private int maxTIndex, maxBetaIndex;
private int imgSizeX, imgSizeY;
public FanBeamBackprojector2D(double focalLength, double deltaT, double deltaBeta ,int imageSizeX, int imageSizeY) {
this.focalLength = (float) focalLength;
this.deltaBeta = (float) deltaBeta;
this.deltaT = (float) deltaT;
this.imgSizeX = imageSizeX;
this.imgSizeY = imageSizeY;
}
public void initSinogramParams(Grid2D sino) {
this.maxTIndex = sino.getSize()[0];
this.deltaT = (float) sino.getSpacing()[0];
this.maxT = maxTIndex * deltaT;
this.maxBetaIndex = sino.getSize()[1];
this.deltaBeta = (float) sino.getSpacing()[1];
this.maxBeta = (maxBetaIndex -1) * deltaBeta;
}
public void initSinogramParams(OpenCLGrid2D sinoCL) {
this.maxTIndex = sinoCL.getSize()[0];
this.deltaT = (float) sinoCL.getSpacing()[0];
this.maxT = maxTIndex * deltaT;
this.maxBetaIndex = sinoCL.getSize()[1];
this.deltaBeta = (float) sinoCL.getSpacing()[1];
this.maxBeta = (maxBetaIndex -1) * deltaBeta;
}
/**
* Ray driven implementation of the Backprojector. This methods still contains a bug.
* ParkerWeighting does not work. Use Pixel-Driven instead.
* @param sino the sinogram for backprojection
* @return the backprojected image
*/
public Grid2D backprojectRayDriven(Grid2D sino) {
this.initSinogramParams(sino);
Grid2D grid = new Grid2D(this.imgSizeX, this.imgSizeY);
// set up image bounding box in WC
Translation trans = new Translation(-grid.getSize()[0]/2.0, -grid.getSize()[1]/2.0, -1);
// set up the inverse transform
Transform inverse = trans.inverse();
Box b = new Box(grid.getSize()[0], grid.getSize()[1], 2);
b.applyTransform(trans);
// iterate over the projection angles
for(int e=0; e<maxBetaIndex; ++e){
// compute beta [rad] and angular functions.
float beta = (float) (deltaBeta * e);
float cosBeta = (float) Math.cos(beta);
float sinBeta = (float) Math.sin(beta);
// We use PointND for points in 3D space and SimpleVector for directions.
// compute source location and the beginning of the detector
PointND a = new PointND(focalLength * cosBeta, focalLength * sinBeta, 0.d);
PointND p0 = new PointND(-maxT / 2.f * sinBeta, maxT / 2.f * cosBeta, 0.d);
// compute the normalized vector along the detector
SimpleVector dirDetector = p0.getAbstractVector().multipliedBy(-1);
dirDetector.normalizeL2();
// iterate over the detector bins
for (int i = 0; i < maxTIndex; ++i) {
// compute bin position
float stepsDirection = (float) (0.5f * deltaT + i * deltaT);
PointND p = new PointND(p0);
p.getAbstractVector().add(dirDetector.multipliedBy(stepsDirection));
// set up line equation
StraightLine line = new StraightLine(a, p);
// compute intersections between bounding box and intersection line.
ArrayList<PointND> points = b.intersect(line);
// only if we have intersections
if (2 != points.size()){
if(points.size() == 0) {
line.getDirection().multiplyBy(-1.f);
points = b.intersect(line);
}
if(points.size() == 0)
continue;
}
PointND start = points.get(0);
PointND end = points.get(1);
// get the normalized increment
SimpleVector increment = new SimpleVector(end.getAbstractVector());
increment.subtract(start.getAbstractVector());
float distance = (float) increment.normL2();
increment.divideBy(distance * samplingRate);
float val = sino.getAtIndex(i, e);
start = inverse.transform(start);
// compute the integral along the line.
for (float t = 0.0f; t < distance * samplingRate; ++t) {
PointND current = new PointND(start);
current.getAbstractVector().add(increment.multipliedBy(t));
if (grid.getSize()[0] <= current.get(0) + 1
|| grid.getSize()[1] <= current.get(1) + 1
|| current.get(0) < 0 || current.get(1) < 0)
continue;
InterpolationOperators.addInterpolateLinear(grid, current.get(0), current.get(1), val);
}
}
}
float normalizationFactor = (float) ((float) samplingRate * maxBetaIndex / deltaT / Math.PI);
NumericPointwiseOperators.divideBy(grid, normalizationFactor);
return grid;
}
public Grid2D backprojectRayDriven1D( Grid1D sinodiff1D, int Betaindex) {
//back project the sino difference at the Betaindex-th angle
/*this.maxTIndex = sinodiff1D.getSize()[0];
this.maxBetaIndex=360;
this.deltaT = (float) sinodiff1D.getSpacing()[0];
this.maxT = maxTIndex * deltaT;*/
Grid2D grid = new Grid2D(this.imgSizeX, this.imgSizeY);
// set up image bounding box in WC
Translation trans = new Translation(-grid.getSize()[0]/2.0, -grid.getSize()[1]/2.0, -1);
// set up the inverse transform
Transform inverse = trans.inverse();
Box b = new Box(grid.getSize()[0], grid.getSize()[1], 2);
b.applyTransform(trans);
// compute beta [rad] and angular functions.
float beta = (float) (deltaBeta * Betaindex);
float cosBeta = (float) Math.cos(beta);
float sinBeta = (float) Math.sin(beta);
// We use PointND for points in 3D space and SimpleVector for directions.
// compute source location and the beginning of the detector
PointND a = new PointND(focalLength * cosBeta, focalLength * sinBeta, 0.d);
PointND p0 = new PointND(-maxT / 2.f * sinBeta, maxT / 2.f * cosBeta, 0.d);
// compute the normalized vector along the detector
SimpleVector dirDetector = p0.getAbstractVector().multipliedBy(-1);
dirDetector.normalizeL2();
// iterate over the detector bins
for (int i = 0; i < maxTIndex; ++i) {
// compute bin position
float stepsDirection = (float) (0.5f * deltaT + i * deltaT);
PointND p = new PointND(p0);
p.getAbstractVector().add(dirDetector.multipliedBy(stepsDirection));
// set up line equation
StraightLine line = new StraightLine(a, p);
// compute intersections between bounding box and intersection line.
ArrayList<PointND> points = b.intersect(line);
// only if we have intersections
if (2 != points.size()){
if(points.size() == 0) {
line.getDirection().multiplyBy(-1.f);
points = b.intersect(line);
}
if(points.size() == 0)
continue;
}
PointND start = points.get(0);
PointND end = points.get(1);
// get the normalized increment
SimpleVector increment = new SimpleVector(end.getAbstractVector());
increment.subtract(start.getAbstractVector());
float distance = (float) increment.normL2();
increment.divideBy(distance * samplingRate);
float val = sinodiff1D.getAtIndex(i);
start = inverse.transform(start);
// compute the integral along the line.
for (float t = 0.0f; t < distance * samplingRate; ++t) {
PointND current = new PointND(start);
current.getAbstractVector().add(increment.multipliedBy(t));
if (grid.getSize()[0] <= current.get(0) + 1
|| grid.getSize()[1] <= current.get(1) + 1
|| current.get(0) < 0 || current.get(1) < 0)
continue;
InterpolationOperators.addInterpolateLinear(grid, current.get(0), current.get(1), val);
}
}
float normalizationFactor = (float) ((float) samplingRate * maxBetaIndex / deltaT / Math.PI);
//float normalizationFactor = (float) ((float) samplingRate / deltaT / Math.PI);
NumericPointwiseOperators.divideBy(grid, normalizationFactor);
return grid;
}
/**
* The pixel driven solution for back-projection.
*
* @param sino
* the sinogram
* @return the image
*/
public Grid2D backprojectPixelDriven(Grid2D sino) {
this.initSinogramParams(sino);
Grid2D grid = new Grid2D(this.imgSizeX, this.imgSizeY);
for(int b=0; b<maxBetaIndex; ++b){
// compute beta [rad] and angular functions.
float beta = (float) (deltaBeta * b);
float cosBeta = (float) Math.cos(beta);
float sinBeta = (float) Math.sin(beta);
PointND a = new PointND(focalLength * cosBeta, focalLength * sinBeta, 0.d);
PointND p0 = new PointND(-maxT / 2.f * sinBeta, maxT / 2.f * cosBeta, 0.d);
SimpleVector dirDetector = p0.getAbstractVector().multipliedBy(-1);
StraightLine detectorLine = new StraightLine(p0, dirDetector);
Grid1D subSino = sino.getSubGrid(b);
for(int x=0; x<imgSizeX; ++x){
float wx = (x - imgSizeX/2.0f);
for (int y = 0; y < imgSizeY; ++y) {
float wy = (y - imgSizeY/2.0f);
// compute two points on the line through t and beta
// We use PointND for points in 3D space and SimpleVector for directions.
final PointND reconstructionPointWorld = new PointND(wx, wy, 0.d);
final StraightLine projectionLine = new StraightLine(a, reconstructionPointWorld);
final PointND detectorPixel = projectionLine.intersect(detectorLine);
float valtemp;
if (detectorPixel != null)
{
final SimpleVector p = SimpleOperators.subtract(
detectorPixel.getAbstractVector(), p0.getAbstractVector()
);
double len = p.normL2();
if((p.getElement(0)*dirDetector.getElement(0)+p.getElement(1)*dirDetector.getElement(1))<0)
len=-len;//*****************************************FIXME
double t = (len-0.5d)/deltaT;
if (subSino.getSize()[0] <= t + 1 || t < 0)
continue;
float val = InterpolationOperators.interpolateLinear(subSino, t);
/*
//DistanceWeighting
float radius = (float) reconstructionPointWorld.getAbstractVector().normL2();
float phi = (float) ((Math.PI/2) + Math.atan2(reconstructionPointWorld.get(1), reconstructionPointWorld.get(0)));
float dWeight = (float) ((focalLength +radius*Math.sin(beta - phi))/focalLength);
valtemp = (float) (val / (dWeight*dWeight));
*/
valtemp=val;
}
else
{
//final PointND detectorPixel2 = projectionLine.intersect(detectorLine);
//if (detectorPixel2 == null) {}
valtemp = 0.f;
//distWeights.setAtIndex(x, y, 100.0f);
}
grid.addAtIndex(x, y, valtemp);
//gridIntermediate.addAtIndex(x, y, valtemp);
}
}
} // end for
float normalizationFactor = (float) (maxBetaIndex / Math.PI);
NumericPointwiseOperators.divideBy(grid, normalizationFactor);
return grid;
}
/**
* The pixel driven solution for back-projection.
*
* @param sino
* the sinogram
* @return the image
*/
public Grid2D backprojectPixelDriven1D(Grid1D sinodiff1D, int Betaindex) {
Grid2D grid = new Grid2D(this.imgSizeX, this.imgSizeY);
float beta = (float) (deltaBeta * Betaindex);
float cosBeta = (float) Math.cos(beta);
float sinBeta = (float) Math.sin(beta);
PointND a = new PointND(focalLength * cosBeta, focalLength * sinBeta, 0.d);
PointND p0 = new PointND(-maxT / 2.f * sinBeta, maxT / 2.f * cosBeta, 0.d);
SimpleVector dirDetector = p0.getAbstractVector().multipliedBy(-1);
StraightLine detectorLine = new StraightLine(p0, dirDetector);
for(int x=0; x<imgSizeX; ++x){
float wx = (x - imgSizeX/2.0f);
for (int y = 0; y < imgSizeY; ++y) {
float wy = (y - imgSizeY/2.0f);
// compute two points on the line through t and beta
// We use PointND for points in 3D space and SimpleVector for directions.
final PointND reconstructionPointWorld = new PointND(wx, wy, 0.d);
final StraightLine projectionLine = new StraightLine(a, reconstructionPointWorld);
final PointND detectorPixel = projectionLine.intersect(detectorLine);
float val=0.f;
if (detectorPixel != null)
{
final SimpleVector p = SimpleOperators.subtract(
detectorPixel.getAbstractVector(), p0.getAbstractVector()
);
double len = p.normL2();
if((p.getElement(0)*dirDetector.getElement(0)+p.getElement(1)*dirDetector.getElement(1))<0)
len=-len;//*****************************************FIXME
double t = (len-0.5d)/deltaT;
if (sinodiff1D.getSize()[0] <= t + 1 || t < 0)
continue;
val = InterpolationOperators.interpolateLinear(sinodiff1D, t);
}
else
{
//final PointND detectorPixel2 = projectionLine.intersect(detectorLine);
//if (detectorPixel2 == null) {}
val = 0.f;
}
grid.addAtIndex(x, y, val);
}
}
//float normalizationFactor = (float) (maxBetaIndex / Math.PI);
float normalizationFactor =2.0f/deltaBeta;//FIXME
NumericPointwiseOperators.divideBy(grid, normalizationFactor);
return grid;
}
public Grid2D backprojectPixelDrivenCL(Grid2D sino) {
this.initSinogramParams(sino);
boolean debug = true;
// 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 sinoSize = maxBetaIndex*maxTIndex;
// Length of arrays to process
int localWorkSize = Math.min(device.getMaxWorkGroupSize(), 16); // Local work size dimensions
int globalWorkSizeX = OpenCLUtil.roundUp(localWorkSize, imgSizeX); // rounded up to the nearest multiple of localWorkSize
int globalWorkSizeY = OpenCLUtil.roundUp(localWorkSize, imgSizeY); // rounded up to the nearest multiple of localWorkSize
// load sources, create and build program
CLProgram program = null;
try {
program = context.createProgram(this.getClass().getResourceAsStream("FanBeamBackProjectorPixel.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> sinoBuffer = context.createFloatBuffer(sinoSize, Mem.READ_ONLY);
for (int i=0;i<sinoSize;++i){
sinoBuffer.getBuffer().put(((Grid2D)sino).getBuffer()[i]);
}
sinoBuffer.getBuffer().rewind();
/// CP
CLImage2d<FloatBuffer> sinoGrid = context.createImage2d(
sinoBuffer.getBuffer(), sino.getSize()[0], sino.getSize()[1],
format,Mem.READ_ONLY);
sinoBuffer.release();
// create memory for output image
CLBuffer<FloatBuffer> imgBuffer = context.createFloatBuffer(imgSizeX*imgSizeY, Mem.WRITE_ONLY);
// copy params
CLKernel kernel = program.createCLKernel("backprojectPixelDriven2DCL");
kernel.putArg(sinoGrid).putArg(imgBuffer)
.putArg(imgSizeX).putArg(imgSizeY)
.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(sinoGrid, true)
.finish()
.put2DRangeKernel(kernel, 0, 0, globalWorkSizeX, globalWorkSizeY,
localWorkSize, localWorkSize)
.finish()
.putReadBuffer(imgBuffer, true)
.finish();
// write grid back to grid2D
Grid2D img = new Grid2D(this.imgSizeX, this.imgSizeY);
//img.setSpacing(pxSzXMM, pxSzYMM);
imgBuffer.getBuffer().rewind();
for (int i = 0; i < imgSizeX*imgSizeY; ++i) {
((Grid2D)img).getBuffer()[i] = imgBuffer.getBuffer().get();
}
queue.release();
imgBuffer.release();
sinoGrid.release();
kernel.release();
program.release();
context.release();
return img;
}
public Grid2D backprojectPixelDriven1DCL(Grid1D projection,int index){
/*this.maxTIndex = projection.getNumberOfElements();
this.deltaT = (float) projection.getSpacing()[0];
this.maxT = maxTIndex * deltaT;*/
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 sinoSize =projection.getNumberOfElements();
// Length of arrays to process
int localWorkSize = Math.min(device.getMaxWorkGroupSize(), 16); // Local work size dimensions
int globalWorkSizeX = OpenCLUtil.roundUp(localWorkSize, imgSizeX); // rounded up to the nearest multiple of localWorkSize
int globalWorkSizeY = OpenCLUtil.roundUp(localWorkSize, imgSizeY); // rounded up to the nearest multiple of localWorkSize
// load sources, create and build program
CLProgram program = null;
try {
program = context.createProgram(this.getClass().getResourceAsStream("FanBeamBackProjectorPixel.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> sinoBuffer = context.createFloatBuffer(sinoSize, Mem.READ_ONLY);
for (int i=0;i<sinoSize;++i){
sinoBuffer.getBuffer().put(projection.getBuffer()[i]);//**********************
}
sinoBuffer.getBuffer().rewind();
/// CP
/*
CLImage2d<FloatBuffer> sinoGrid = context.createImage2d(
sinoBuffer.getBuffer(), sino.getSize()[0], sino.getSize()[1],
format);
*/
CLImage2d<FloatBuffer> sinoGrid = context.createImage2d(
sinoBuffer.getBuffer(), sinoSize, 1,
format,Mem.READ_ONLY);
sinoBuffer.release();
// create memory for output image
CLBuffer<FloatBuffer> imgBuffer = context.createFloatBuffer(imgSizeX*imgSizeY, Mem.WRITE_ONLY);
// copy params
CLKernel kernel = program.createCLKernel("backprojectPixelDriven1DCL");
kernel.putArg(sinoGrid).putArg(imgBuffer)
.putArg(imgSizeX).putArg(imgSizeY)
.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(sinoGrid, true)
.finish()
.put2DRangeKernel(kernel, 0, 0, globalWorkSizeX, globalWorkSizeY,
localWorkSize, localWorkSize)
.finish()
.putReadBuffer(imgBuffer, true)
.finish();
// write grid back to grid2D
Grid2D img = new Grid2D(this.imgSizeX, this.imgSizeY);
//img.setSpacing(pxSzXMM, pxSzYMM);
imgBuffer.getBuffer().rewind();
for (int i = 0; i < imgSizeX*imgSizeY; ++i) {
((Grid2D)img).getBuffer()[i] = imgBuffer.getBuffer().get();
}
queue.release();
imgBuffer.release();
sinoGrid.release();
kernel.release();
program.release();
context.release();
return img;
}
public void fastBackprojectPixelDrivenCL(OpenCLGrid2D sinoCL, OpenCLGrid2D gridCL) {
sinoCL.getDelegate().prepareForDeviceOperation();
gridCL.getDelegate().prepareForDeviceOperation();
this.initSinogramParams(sinoCL);
CLContext context = OpenCLUtil.getStaticContext();
CLDevice device = context.getMaxFlopsDevice();
int localWorkSize = Math.min(device.getMaxWorkGroupSize(), 16); // Local work size dimensions
int globalWorkSizeX = OpenCLUtil.roundUp(localWorkSize, imgSizeX); // rounded up to the nearest multiple of localWorkSize
int globalWorkSizeY = OpenCLUtil.roundUp(localWorkSize, imgSizeY); // rounded up to the nearest multiple of localWorkSize
// load sources, create and build program
CLProgram program = null;
try {
program = context.createProgram(this.getClass().getResourceAsStream("FanBeamBackProjectorPixel.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> sinoGrid = context.createImage2d(sinoCL.getDelegate().getCLBuffer().getBuffer(), sinoCL.getSize()[0], sinoCL.getSize()[1],format,Mem.READ_ONLY);
// copy params
CLKernel kernel = program.createCLKernel("backprojectPixelDriven2DCL");
kernel.putArg(sinoGrid).putArg(gridCL.getDelegate().getCLBuffer())
.putArg(imgSizeX).putArg(imgSizeY)
.putArg((float)maxT).putArg((float)deltaT)
.putArg((float)maxBeta).putArg((float)deltaBeta)
.putArg((float)focalLength).putArg(maxTIndex).putArg(maxBetaIndex);
System.out.println(globalWorkSizeY +" "+globalWorkSizeX+ " "+localWorkSize);
// createCommandQueue
CLCommandQueue queue = device.createCommandQueue();
queue
.putCopyBufferToImage(sinoCL.getDelegate().getCLBuffer(),sinoGrid).finish()
.put2DRangeKernel(kernel, 0, 0, globalWorkSizeX, globalWorkSizeY,localWorkSize, localWorkSize).putBarrier()
.finish();
gridCL.getDelegate().notifyDeviceChange();
sinoGrid.release();
kernel.rewind();
queue.release();
kernel.release();
program.release();
}
public void fastBackprojectRayDrivenCL(OpenCLGrid2D sinoCL, OpenCLGrid2D gridCL) {
sinoCL.getDelegate().prepareForDeviceOperation();
gridCL.getDelegate().prepareForDeviceOperation();
this.initSinogramParams(sinoCL);
CLContext context = OpenCLUtil.getStaticContext();
CLDevice device = context.getMaxFlopsDevice();
// Length of arrays to process
int localWorkSize = Math.min(device.getMaxWorkGroupSize(), 16); // 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("FanBeamBackProjectorRay.cl"))
.build();
} catch (IOException e) {
// TODO Auto-generated catch block
e.printStackTrace();
System.exit(-1);
}
// copy params
CLKernel kernel = program.createCLKernel("backprojectRayDriven2DCL");
kernel.putArg(sinoCL.getDelegate().getCLBuffer()).putArg(gridCL.getDelegate().getCLBuffer())
.putArg(imgSizeX).putArg(imgSizeY)
.putArg((float)maxT).putArg((float)deltaT)
.putArg((float)maxBeta).putArg((float)deltaBeta)
.putArg((float)focalLength).putArg(maxTIndex).putArg(maxBetaIndex);
// createCommandQueue
CLCommandQueue queue = device.createCommandQueue();
queue
.put2DRangeKernel(kernel, 0, 0, globalWorkSizeBeta, globalWorkSizeT,localWorkSize, localWorkSize).putBarrier()
.finish();
gridCL.getDelegate().notifyDeviceChange();
queue.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).
*/