4 files changed, 59 insertions, 0 deletions

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);