From 3eb6f61e48f7ab6ddaa0e78d9140a4322dbb92fb Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Thu, 2 Oct 2014 13:50:59 +0200
Subject: Add CUDA SIRT::doSlabCorrections() function
This function optionally compensates for effectively infinitely large slab-like
objects of finite thickness 1.
---
cuda/2d/arith.cu | 9 +++++++++
cuda/2d/arith.h | 1 +
cuda/2d/sirt.cu | 45 +++++++++++++++++++++++++++++++++++++++++++++
cuda/2d/sirt.h | 4 ++++
4 files changed, 59 insertions(+)
(limited to 'cuda/2d')
diff --git a/cuda/2d/arith.cu b/cuda/2d/arith.cu
index d1b189c..0f84e52 100644
--- a/cuda/2d/arith.cu
+++ b/cuda/2d/arith.cu
@@ -68,6 +68,14 @@ struct opMul {
out *= in;
}
};
+struct opDiv {
+ __device__ void operator()(float& out, const float in) {
+ if (in > 0.000001f) // out is assumed to be positive
+ out /= in;
+ else
+ out = 0.0f;
+ }
+};
struct opMul2 {
__device__ void operator()(float& out, const float in1, const float in2) {
out *= in1 * in2;
@@ -682,6 +690,7 @@ INST_DtoD(opDividedBy)
INST_toD(opInvert)
INST_FtoD(opSet)
INST_FtoD(opMul)
+INST_DtoD(opDiv)
INST_DFtoD(opMulMask)
INST_FtoD(opAdd)
INST_FtoD(opClampMin)
diff --git a/cuda/2d/arith.h b/cuda/2d/arith.h
index f87db99..31ecc17 100644
--- a/cuda/2d/arith.h
+++ b/cuda/2d/arith.h
@@ -41,6 +41,7 @@ struct opAddMul;
struct opAdd;
struct opAdd2;
struct opMul;
+struct opDiv;
struct opMul2;
struct opDividedBy;
struct opInvert;
diff --git a/cuda/2d/sirt.cu b/cuda/2d/sirt.cu
index c7fe219..8dcd553 100644
--- a/cuda/2d/sirt.cu
+++ b/cuda/2d/sirt.cu
@@ -142,6 +142,51 @@ bool SIRT::precomputeWeights()
return true;
}
+bool SIRT::doSlabCorrections()
+{
+ // This function compensates for effectively infinitely large slab-like
+ // objects of finite thickness 1.
+
+ // Each ray through the object has an intersection of length d/cos(alpha).
+ // The length of the ray actually intersecting the reconstruction volume is
+ // given by D_lineWeight. By dividing by 1/cos(alpha) and multiplying by the
+ // lineweights, we correct for this missing attenuation outside of the
+ // reconstruction volume, assuming the object is homogeneous.
+
+ // This effectively scales the output values by assuming the thickness d
+ // is 1 unit.
+
+
+ // This function in its current implementation only works if there are no masks.
+ // In this case, init() will also have already called precomputeWeights(),
+ // so we can use D_lineWeight.
+ if (useVolumeMask || useSinogramMask)
+ return false;
+
+ // multiply by line weights
+ processSino<opDiv>(D_sinoData, D_lineWeight, projPitch, dims);
+
+ SDimensions subdims = dims;
+ subdims.iProjAngles = 1;
+
+ // divide by 1/cos(angle)
+ // ...but limit the correction to -80/+80 degrees.
+ float bound = cosf(1.3963f);
+ float* t = (float*)D_sinoData;
+ for (int i = 0; i < dims.iProjAngles; ++i) {
+ float f = fabs(cosf(angles[i]));
+
+ if (f < bound)
+ f = bound;
+
+ processSino<opMul>(t, f, sinoPitch, subdims);
+ t += sinoPitch;
+ }
+
+ return true;
+}
+
+
bool SIRT::setMinMaxMasks(float* D_minMaskData_, float* D_maxMaskData_,
unsigned int iPitch)
{
diff --git a/cuda/2d/sirt.h b/cuda/2d/sirt.h
index 1dbf675..0edddab 100644
--- a/cuda/2d/sirt.h
+++ b/cuda/2d/sirt.h
@@ -41,6 +41,10 @@ public:
virtual bool init();
+ // Do optional long-object compensation. See the comments in sirt.cu.
+ // Call this after init(). It can not be used in combination with masks.
+ bool doSlabCorrections();
+
// Set min/max masks to existing GPU memory buffers
bool setMinMaxMasks(float* D_minMaskData, float* D_maxMaskData,
unsigned int pitch);
--
cgit v1.2.3