summaryrefslogtreecommitdiffstats
path: root/cuda
diff options
context:
space:
mode:
authorWillem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>2014-10-02 13:50:59 +0200
committerWillem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>2014-10-02 13:51:54 +0200
commit12fd6925502d3ead9b73498f44f9b2f9abbc9bea (patch)
tree10970caff57e8f6b127dc223f73270eb3a490dc8 /cuda
parent8188b0697751d7c7ec6fe3069d0fceb50ebd0c9e (diff)
downloadastra-12fd6925502d3ead9b73498f44f9b2f9abbc9bea.tar.gz
astra-12fd6925502d3ead9b73498f44f9b2f9abbc9bea.tar.bz2
astra-12fd6925502d3ead9b73498f44f9b2f9abbc9bea.tar.xz
astra-12fd6925502d3ead9b73498f44f9b2f9abbc9bea.zip
Add CUDA SIRT::doSlabCorrections() function
This function optionally compensates for effectively infinitely large slab-like objects of finite thickness 1.
Diffstat (limited to 'cuda')
-rw-r--r--cuda/2d/arith.cu9
-rw-r--r--cuda/2d/arith.h1
-rw-r--r--cuda/2d/sirt.cu45
-rw-r--r--cuda/2d/sirt.h4
4 files changed, 59 insertions, 0 deletions
diff --git a/cuda/2d/arith.cu b/cuda/2d/arith.cu
index 04d4de9..93b1ee0 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 f730e2f..5a7d034 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 d34a180..98c7f09 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 5592616..d44a20a 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);