Skip to content

Mls deformation #149

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jenniferMaier wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Mls deformation #149

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
b6cc56e
Added option for deformed back-projection based on MLS deformation.
jenniferMaier Nov 23, 2021
493f19b
Updated MLS deformation to be able to process arbitrary number of poi…
jenniferMaier Nov 30, 2021
8d56d53
removed console print out
jenniferMaier Nov 30, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
160 changes: 130 additions & 30 deletions src/edu/stanford/rsl/conrad/opencl/OpenCLBackProjector.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,12 +26,23 @@
import edu.stanford.rsl.conrad.data.numeric.Grid3D;
import edu.stanford.rsl.conrad.io.ImagePlusDataSink;
import edu.stanford.rsl.conrad.numerics.SimpleMatrix;
import edu.stanford.rsl.conrad.numerics.SimpleVector;
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;
import edu.stanford.rsl.conrad.utils.ImageUtil;

/**
* Backprojection on GPU using OpenCL.
*
* If Moving Least Squares deformation should be included in the reconstruction, the deformation parameters (points
* before and after) need to be set using the method setDeformParameters() before the init() method is called.
*
* @author Martin Berger, Jennifer Maier
*
*/

public class OpenCLBackProjector extends VOIBasedReconstructionFilter implements Runnable, Citeable{

/**
Expand Down Expand Up @@ -77,6 +88,13 @@ public class OpenCLBackProjector extends VOIBasedReconstructionFilter implements
*/
protected float h_volume[];

/**
* Deformed 3D coordinates
*/
protected SimpleVector[] p = null;
protected SimpleVector[][] qs = null;
private CLBuffer<FloatBuffer> pBuff;
private CLBuffer<FloatBuffer> qBuff;

/**
* The global variable of the module which stores the
Expand All @@ -95,6 +113,8 @@ public class OpenCLBackProjector extends VOIBasedReconstructionFilter implements

private boolean initialized = false;



public OpenCLBackProjector () {
super();
}
Expand Down Expand Up @@ -187,8 +207,11 @@ protected void init(){
}

// create the computing kernel
kernelFunction = program.createCLKernel("backprojectKernel");

if (this.p == null) {
kernelFunction = program.createCLKernel("backprojectKernel");
} else {
kernelFunction = program.createCLKernel("backprojectKernelDeformed");
}
// create the reconstruction volume;
// createFloatBuffer uses a byteBuffer internally --> h_volume.length cannot be > 2^31/4 = 2^31/2^2 = 2^29
// Thus, 2^29 would already cause a overflow (negative sign) of the integer in the byte buffer! Maximum length is (2^29-1)
Expand All @@ -210,7 +233,7 @@ protected void init(){

commandQueue.
putWriteBuffer(volumePointer, true).
finish();
finish();

initialized = true;
}
Expand Down Expand Up @@ -326,15 +349,15 @@ protected synchronized void initProjectionData(Grid2D projection){
// 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(projection.getBuffer());
projectionArray.getBuffer().rewind();

if(projectionTex != null && !projectionTex.isReleased()){
projectionTex.release();
}

// set the texture
CLImageFormat format = new CLImageFormat(ChannelOrder.INTENSITY, ChannelType.FLOAT);
projectionTex = context.createImage2d(projectionArray.getBuffer(), projection.getWidth(), projection.getHeight(), format, Mem.READ_ONLY);
Expand Down Expand Up @@ -399,6 +422,11 @@ protected synchronized void projectSingleProjection(int projectionNumber, int di
if (!largeVolumeMode) {
projections.remove(projectionNumber);
}

if (this.p != null) {
initDeformedCoordinates(projectionNumber);
}

// backproject for each slice
// OpenCL Grids are only two dimensional!
int reconDimensionZ = dimz;
Expand All @@ -407,28 +435,51 @@ protected synchronized void projectSingleProjection(int projectionNumber, int di
double voxelSpacingZ = getGeometry().getVoxelSpacingZ();

// write kernel parameters
kernelFunction.rewind();
kernelFunction
.putArg(volumePointer)
.putArg(getGeometry().getReconDimensionX())
.putArg(getGeometry().getReconDimensionY())
.putArg(reconDimensionZ)
.putArg((int) lineOffset)
.putArg((float) voxelSpacingX)
.putArg((float) voxelSpacingY)
.putArg((float) voxelSpacingZ)
.putArg((float) offsetX)
.putArg((float) offsetY)
.putArg((float) offsetZ)
.putArg(projectionTex)
.putArg(projectionMatrix)
.putArg(projectionMultiplier);
// deformed coordinates
if (this.pBuff == null) {
kernelFunction.rewind();
kernelFunction
.putArg(volumePointer)
.putArg(getGeometry().getReconDimensionX())
.putArg(getGeometry().getReconDimensionY())
.putArg(reconDimensionZ)
.putArg((int) lineOffset)
.putArg((float) voxelSpacingX)
.putArg((float) voxelSpacingY)
.putArg((float) voxelSpacingZ)
.putArg((float) offsetX)
.putArg((float) offsetY)
.putArg((float) offsetZ)
.putArg(projectionTex)
.putArg(projectionMatrix)
.putArg(projectionMultiplier);
} else {
kernelFunction.rewind();
kernelFunction
.putArg(volumePointer)
.putArg(getGeometry().getReconDimensionX())
.putArg(getGeometry().getReconDimensionY())
.putArg(reconDimensionZ)
.putArg((int) lineOffset)
.putArg((float) voxelSpacingX)
.putArg((float) voxelSpacingY)
.putArg((float) voxelSpacingZ)
.putArg((float) offsetX)
.putArg((float) offsetY)
.putArg((float) offsetZ)
.putArg(projectionTex)
.putArg(projectionMatrix)
.putArg(projectionMultiplier)
.putArg(pBuff)
.putArg(qBuff)
.putArg((int) p.length);
}



int[] realLocalSize = new int[2];
realLocalSize[0] = Math.min(device.getMaxWorkGroupSize(),bpBlockSize[0]);
realLocalSize[1] = Math.max(1, Math.min(device.getMaxWorkGroupSize()/realLocalSize[0], bpBlockSize[1]));

// rounded up to the nearest multiple of localWorkSize
int[] globalWorkSize = {getGeometry().getReconDimensionX(), getGeometry().getReconDimensionY()};
if ((globalWorkSize[0] % realLocalSize[0] ) != 0){
Expand All @@ -439,12 +490,47 @@ protected synchronized void projectSingleProjection(int projectionNumber, int di
}

// Call the OpenCL kernel, writing the results into the volume which is pointed at
commandQueue
.putWriteImage(projectionTex, true)
.put2DRangeKernel(kernelFunction, 0, 0, globalWorkSize[0], globalWorkSize[1], realLocalSize[0], realLocalSize[1])
//.finish()
//.putReadBuffer(dOut, true)
.finish();
if (this.pBuff == null) {
commandQueue
.putWriteImage(projectionTex, true)
.put2DRangeKernel(kernelFunction, 0, 0, globalWorkSize[0], globalWorkSize[1], realLocalSize[0], realLocalSize[1])
.finish();
} else {
commandQueue
.putWriteImage(projectionTex, true)
.putWriteBuffer(pBuff, true)
.putWriteBuffer(qBuff, true)
.put2DRangeKernel(kernelFunction, 0, 0, globalWorkSize[0], globalWorkSize[1], realLocalSize[0], realLocalSize[1])
.finish();
}


}

private void initDeformedCoordinates(int projectionNumber) {

SimpleVector[] q = this.qs[projectionNumber];

float[] p_floats = new float[3*p.length];
float[] q_floats = new float[3*q.length];
for (int i = 0; i < p.length; i++) {
p_floats[3*i] = (float)p[i].getElement(0);
p_floats[3*i+1] = (float)p[i].getElement(1);
p_floats[3*i+2] = (float)p[i].getElement(2);
q_floats[3*i] = (float)q[i].getElement(0);
q_floats[3*i+1] = (float)q[i].getElement(1);
q_floats[3*i+2] = (float)q[i].getElement(2);
}

CLBuffer<FloatBuffer> pBuffer = context.createFloatBuffer(p.length*p[0].getLen(), Mem.READ_ONLY);
pBuffer.getBuffer().put(p_floats);
pBuffer.getBuffer().rewind();
this.pBuff = pBuffer;

CLBuffer<FloatBuffer> qBuffer = context.createFloatBuffer(q.length*q[0].getLen(), Mem.READ_ONLY);
qBuffer.getBuffer().put(q_floats);
qBuffer.getBuffer().rewind();
this.qBuff = qBuffer;
}

public void loadInputQueue(Grid3D input) throws IOException {
Expand Down Expand Up @@ -644,8 +730,22 @@ public String getToolName(){
return "OpenCL Backprojector";
}

/**
* set parameters used for deformation
*
* @param p original positions of points used for deformation. Contains n 3D-vectors if there are n points used
* for deformation
* @param qs deformed positions of points used for deformation. If there are p projections used for backprojection
* and n points used for deformation, the array contains p x n 3D-vectors.
*/
public void setDeformParameters(SimpleVector[] p, SimpleVector[][] qs) {
this.p = p;
this.qs = qs;
}


}
/*
* Copyright (C) 2010-2014 Martin Berger
* Copyright (C) 2010-2021 Martin Berger, Jennifer Maier
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/
8 changes: 7 additions & 1 deletion src/edu/stanford/rsl/conrad/opencl/OpenCLForwardProjector.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -227,7 +227,7 @@ public void computeCanonicalProjectionMatrix(CLBuffer<FloatBuffer> detectorDirec
originShift = getOriginTransform();

// New srcPoint in the Canonical coord sys
SimpleVector srcPtW = proj.computeCameraCenter().negated();//computeSrcPt(projectionMatrix, invARmatrixMat);
SimpleVector srcPtW = proj.computeCameraCenter().negated();//computeSrcPt(proj.computeP(), invARmatrixMat); //
srcPoint.getBuffer().put((float) -(-0.5 * (volumeSize[0] -1.0) + originShift.getElement(0)*invVoxelScale.getElement(0,0) + invVoxelScale.getElement(0,0) * srcPtW.getElement(0)));
srcPoint.getBuffer().put((float) -(-0.5 * (volumeSize[1] -1.0) + originShift.getElement(1)*invVoxelScale.getElement(1,1) + invVoxelScale.getElement(1,1) * srcPtW.getElement(1)));
srcPoint.getBuffer().put((float) -(-0.5 * (volumeSize[2] -1.0) + originShift.getElement(2)*invVoxelScale.getElement(2,2) + invVoxelScale.getElement(2,2) * srcPtW.getElement(2)));
Expand Down Expand Up @@ -281,6 +281,12 @@ public static Jama.Matrix computeSrcPt(Jama.Matrix projectionMatrix, Jama.Matrix
return invertedProjMatrix.times(at);
}

public static SimpleVector computeSrcPt(SimpleMatrix projectionMatrix, SimpleMatrix invertedProjMatrix) {
SimpleVector at = projectionMatrix.getSubCol(0, 3, 3);//.getMatrix(0, 2, 3, 3);
//at = at.times(-1.0);
return SimpleOperators.multiply(invertedProjMatrix, at);
}

protected SimpleVector getOriginTransform(){
SimpleVector currOrigin = new SimpleVector(this.origin);
// compute centered origin as assumed by forward projector
Expand Down
Loading