From b6cc56e921e301cb7fb9d7e0dadd601217eadbf8 Mon Sep 17 00:00:00 2001 From: Jennifer Date: Tue, 23 Nov 2021 14:44:56 +0100 Subject: [PATCH 1/3] Added option for deformed back-projection based on MLS deformation. Added source point computation based on SimpleMatrix instead of Jama.Matrix. --- .../conrad/opencl/OpenCLBackProjector.java | 194 +++++- .../conrad/opencl/OpenCLForwardProjector.java | 8 +- .../rsl/conrad/opencl/backprojectCL.cl | 589 ++++++++++++++++++ 3 files changed, 768 insertions(+), 23 deletions(-) diff --git a/src/edu/stanford/rsl/conrad/opencl/OpenCLBackProjector.java b/src/edu/stanford/rsl/conrad/opencl/OpenCLBackProjector.java index 488a84f3..1879b6cc 100644 --- a/src/edu/stanford/rsl/conrad/opencl/OpenCLBackProjector.java +++ b/src/edu/stanford/rsl/conrad/opencl/OpenCLBackProjector.java @@ -21,11 +21,17 @@ import com.jogamp.opencl.CLImageFormat.ChannelType; import com.jogamp.opencl.CLMemory.Mem; +import Jama.Matrix; +import Jama.SingularValueDecomposition; import edu.stanford.rsl.apps.gui.Citeable; import edu.stanford.rsl.conrad.data.numeric.Grid2D; import edu.stanford.rsl.conrad.data.numeric.Grid3D; +import edu.stanford.rsl.conrad.geometry.trajectories.Trajectory; import edu.stanford.rsl.conrad.io.ImagePlusDataSink; 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.numerics.SimpleVector.VectorNormType; import edu.stanford.rsl.conrad.reconstruction.VOIBasedReconstructionFilter; import edu.stanford.rsl.conrad.utils.CONRAD; import edu.stanford.rsl.conrad.utils.Configuration; @@ -77,6 +83,13 @@ public class OpenCLBackProjector extends VOIBasedReconstructionFilter implements */ protected float h_volume[]; + /** + * Deformed 3D coordinates + */ + protected SimpleVector[] p = null; + protected SimpleVector[][] qs = null; + private CLBuffer pBuff; + private CLBuffer qBuff; /** * The global variable of the module which stores the @@ -95,6 +108,8 @@ public class OpenCLBackProjector extends VOIBasedReconstructionFilter implements private boolean initialized = false; + + public OpenCLBackProjector () { super(); } @@ -187,8 +202,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) @@ -326,7 +344,7 @@ 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(); @@ -334,7 +352,7 @@ protected synchronized void initProjectionData(Grid2D projection){ 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); @@ -399,6 +417,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; @@ -407,28 +430,69 @@ 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 { +// System.out.println("pBuff = " + pBuff.getBuffer().get(0) + ", " + pBuff.getBuffer().get(1) + ", " + pBuff.getBuffer().get(2)); +// System.out.println("qBuff = " + qBuff.getBuffer().get(0) + ", " + qBuff.getBuffer().get(1) + ", " + qBuff.getBuffer().get(2)); + 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((float) p[0].getElement(0)) + .putArg((float) p[0].getElement(1)) + .putArg((float) p[0].getElement(2)) + .putArg((float) p[1].getElement(0)) + .putArg((float) p[1].getElement(1)) + .putArg((float) p[1].getElement(2)) + .putArg((float) p[2].getElement(0)) + .putArg((float) p[2].getElement(1)) + .putArg((float) p[2].getElement(2)) + .putArg((float) qs[projectionNumber][0].getElement(0)) + .putArg((float) qs[projectionNumber][0].getElement(1)) + .putArg((float) qs[projectionNumber][0].getElement(2)) + .putArg((float) qs[projectionNumber][1].getElement(0)) + .putArg((float) qs[projectionNumber][1].getElement(1)) + .putArg((float) qs[projectionNumber][1].getElement(2)) + .putArg((float) qs[projectionNumber][2].getElement(0)) + .putArg((float) qs[projectionNumber][2].getElement(1)) + .putArg((float) qs[projectionNumber][2].getElement(2)); + } + + - 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){ @@ -447,6 +511,85 @@ protected synchronized void projectSingleProjection(int projectionNumber, int di .finish(); } + private void initDeformedCoordinates(int projectionNumber) { + + SimpleVector[] q = this.qs[projectionNumber]; + + float[] p_floats = { + (float)p[0].getElement(0), (float)p[0].getElement(1), (float)p[0].getElement(2), + (float)p[1].getElement(0), (float)p[1].getElement(1), (float)p[1].getElement(2), + (float)p[2].getElement(0), (float)p[2].getElement(1), (float)p[2].getElement(2)}; + + float[] q_floats = { + (float)q[0].getElement(0), (float)q[0].getElement(1), (float)q[0].getElement(2), + (float)q[1].getElement(0), (float)q[1].getElement(1), (float)q[1].getElement(2), + (float)q[2].getElement(0), (float)q[2].getElement(1), (float)q[2].getElement(2)}; + + CLBuffer pBuffer = context.createFloatBuffer(p.length*p[0].getLen(), Mem.READ_ONLY); + pBuffer.getBuffer().put(p_floats); + pBuffer.getBuffer().rewind(); + this.pBuff = pBuffer; + + CLBuffer qBuffer = context.createFloatBuffer(q.length*q[0].getLen(), Mem.READ_ONLY); + qBuffer.getBuffer().put(q_floats); + qBuffer.getBuffer().rewind(); + this.qBuff = qBuffer; + } + + public static SimpleVector mls_rigid_deformation3D(Grid3D volume, SimpleVector coords, SimpleVector[] pIn, SimpleVector[] qIn, double alpha) { + // Rigid deformation + // Params: + // image - ndarray: original image + // p - ndarray: an array with size [n, 2], original control points + // q - ndarray: an array with size [n, 2], final control points + // alpha - float: parameter used by weights + // density - float: density of the grids + // Return: + // A deformed image. + + + SimpleVector[] p = pIn.clone(); + SimpleVector[] q = qIn.clone(); + SimpleVector v = coords.clone(); + + // wi + double[] w = new double[q.length]; + double wSum = 0.0; + for (int i = 0; i < w.length; i++) { + w[i] = 1.0/Math.pow(SimpleOperators.subtract(p[i], v).norm(VectorNormType.VEC_NORM_L2), 2.0*alpha); + wSum += w[i]; + } + + // q* and p* + SimpleVector qStar = new SimpleVector(q[0].getLen()); + SimpleVector pStar = new SimpleVector(p[0].getLen()); + for (int i = 0; i < w.length; i++) { + qStar.add(q[i].multipliedBy(w[i])); + pStar.add(p[i].multipliedBy(w[i])); + } + qStar.divideBy(wSum); + pStar.divideBy(wSum); + + // PQ' + SimpleMatrix PQtrans = new SimpleMatrix(pIn[0].getLen(), qIn[0].getLen()); + for (int i = 0; i < pIn.length; i++) { + SimpleMatrix pi_min_pStar = new SimpleMatrix(new double[][] {SimpleOperators.subtract(p[i], pStar).copyAsDoubleArray()}); + SimpleMatrix qi_min_qStar = new SimpleMatrix(new double[][] {SimpleOperators.subtract(q[i], qStar).copyAsDoubleArray()}); + PQtrans.add(SimpleOperators.multiplyMatrices(pi_min_pStar.transposed(), qi_min_qStar).multipliedBy(w[i])); + } + SingularValueDecomposition svd = new SingularValueDecomposition(new Matrix(PQtrans.copyAsDoubleArray())); + SimpleMatrix U = new SimpleMatrix(svd.getU().getArrayCopy()); + SimpleMatrix V = new SimpleMatrix(svd.getV().getArrayCopy()); + + // transform input coords + SimpleVector v_min_pStar = SimpleOperators.subtract(v, pStar); + SimpleMatrix VUtrans = SimpleOperators.multiplyMatrices(V,U.transposed()); + SimpleVector v_transformed = SimpleOperators.add(SimpleOperators.multiply(VUtrans, v_min_pStar), qStar); + + return v_transformed; + + } + public void loadInputQueue(Grid3D input) throws IOException { ImageGridBuffer buf = new ImageGridBuffer(); buf.set(input); @@ -566,6 +709,7 @@ public void OpenCLRun() { private synchronized void workOnProjectionData(){ if (projectionsAvailable.size() > 0){ + System.out.println(projectionsAvailable.size()); Integer current = projectionsAvailable.get(0); projectionsAvailable.remove(0); projectSingleProjection(current.intValue(), @@ -644,6 +788,12 @@ public String getToolName(){ return "OpenCL Backprojector"; } + public void setDeformParameters(SimpleVector[] p, SimpleVector[][] qs) { + this.p = p; + this.qs = qs; + } + + } /* * Copyright (C) 2010-2014 Martin Berger diff --git a/src/edu/stanford/rsl/conrad/opencl/OpenCLForwardProjector.java b/src/edu/stanford/rsl/conrad/opencl/OpenCLForwardProjector.java index cc7010c5..8731757a 100644 --- a/src/edu/stanford/rsl/conrad/opencl/OpenCLForwardProjector.java +++ b/src/edu/stanford/rsl/conrad/opencl/OpenCLForwardProjector.java @@ -227,7 +227,7 @@ public void computeCanonicalProjectionMatrix(CLBuffer 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))); @@ -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 diff --git a/src/edu/stanford/rsl/conrad/opencl/backprojectCL.cl b/src/edu/stanford/rsl/conrad/opencl/backprojectCL.cl index 9da47d27..f0dbdc26 100644 --- a/src/edu/stanford/rsl/conrad/opencl/backprojectCL.cl +++ b/src/edu/stanford/rsl/conrad/opencl/backprojectCL.cl @@ -96,6 +96,595 @@ __kernel void backprojectKernel( return; } +/* -------------------------------------------------------------------------- + * + * + * Voxel-based BP algorithm implementation in OpenCL kernel programming + * WITH DEFORMATION! + * + * + * -------------------------------------------------------------------------- */ +float4 mls_deformation(float v0, float v1, float v2, float p00, float p01, float p02, float p10, float p11, float p12, float p20, float p21, float p22, float q00, float q01, float q02, float q10, float q11, float q12, float q20, float q21, float q22) +{ + // wi + float w0 = 1/((p00 - v0) * (p00 - v0) + (p01 - v1) * (p01 - v1) + (p02 - v2) * (p02 - v2)); + float w1 = 1/((p10 - v0) * (p10 - v0) + (p11 - v1) * (p11 - v1) + (p12 - v2) * (p12 - v2)); + float w2 = 1/((p20 - v0) * (p20 - v0) + (p21 - v1) * (p21 - v1) + (p22 - v2) * (p22 - v2)); + float wSum = w0 + w1 + w2; + + // q* and p* + float qStar0 = (q00 * w0 + q10 * w1 + q20 * w2) / wSum; + float qStar1 = (q01 * w0 + q11 * w1 + q21 * w2) / wSum; + float qStar2 = (q02 * w0 + q12 * w1 + q22 * w2) / wSum; + + float pStar0 = (p00 * w0 + p10 * w1 + p20 * w2) / wSum; + float pStar1 = (p01 * w0 + p11 * w1 + p21 * w2) / wSum; + float pStar2 = (p02 * w0 + p12 * w1 + p22 * w2) / wSum; + + + // PQ' + float p0_min_pStar0 = p00 - pStar0; + float p0_min_pStar1 = p01 - pStar1; + float p0_min_pStar2 = p02 - pStar2; + + float p1_min_pStar0 = p10 - pStar0; + float p1_min_pStar1 = p11 - pStar1; + float p1_min_pStar2 = p12 - pStar2; + + float p2_min_pStar0 = p20 - pStar0; + float p2_min_pStar1 = p21 - pStar1; + float p2_min_pStar2 = p22 - pStar2; + + float q0_min_qStar0 = q00 - qStar0; + float q0_min_qStar1 = q01 - qStar1; + float q0_min_qStar2 = q02 - qStar2; + + float q1_min_qStar0 = q10 - qStar0; + float q1_min_qStar1 = q11 - qStar1; + float q1_min_qStar2 = q12 - qStar2; + + float q2_min_qStar0 = q20 - qStar0; + float q2_min_qStar1 = q21 - qStar1; + float q2_min_qStar2 = q22 - qStar2; + + float PQtrans0 = w0*p0_min_pStar0*q0_min_qStar0 + w1*p1_min_pStar0*q1_min_qStar0 + w2*p2_min_pStar0*q2_min_qStar0; + float PQtrans3 = w0*p0_min_pStar0*q0_min_qStar1 + w1*p1_min_pStar0*q1_min_qStar1 + w2*p2_min_pStar0*q2_min_qStar1; + float PQtrans6 = w0*p0_min_pStar0*q0_min_qStar2 + w1*p1_min_pStar0*q1_min_qStar2 + w2*p2_min_pStar0*q2_min_qStar2; + + float PQtrans1 = w0*p0_min_pStar1*q0_min_qStar0 + w1*p1_min_pStar1*q1_min_qStar0 + w2*p2_min_pStar1*q2_min_qStar0; + float PQtrans4 = w0*p0_min_pStar1*q0_min_qStar1 + w1*p1_min_pStar1*q1_min_qStar1 + w2*p2_min_pStar1*q2_min_qStar1; + float PQtrans7 = w0*p0_min_pStar1*q0_min_qStar2 + w1*p1_min_pStar1*q1_min_qStar2 + w2*p2_min_pStar1*q2_min_qStar2; + + float PQtrans2 = w0*p0_min_pStar2*q0_min_qStar0 + w1*p1_min_pStar2*q1_min_qStar0 + w2*p2_min_pStar2*q2_min_qStar0; + float PQtrans5 = w0*p0_min_pStar2*q0_min_qStar1 + w1*p1_min_pStar2*q1_min_qStar1 + w2*p2_min_pStar2*q2_min_qStar1; + float PQtrans8 = w0*p0_min_pStar2*q0_min_qStar2 + w1*p1_min_pStar2*q1_min_qStar2 + w2*p2_min_pStar2*q2_min_qStar2; + + // svd + int n = 3; + + double tol = 0.000000000000001f; + int loopmax = 100 * n; + int loopcount = 0; + + float U0 = 1; + float U1 = 0; + float U2 = 0; + float U3 = 0; + float U4 = 1; + float U5 = 0; + float U6 = 0; + float U7 = 0; + float U8 = 1; + + float V0 = 1; + float V1 = 0; + float V2 = 0; + float V3 = 0; + float V4 = 1; + float V5 = 0; + float V6 = 0; + float V7 = 0; + float V8 = 1; + + float S0 = PQtrans0; + float S1 = PQtrans3; + float S2 = PQtrans6; + float S3 = PQtrans1; + float S4 = PQtrans4; + float S5 = PQtrans7; + float S6 = PQtrans2; + float S7 = PQtrans5; + float S8 = PQtrans8; + + float err = 10000; + + while (err > tol && loopcount < loopmax) { + + // qr ////////////////////////// + int m = 3; + + float R0 = S0; + float R1 = S3; + float R2 = S6; + float R3 = S1; + float R4 = S4; + float R5 = S7; + float R6 = S2; + float R7 = S5; + float R8 = S8; + + float Q0 = 1; + float Q1 = 0; + float Q2 = 0; + float Q3 = 0; + float Q4 = 1; + float Q5 = 0; + float Q6 = 0; + float Q7 = 0; + float Q8 = 1; + + // k = 0; + float x0 = R0; + float x1 = R1; + float x2 = R2; + + float norm_x = sqrt(x0*x0 + x1*x1 + x2*x2); + float g = -sign(x0) * norm_x; + x0 = x0 - g; + + // Orthogonal transformation matrix that eliminates one element + // below the diagonal of the matrix it is post-multiplying: + norm_x = sqrt(x0*x0 + x1*x1 + x2*x2); + + if (norm_x != 0) { + + x0 = x0/norm_x; + x1 = x1/norm_x; + x2 = x2/norm_x; + + float u1 = 2 * (R0*x0 + R1*x1 + R2*x2); + float u2 = 2 * (R3*x0 + R4*x1 + R5*x2); + float u3 = 2 * (R6*x0 + R7*x1 + R8*x2); + + R0 = R0 - x0*u1; + R3 = R3 - x0*u2; + R6 = R6 - x0*u3; + R1 = R1 - x1*u1; + R4 = R4 - x1*u2; + R7 = R7 - x1*u3; + R2 = R2 - x2*u1; + R5 = R5 - x2*u2; + R8 = R8 - x2*u3; + + float Q0tmp = Q0 - 2*(Q0*x0*x0 + Q3*x1*x0 + Q6*x2*x0); + float Q3tmp = Q3 - 2*(Q0*x0*x1 + Q3*x1*x1 + Q6*x2*x1); + float Q6tmp = Q6 - 2*(Q0*x0*x2 + Q3*x1*x2 + Q6*x2*x2); + float Q1tmp = Q1 - 2*(Q1*x0*x0 + Q4*x1*x0 + Q7*x2*x0); + float Q4tmp = Q4 - 2*(Q1*x0*x1 + Q4*x1*x1 + Q7*x2*x1); + float Q7tmp = Q7 - 2*(Q1*x0*x2 + Q4*x1*x2 + Q7*x2*x2); + float Q2tmp = Q2 - 2*(Q2*x0*x0 + Q5*x1*x0 + Q8*x2*x0); + float Q5tmp = Q5 - 2*(Q2*x0*x1 + Q5*x1*x1 + Q8*x2*x1); + float Q8tmp = Q8 - 2*(Q2*x0*x2 + Q5*x1*x2 + Q8*x2*x2); + + Q0 = Q0tmp; + Q1 = Q1tmp; + Q2 = Q2tmp; + Q3 = Q3tmp; + Q4 = Q4tmp; + Q5 = Q5tmp; + Q6 = Q6tmp; + Q7 = Q7tmp; + Q8 = Q8tmp; + } + + // k = 1; + x0 = 0; + x1 = R4; + x2 = R5; + + norm_x = sqrt(x0*x0 + x1*x1 + x2*x2); + g = -sign(x1) * norm_x; + x1 = x1 - g; + + // Orthogonal transformation matrix that eliminates one element + // below the diagonal of the matrix it is post-multiplying: + norm_x = sqrt(x0*x0 + x1*x1 + x2*x2); + + if (norm_x != 0) { + + x0 = x0/norm_x; + x1 = x1/norm_x; + x2 = x2/norm_x; + + float u1 = 2 * (R0*x0 + R1*x1 + R2*x2); + float u2 = 2 * (R3*x0 + R4*x1 + R5*x2); + float u3 = 2 * (R6*x0 + R7*x1 + R8*x2); + + R0 = R0 - x0*u1; + R3 = R3 - x0*u2; + R6 = R6 - x0*u3; + R1 = R1 - x1*u1; + R4 = R4 - x1*u2; + R7 = R7 - x1*u3; + R2 = R2 - x2*u1; + R5 = R5 - x2*u2; + R8 = R8 - x2*u3; + + float Q0tmp = Q0 - 2*(Q0*x0*x0 + Q3*x1*x0 + Q6*x2*x0); + float Q3tmp = Q3 - 2*(Q0*x0*x1 + Q3*x1*x1 + Q6*x2*x1); + float Q6tmp = Q6 - 2*(Q0*x0*x2 + Q3*x1*x2 + Q6*x2*x2); + float Q1tmp = Q1 - 2*(Q1*x0*x0 + Q4*x1*x0 + Q7*x2*x0); + float Q4tmp = Q4 - 2*(Q1*x0*x1 + Q4*x1*x1 + Q7*x2*x1); + float Q7tmp = Q7 - 2*(Q1*x0*x2 + Q4*x1*x2 + Q7*x2*x2); + float Q2tmp = Q2 - 2*(Q2*x0*x0 + Q5*x1*x0 + Q8*x2*x0); + float Q5tmp = Q5 - 2*(Q2*x0*x1 + Q5*x1*x1 + Q8*x2*x1); + float Q8tmp = Q8 - 2*(Q2*x0*x2 + Q5*x1*x2 + Q8*x2*x2); + + Q0 = Q0tmp; + Q1 = Q1tmp; + Q2 = Q2tmp; + Q3 = Q3tmp; + Q4 = Q4tmp; + Q5 = Q5tmp; + Q6 = Q6tmp; + Q7 = Q7tmp; + Q8 = Q8tmp; + } + + ////////////////////////////////////// + + S0 = R0; + S1 = R1; + S2 = R2; + S3 = R3; + S4 = R4; + S5 = R5; + S6 = R6; + S7 = R7; + S8 = R8; + + float U0tmp = U0 * Q0 + U3 * Q1 + U6 * Q2; + float U1tmp = U1 * Q0 + U4 * Q1 + U7 * Q2; + float U2tmp = U2 * Q0 + U5 * Q1 + U8 * Q2; + float U3tmp = U0 * Q3 + U3 * Q4 + U6 * Q5; + float U4tmp = U1 * Q3 + U4 * Q4 + U7 * Q5; + float U5tmp = U2 * Q3 + U5 * Q4 + U8 * Q5; + float U6tmp = U0 * Q6 + U3 * Q7 + U6 * Q8; + float U7tmp = U1 * Q6 + U4 * Q7 + U7 * Q8; + float U8tmp = U2 * Q6 + U5 * Q7 + U8 * Q8; + U0 = U0tmp; + U1 = U1tmp; + U2 = U2tmp; + U3 = U3tmp; + U4 = U4tmp; + U5 = U5tmp; + U6 = U6tmp; + U7 = U7tmp; + U8 = U8tmp; + + // qr ////////////////////////// + R0 = S0; + R1 = S3; + R2 = S6; + R3 = S1; + R4 = S4; + R5 = S7; + R6 = S2; + R7 = S5; + R8 = S8; + + Q0 = 1; + Q1 = 0; + Q2 = 0; + Q3 = 0; + Q4 = 1; + Q5 = 0; + Q6 = 0; + Q7 = 0; + Q8 = 1; + + // k = 0; + x0 = R0; + x1 = R1; + x2 = R2; + + norm_x = sqrt(x0*x0 + x1*x1 + x2*x2); + g = -sign(x0) * norm_x; + x0 = x0 - g; + + // Orthogonal transformation matrix that eliminates one element + // below the diagonal of the matrix it is post-multiplying: + norm_x = sqrt(x0*x0 + x1*x1 + x2*x2); + + if (norm_x != 0) { + + x0 = x0/norm_x; + x1 = x1/norm_x; + x2 = x2/norm_x; + + float u1 = 2 * (R0*x0 + R1*x1 + R2*x2); + float u2 = 2 * (R3*x0 + R4*x1 + R5*x2); + float u3 = 2 * (R6*x0 + R7*x1 + R8*x2); + + R0 = R0 - x0*u1; + R3 = R3 - x0*u2; + R6 = R6 - x0*u3; + R1 = R1 - x1*u1; + R4 = R4 - x1*u2; + R7 = R7 - x1*u3; + R2 = R2 - x2*u1; + R5 = R5 - x2*u2; + R8 = R8 - x2*u3; + + float Q0tmp = Q0 - 2*(Q0*x0*x0 + Q3*x1*x0 + Q6*x2*x0); + float Q3tmp = Q3 - 2*(Q0*x0*x1 + Q3*x1*x1 + Q6*x2*x1); + float Q6tmp = Q6 - 2*(Q0*x0*x2 + Q3*x1*x2 + Q6*x2*x2); + float Q1tmp = Q1 - 2*(Q1*x0*x0 + Q4*x1*x0 + Q7*x2*x0); + float Q4tmp = Q4 - 2*(Q1*x0*x1 + Q4*x1*x1 + Q7*x2*x1); + float Q7tmp = Q7 - 2*(Q1*x0*x2 + Q4*x1*x2 + Q7*x2*x2); + float Q2tmp = Q2 - 2*(Q2*x0*x0 + Q5*x1*x0 + Q8*x2*x0); + float Q5tmp = Q5 - 2*(Q2*x0*x1 + Q5*x1*x1 + Q8*x2*x1); + float Q8tmp = Q8 - 2*(Q2*x0*x2 + Q5*x1*x2 + Q8*x2*x2); + + Q0 = Q0tmp; + Q1 = Q1tmp; + Q2 = Q2tmp; + Q3 = Q3tmp; + Q4 = Q4tmp; + Q5 = Q5tmp; + Q6 = Q6tmp; + Q7 = Q7tmp; + Q8 = Q8tmp; + } + + // k = 1; + x0 = 0; + x1 = R4; + x2 = R5; + + norm_x = sqrt(x0*x0 + x1*x1 + x2*x2); + g = -sign(x1) * norm_x; + x1 = x1 - g; + + // Orthogonal transformation matrix that eliminates one element + // below the diagonal of the matrix it is post-multiplying: + norm_x = sqrt(x0*x0 + x1*x1 + x2*x2); + + if (norm_x != 0) { + + x0 = x0/norm_x; + x1 = x1/norm_x; + x2 = x2/norm_x; + + float u1 = 2 * (R0*x0 + R1*x1 + R2*x2); + float u2 = 2 * (R3*x0 + R4*x1 + R5*x2); + float u3 = 2 * (R6*x0 + R7*x1 + R8*x2); + + R0 = R0 - x0*u1; + R3 = R3 - x0*u2; + R6 = R6 - x0*u3; + R1 = R1 - x1*u1; + R4 = R4 - x1*u2; + R7 = R7 - x1*u3; + R2 = R2 - x2*u1; + R5 = R5 - x2*u2; + R8 = R8 - x2*u3; + + float Q0tmp = Q0 - 2*(Q0*x0*x0 + Q3*x1*x0 + Q6*x2*x0); + float Q3tmp = Q3 - 2*(Q0*x0*x1 + Q3*x1*x1 + Q6*x2*x1); + float Q6tmp = Q6 - 2*(Q0*x0*x2 + Q3*x1*x2 + Q6*x2*x2); + float Q1tmp = Q1 - 2*(Q1*x0*x0 + Q4*x1*x0 + Q7*x2*x0); + float Q4tmp = Q4 - 2*(Q1*x0*x1 + Q4*x1*x1 + Q7*x2*x1); + float Q7tmp = Q7 - 2*(Q1*x0*x2 + Q4*x1*x2 + Q7*x2*x2); + float Q2tmp = Q2 - 2*(Q2*x0*x0 + Q5*x1*x0 + Q8*x2*x0); + float Q5tmp = Q5 - 2*(Q2*x0*x1 + Q5*x1*x1 + Q8*x2*x1); + float Q8tmp = Q8 - 2*(Q2*x0*x2 + Q5*x1*x2 + Q8*x2*x2); + + Q0 = Q0tmp; + Q1 = Q1tmp; + Q2 = Q2tmp; + Q3 = Q3tmp; + Q4 = Q4tmp; + Q5 = Q5tmp; + Q6 = Q6tmp; + Q7 = Q7tmp; + Q8 = Q8tmp; + } + + ////////////////////////////////////// + + S0 = R0; + S1 = R1; + S2 = R2; + S3 = R3; + S4 = R4; + S5 = R5; + S6 = R6; + S7 = R7; + S8 = R8; + + float V0tmp = V0 * Q0 + V3 * Q1 + V6 * Q2; + float V1tmp = V1 * Q0 + V4 * Q1 + V7 * Q2; + float V2tmp = V2 * Q0 + V5 * Q1 + V8 * Q2; + float V3tmp = V0 * Q3 + V3 * Q4 + V6 * Q5; + float V4tmp = V1 * Q3 + V4 * Q4 + V7 * Q5; + float V5tmp = V2 * Q3 + V5 * Q4 + V8 * Q5; + float V6tmp = V0 * Q6 + V3 * Q7 + V6 * Q8; + float V7tmp = V1 * Q6 + V4 * Q7 + V7 * Q8; + float V8tmp = V2 * Q6 + V5 * Q7 + V8 * Q8; + V0 = V0tmp; + V1 = V1tmp; + V2 = V2tmp; + V3 = V3tmp; + V4 = V4tmp; + V5 = V5tmp; + V6 = V6tmp; + V7 = V7tmp; + V8 = V8tmp; + + // exit when we get "close" + float E = sqrt(S3*S3 + S6*S6 + S7*S7); + float F = sqrt(S0*S0 + S4*S4 + S8*S8); + if (F == 0) { + F = 1; + } + err = E/F; + loopcount = loopcount + 1; + } + + // fix the signs in S + float SS0 = S0; + float SS1 = S4; + float SS2 = S8; + + if (SS0 < -0.000000001) { + U0 = -U0; + U1 = -U1; + U2 = -U2; + } + + if (SS1 < -0.000000001) { + U3 = -U3; + U4 = -U4; + U5 = -U5; + } + + if (SS2 < -0.000000001) { + U6 = -U6; + U7 = -U7; + U8 = -U8; + } + + // transform input coords + float v_min_pStar0 = v0 - pStar0; + float v_min_pStar1 = v1 - pStar1; + float v_min_pStar2 = v2 - pStar2; + + float VUtrans0 = V0*U0 + V3*U3 + V6*U6; + float VUtrans1 = V1*U0 + V4*U3 + V7*U6; + float VUtrans2 = V2*U0 + V5*U3 + V8*U6; + float VUtrans3 = V0*U1 + V3*U4 + V6*U7; + float VUtrans4 = V1*U1 + V4*U4 + V7*U7; + float VUtrans5 = V2*U1 + V5*U4 + V8*U7; + float VUtrans6 = V0*U2 + V3*U5 + V6*U8; + float VUtrans7 = V1*U2 + V4*U5 + V7*U8; + float VUtrans8 = V2*U2 + V5*U5 + V8*U8; + + if (VUtrans0*VUtrans4*VUtrans8 + VUtrans3*VUtrans7*VUtrans2 + VUtrans6*VUtrans1*VUtrans5 - VUtrans2*VUtrans4*VUtrans6 - VUtrans5*VUtrans7*VUtrans0 - VUtrans8*VUtrans1*VUtrans3 <= 0) + { + V6 = -V6; + V7 = -V7; + V8 = -V8; + VUtrans0 = V0*U0 + V3*U3 + V6*U6; + VUtrans1 = V1*U0 + V4*U3 + V7*U6; + VUtrans2 = V2*U0 + V5*U3 + V8*U6; + VUtrans3 = V0*U1 + V3*U4 + V6*U7; + VUtrans4 = V1*U1 + V4*U4 + V7*U7; + VUtrans5 = V2*U1 + V5*U4 + V8*U7; + VUtrans6 = V0*U2 + V3*U5 + V6*U8; + VUtrans7 = V1*U2 + V4*U5 + V7*U8; + VUtrans8 = V2*U2 + V5*U5 + V8*U8; + } + + float v_transformed0 = VUtrans0*v_min_pStar0 + VUtrans3*v_min_pStar1 + VUtrans6*v_min_pStar2 + qStar0; + float v_transformed1 = VUtrans1*v_min_pStar0 + VUtrans4*v_min_pStar1 + VUtrans7*v_min_pStar2 + qStar1; + float v_transformed2 = VUtrans2*v_min_pStar0 + VUtrans5*v_min_pStar1 + VUtrans8*v_min_pStar2 + qStar2; + + float4 v_transformed = (float4)(v_transformed0, v_transformed1, v_transformed2, 0.0f); + return v_transformed; + +} + + +__kernel void backprojectKernelDeformed( + __global TvoxelValue* pVolume, + int sizeX, + int sizeY, + int recoSizeZ, + int offset, + float voxelSpacingX, + float voxelSpacingY, + float voxelSpacingZ, + float offsetX, + float offsetY, + float offsetZ, + __read_only image2d_t gTex2D, + __constant Tcoord_dev* gProjMatrix, + float projMultiplier, + float p00, + float p01, + float p02, + float p10, + float p11, + float p12, + float p20, + float p21, + float p22, + float q00, + float q01, + float q02, + float q10, + float q11, + float q12, + float q20, + float q21, + float q22) +{ + int gidx = get_group_id(0); + int gidy = get_group_id(1); + int lidx = get_local_id(0); + int lidy = get_local_id(1); + + int locSizex = get_local_size(0); + int locSizey = get_local_size(1); + + unsigned int x = gidx*locSizex+lidx; + unsigned int y = gidy*locSizey+lidy; + + if (x >= sizeX || y >= sizeY) + return; + + unsigned int zStride = sizeX*sizeY; + unsigned int yStride = sizeX; + + for (unsigned int z = 0; z < recoSizeZ; z++){ + + float4 coords_new = mls_deformation((x * voxelSpacingX) - offsetX, (y * voxelSpacingY) - offsetY, (z * voxelSpacingZ) - offsetZ, p00, p01, p02, p10, p11, p12, p20, p21, p22, q00, q01, q02, q10, q11, q12, q20, q21, q22); + float xcoord = coords_new.x; + float ycoord = coords_new.y; + float zcoord = coords_new.z; + //float xcoord = (x * voxelSpacingX) - offsetX + (q10-p10); + //float ycoord = (y * voxelSpacingY) - offsetY + (q11-p11); + //float zcoord = (z * voxelSpacingZ) - offsetZ + (q12-p12); + + // precompute everything except z from the matrix multiplication + float precomputeR = gProjMatrix[3] + ycoord * gProjMatrix[1] + xcoord * gProjMatrix[0]; + float precomputeS = gProjMatrix[7] + ycoord * gProjMatrix[5] + xcoord * gProjMatrix[4]; + float precomputeT = gProjMatrix[11] + ycoord * gProjMatrix[9] + xcoord * gProjMatrix[8]; + + // compute homogeneous coordinates + float r = zcoord * gProjMatrix[2] + precomputeR; + float s = zcoord * gProjMatrix[6] + precomputeS; + float t = zcoord * gProjMatrix[10] + precomputeT; + + // compute projection coordinates + float denom = 1.0f / t; + float fu = r * denom; + float fv = s * denom; + + float proj_val = read_imagef(gTex2D, sampler, (float2)(fu + 0.5f + offset, fv + 0.5f)).x; + + // compute volume index for x,y,z + unsigned long idx = z*zStride + y*yStride + x; + + // Distance weighting by 1/t^2 + pVolume[idx] += proj_val * denom * denom * projMultiplier; + + } + + return; +} + /* * Copyright (C) 2010-2014 Martin Berger From 493f19b27cca2f410db6fc690dc49c24c937d559 Mon Sep 17 00:00:00 2001 From: Jennifer Date: Tue, 30 Nov 2021 10:46:40 +0100 Subject: [PATCH 2/3] Updated MLS deformation to be able to process arbitrary number of points for deformation. Added comments. --- .../conrad/opencl/OpenCLBackProjector.java | 145 +++++---------- .../rsl/conrad/opencl/backprojectCL.cl | 166 +++++++++--------- 2 files changed, 128 insertions(+), 183 deletions(-) diff --git a/src/edu/stanford/rsl/conrad/opencl/OpenCLBackProjector.java b/src/edu/stanford/rsl/conrad/opencl/OpenCLBackProjector.java index 1879b6cc..00b529af 100644 --- a/src/edu/stanford/rsl/conrad/opencl/OpenCLBackProjector.java +++ b/src/edu/stanford/rsl/conrad/opencl/OpenCLBackProjector.java @@ -21,23 +21,28 @@ import com.jogamp.opencl.CLImageFormat.ChannelType; import com.jogamp.opencl.CLMemory.Mem; -import Jama.Matrix; -import Jama.SingularValueDecomposition; import edu.stanford.rsl.apps.gui.Citeable; import edu.stanford.rsl.conrad.data.numeric.Grid2D; import edu.stanford.rsl.conrad.data.numeric.Grid3D; -import edu.stanford.rsl.conrad.geometry.trajectories.Trajectory; import edu.stanford.rsl.conrad.io.ImagePlusDataSink; 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.numerics.SimpleVector.VectorNormType; 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{ /** @@ -228,7 +233,7 @@ protected void init(){ commandQueue. putWriteBuffer(volumePointer, true). - finish(); + finish(); initialized = true; } @@ -449,8 +454,6 @@ protected synchronized void projectSingleProjection(int projectionNumber, int di .putArg(projectionMatrix) .putArg(projectionMultiplier); } else { -// System.out.println("pBuff = " + pBuff.getBuffer().get(0) + ", " + pBuff.getBuffer().get(1) + ", " + pBuff.getBuffer().get(2)); -// System.out.println("qBuff = " + qBuff.getBuffer().get(0) + ", " + qBuff.getBuffer().get(1) + ", " + qBuff.getBuffer().get(2)); kernelFunction.rewind(); kernelFunction .putArg(volumePointer) @@ -467,28 +470,12 @@ protected synchronized void projectSingleProjection(int projectionNumber, int di .putArg(projectionTex) .putArg(projectionMatrix) .putArg(projectionMultiplier) - .putArg((float) p[0].getElement(0)) - .putArg((float) p[0].getElement(1)) - .putArg((float) p[0].getElement(2)) - .putArg((float) p[1].getElement(0)) - .putArg((float) p[1].getElement(1)) - .putArg((float) p[1].getElement(2)) - .putArg((float) p[2].getElement(0)) - .putArg((float) p[2].getElement(1)) - .putArg((float) p[2].getElement(2)) - .putArg((float) qs[projectionNumber][0].getElement(0)) - .putArg((float) qs[projectionNumber][0].getElement(1)) - .putArg((float) qs[projectionNumber][0].getElement(2)) - .putArg((float) qs[projectionNumber][1].getElement(0)) - .putArg((float) qs[projectionNumber][1].getElement(1)) - .putArg((float) qs[projectionNumber][1].getElement(2)) - .putArg((float) qs[projectionNumber][2].getElement(0)) - .putArg((float) qs[projectionNumber][2].getElement(1)) - .putArg((float) qs[projectionNumber][2].getElement(2)); + .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])); @@ -503,27 +490,37 @@ 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 = { - (float)p[0].getElement(0), (float)p[0].getElement(1), (float)p[0].getElement(2), - (float)p[1].getElement(0), (float)p[1].getElement(1), (float)p[1].getElement(2), - (float)p[2].getElement(0), (float)p[2].getElement(1), (float)p[2].getElement(2)}; - - float[] q_floats = { - (float)q[0].getElement(0), (float)q[0].getElement(1), (float)q[0].getElement(2), - (float)q[1].getElement(0), (float)q[1].getElement(1), (float)q[1].getElement(2), - (float)q[2].getElement(0), (float)q[2].getElement(1), (float)q[2].getElement(2)}; + 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 pBuffer = context.createFloatBuffer(p.length*p[0].getLen(), Mem.READ_ONLY); pBuffer.getBuffer().put(p_floats); @@ -536,60 +533,6 @@ private void initDeformedCoordinates(int projectionNumber) { this.qBuff = qBuffer; } - public static SimpleVector mls_rigid_deformation3D(Grid3D volume, SimpleVector coords, SimpleVector[] pIn, SimpleVector[] qIn, double alpha) { - // Rigid deformation - // Params: - // image - ndarray: original image - // p - ndarray: an array with size [n, 2], original control points - // q - ndarray: an array with size [n, 2], final control points - // alpha - float: parameter used by weights - // density - float: density of the grids - // Return: - // A deformed image. - - - SimpleVector[] p = pIn.clone(); - SimpleVector[] q = qIn.clone(); - SimpleVector v = coords.clone(); - - // wi - double[] w = new double[q.length]; - double wSum = 0.0; - for (int i = 0; i < w.length; i++) { - w[i] = 1.0/Math.pow(SimpleOperators.subtract(p[i], v).norm(VectorNormType.VEC_NORM_L2), 2.0*alpha); - wSum += w[i]; - } - - // q* and p* - SimpleVector qStar = new SimpleVector(q[0].getLen()); - SimpleVector pStar = new SimpleVector(p[0].getLen()); - for (int i = 0; i < w.length; i++) { - qStar.add(q[i].multipliedBy(w[i])); - pStar.add(p[i].multipliedBy(w[i])); - } - qStar.divideBy(wSum); - pStar.divideBy(wSum); - - // PQ' - SimpleMatrix PQtrans = new SimpleMatrix(pIn[0].getLen(), qIn[0].getLen()); - for (int i = 0; i < pIn.length; i++) { - SimpleMatrix pi_min_pStar = new SimpleMatrix(new double[][] {SimpleOperators.subtract(p[i], pStar).copyAsDoubleArray()}); - SimpleMatrix qi_min_qStar = new SimpleMatrix(new double[][] {SimpleOperators.subtract(q[i], qStar).copyAsDoubleArray()}); - PQtrans.add(SimpleOperators.multiplyMatrices(pi_min_pStar.transposed(), qi_min_qStar).multipliedBy(w[i])); - } - SingularValueDecomposition svd = new SingularValueDecomposition(new Matrix(PQtrans.copyAsDoubleArray())); - SimpleMatrix U = new SimpleMatrix(svd.getU().getArrayCopy()); - SimpleMatrix V = new SimpleMatrix(svd.getV().getArrayCopy()); - - // transform input coords - SimpleVector v_min_pStar = SimpleOperators.subtract(v, pStar); - SimpleMatrix VUtrans = SimpleOperators.multiplyMatrices(V,U.transposed()); - SimpleVector v_transformed = SimpleOperators.add(SimpleOperators.multiply(VUtrans, v_min_pStar), qStar); - - return v_transformed; - - } - public void loadInputQueue(Grid3D input) throws IOException { ImageGridBuffer buf = new ImageGridBuffer(); buf.set(input); @@ -788,6 +731,14 @@ 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; @@ -796,6 +747,6 @@ public void setDeformParameters(SimpleVector[] p, SimpleVector[][] 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). */ \ No newline at end of file diff --git a/src/edu/stanford/rsl/conrad/opencl/backprojectCL.cl b/src/edu/stanford/rsl/conrad/opencl/backprojectCL.cl index f0dbdc26..c9798f0d 100644 --- a/src/edu/stanford/rsl/conrad/opencl/backprojectCL.cl +++ b/src/edu/stanford/rsl/conrad/opencl/backprojectCL.cl @@ -104,60 +104,70 @@ __kernel void backprojectKernel( * * * -------------------------------------------------------------------------- */ -float4 mls_deformation(float v0, float v1, float v2, float p00, float p01, float p02, float p10, float p11, float p12, float p20, float p21, float p22, float q00, float q01, float q02, float q10, float q11, float q12, float q20, float q21, float q22) +float4 mls_deformation(float x, float y, float z, global float* p, global float* q, int noPoints) { - // wi - float w0 = 1/((p00 - v0) * (p00 - v0) + (p01 - v1) * (p01 - v1) + (p02 - v2) * (p02 - v2)); - float w1 = 1/((p10 - v0) * (p10 - v0) + (p11 - v1) * (p11 - v1) + (p12 - v2) * (p12 - v2)); - float w2 = 1/((p20 - v0) * (p20 - v0) + (p21 - v1) * (p21 - v1) + (p22 - v2) * (p22 - v2)); - float wSum = w0 + w1 + w2; - + // q* and p* - float qStar0 = (q00 * w0 + q10 * w1 + q20 * w2) / wSum; - float qStar1 = (q01 * w0 + q11 * w1 + q21 * w2) / wSum; - float qStar2 = (q02 * w0 + q12 * w1 + q22 * w2) / wSum; - - float pStar0 = (p00 * w0 + p10 * w1 + p20 * w2) / wSum; - float pStar1 = (p01 * w0 + p11 * w1 + p21 * w2) / wSum; - float pStar2 = (p02 * w0 + p12 * w1 + p22 * w2) / wSum; - - + float qStar0 = 0; + float qStar1 = 0; + float qStar2 = 0; + float pStar0 = 0; + float pStar1 = 0; + float pStar2 = 0; + float wSum = 0; + + for (unsigned int idx = 0; idx < noPoints; idx++) { + float w = 1/((p[idx*3] - x) * (p[idx*3] - x) + (p[idx*3+1] - y) * (p[idx*3+1] - y) + (p[idx*3+2] - z) * (p[idx*3+2] - z)); + wSum += w; + qStar0 += q[idx*3] * w; + qStar1 += q[idx*3+1] * w; + qStar2 += q[idx*3+2] * w; + pStar0 += p[idx*3] * w; + pStar1 += p[idx*3+1] * w; + pStar2 += p[idx*3+2] * w; + } + qStar0 /= wSum; + qStar1 /= wSum; + qStar2 /= wSum; + pStar0 /= wSum; + pStar1 /= wSum; + pStar2 /= wSum; + + // PQ' - float p0_min_pStar0 = p00 - pStar0; - float p0_min_pStar1 = p01 - pStar1; - float p0_min_pStar2 = p02 - pStar2; - - float p1_min_pStar0 = p10 - pStar0; - float p1_min_pStar1 = p11 - pStar1; - float p1_min_pStar2 = p12 - pStar2; - - float p2_min_pStar0 = p20 - pStar0; - float p2_min_pStar1 = p21 - pStar1; - float p2_min_pStar2 = p22 - pStar2; - - float q0_min_qStar0 = q00 - qStar0; - float q0_min_qStar1 = q01 - qStar1; - float q0_min_qStar2 = q02 - qStar2; - - float q1_min_qStar0 = q10 - qStar0; - float q1_min_qStar1 = q11 - qStar1; - float q1_min_qStar2 = q12 - qStar2; - - float q2_min_qStar0 = q20 - qStar0; - float q2_min_qStar1 = q21 - qStar1; - float q2_min_qStar2 = q22 - qStar2; - - float PQtrans0 = w0*p0_min_pStar0*q0_min_qStar0 + w1*p1_min_pStar0*q1_min_qStar0 + w2*p2_min_pStar0*q2_min_qStar0; - float PQtrans3 = w0*p0_min_pStar0*q0_min_qStar1 + w1*p1_min_pStar0*q1_min_qStar1 + w2*p2_min_pStar0*q2_min_qStar1; - float PQtrans6 = w0*p0_min_pStar0*q0_min_qStar2 + w1*p1_min_pStar0*q1_min_qStar2 + w2*p2_min_pStar0*q2_min_qStar2; - - float PQtrans1 = w0*p0_min_pStar1*q0_min_qStar0 + w1*p1_min_pStar1*q1_min_qStar0 + w2*p2_min_pStar1*q2_min_qStar0; - float PQtrans4 = w0*p0_min_pStar1*q0_min_qStar1 + w1*p1_min_pStar1*q1_min_qStar1 + w2*p2_min_pStar1*q2_min_qStar1; - float PQtrans7 = w0*p0_min_pStar1*q0_min_qStar2 + w1*p1_min_pStar1*q1_min_qStar2 + w2*p2_min_pStar1*q2_min_qStar2; - - float PQtrans2 = w0*p0_min_pStar2*q0_min_qStar0 + w1*p1_min_pStar2*q1_min_qStar0 + w2*p2_min_pStar2*q2_min_qStar0; - float PQtrans5 = w0*p0_min_pStar2*q0_min_qStar1 + w1*p1_min_pStar2*q1_min_qStar1 + w2*p2_min_pStar2*q2_min_qStar1; - float PQtrans8 = w0*p0_min_pStar2*q0_min_qStar2 + w1*p1_min_pStar2*q1_min_qStar2 + w2*p2_min_pStar2*q2_min_qStar2; + float PQtrans0 = 0; + float PQtrans1 = 0; + float PQtrans2 = 0; + float PQtrans3 = 0; + float PQtrans4 = 0; + float PQtrans5 = 0; + float PQtrans6 = 0; + float PQtrans7 = 0; + float PQtrans8 = 0; + for (unsigned int idx = 0; idx < noPoints; idx++) { + float w = 1/((p[idx*3] - x) * (p[idx*3] - x) + (p[idx*3+1] - y) * (p[idx*3+1] - y) + (p[idx*3+2] - z) * (p[idx*3+2] - z)); + + float p_min_pStar0 = p[idx*3] - pStar0; + float p_min_pStar1 = p[idx*3+1] - pStar1; + float p_min_pStar2 = p[idx*3+2] - pStar2; + + float q_min_qStar0 = q[idx*3] - qStar0; + float q_min_qStar1 = q[idx*3+1] - qStar1; + float q_min_qStar2 = q[idx*3+2] - qStar2; + + PQtrans0 += w*p_min_pStar0*q_min_qStar0; + PQtrans3 += w*p_min_pStar0*q_min_qStar1; + PQtrans6 += w*p_min_pStar0*q_min_qStar2; + + PQtrans1 += w*p_min_pStar1*q_min_qStar0; + PQtrans4 += w*p_min_pStar1*q_min_qStar1; + PQtrans7 += w*p_min_pStar1*q_min_qStar2; + + PQtrans2 += w*p_min_pStar2*q_min_qStar0; + PQtrans5 += w*p_min_pStar2*q_min_qStar1; + PQtrans8 += w*p_min_pStar2*q_min_qStar2; + + } // svd int n = 3; @@ -556,9 +566,9 @@ float4 mls_deformation(float v0, float v1, float v2, float p00, float p01, float } // transform input coords - float v_min_pStar0 = v0 - pStar0; - float v_min_pStar1 = v1 - pStar1; - float v_min_pStar2 = v2 - pStar2; + float x_min_pStar = x - pStar0; + float y_min_pStar = y - pStar1; + float z_min_pStar = z - pStar2; float VUtrans0 = V0*U0 + V3*U3 + V6*U6; float VUtrans1 = V1*U0 + V4*U3 + V7*U6; @@ -586,11 +596,11 @@ float4 mls_deformation(float v0, float v1, float v2, float p00, float p01, float VUtrans8 = V2*U2 + V5*U5 + V8*U8; } - float v_transformed0 = VUtrans0*v_min_pStar0 + VUtrans3*v_min_pStar1 + VUtrans6*v_min_pStar2 + qStar0; - float v_transformed1 = VUtrans1*v_min_pStar0 + VUtrans4*v_min_pStar1 + VUtrans7*v_min_pStar2 + qStar1; - float v_transformed2 = VUtrans2*v_min_pStar0 + VUtrans5*v_min_pStar1 + VUtrans8*v_min_pStar2 + qStar2; - - float4 v_transformed = (float4)(v_transformed0, v_transformed1, v_transformed2, 0.0f); + float x_transformed = VUtrans0*x_min_pStar + VUtrans3*y_min_pStar + VUtrans6*z_min_pStar + qStar0; + float y_transformed = VUtrans1*x_min_pStar + VUtrans4*y_min_pStar + VUtrans7*z_min_pStar + qStar1; + float z_transformed = VUtrans2*x_min_pStar + VUtrans5*y_min_pStar + VUtrans8*z_min_pStar + qStar2; + + float4 v_transformed = (float4)(x_transformed, y_transformed, z_transformed, 0.0f); return v_transformed; } @@ -611,25 +621,11 @@ __kernel void backprojectKernelDeformed( __read_only image2d_t gTex2D, __constant Tcoord_dev* gProjMatrix, float projMultiplier, - float p00, - float p01, - float p02, - float p10, - float p11, - float p12, - float p20, - float p21, - float p22, - float q00, - float q01, - float q02, - float q10, - float q11, - float q12, - float q20, - float q21, - float q22) + global float* p, + global float* q, + int noPoints) { + int gidx = get_group_id(0); int gidy = get_group_id(1); int lidx = get_local_id(0); @@ -649,14 +645,11 @@ __kernel void backprojectKernelDeformed( for (unsigned int z = 0; z < recoSizeZ; z++){ - float4 coords_new = mls_deformation((x * voxelSpacingX) - offsetX, (y * voxelSpacingY) - offsetY, (z * voxelSpacingZ) - offsetZ, p00, p01, p02, p10, p11, p12, p20, p21, p22, q00, q01, q02, q10, q11, q12, q20, q21, q22); + float4 coords_new = mls_deformation((x * voxelSpacingX) - offsetX, (y * voxelSpacingY) - offsetY, (z * voxelSpacingZ) - offsetZ, p, q, noPoints); float xcoord = coords_new.x; float ycoord = coords_new.y; float zcoord = coords_new.z; - //float xcoord = (x * voxelSpacingX) - offsetX + (q10-p10); - //float ycoord = (y * voxelSpacingY) - offsetY + (q11-p11); - //float zcoord = (z * voxelSpacingZ) - offsetZ + (q12-p12); - + // precompute everything except z from the matrix multiplication float precomputeR = gProjMatrix[3] + ycoord * gProjMatrix[1] + xcoord * gProjMatrix[0]; float precomputeS = gProjMatrix[7] + ycoord * gProjMatrix[5] + xcoord * gProjMatrix[4]; @@ -666,14 +659,14 @@ __kernel void backprojectKernelDeformed( float r = zcoord * gProjMatrix[2] + precomputeR; float s = zcoord * gProjMatrix[6] + precomputeS; float t = zcoord * gProjMatrix[10] + precomputeT; - + // compute projection coordinates float denom = 1.0f / t; float fu = r * denom; float fv = s * denom; - + float proj_val = read_imagef(gTex2D, sampler, (float2)(fu + 0.5f + offset, fv + 0.5f)).x; - + // compute volume index for x,y,z unsigned long idx = z*zStride + y*yStride + x; @@ -686,7 +679,8 @@ __kernel void backprojectKernelDeformed( } + /* - * 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). */ \ No newline at end of file From 8d56d532a5af9a745293197ca000411418c57867 Mon Sep 17 00:00:00 2001 From: Jennifer Date: Tue, 30 Nov 2021 12:00:48 +0100 Subject: [PATCH 3/3] removed console print out --- src/edu/stanford/rsl/conrad/opencl/OpenCLBackProjector.java | 1 - 1 file changed, 1 deletion(-) diff --git a/src/edu/stanford/rsl/conrad/opencl/OpenCLBackProjector.java b/src/edu/stanford/rsl/conrad/opencl/OpenCLBackProjector.java index 00b529af..ff01bbbb 100644 --- a/src/edu/stanford/rsl/conrad/opencl/OpenCLBackProjector.java +++ b/src/edu/stanford/rsl/conrad/opencl/OpenCLBackProjector.java @@ -652,7 +652,6 @@ public void OpenCLRun() { private synchronized void workOnProjectionData(){ if (projectionsAvailable.size() > 0){ - System.out.println(projectionsAvailable.size()); Integer current = projectionsAvailable.get(0); projectionsAvailable.remove(0); projectSingleProjection(current.intValue(),