Add relaxation parameters to SIRT, SART

4 files changed, 21 insertions, 2 deletions

diff --git a/cuda/2d/sart.cu b/cuda/2d/sart.cu
index e5cb5bb..c8608a3 100644
--- a/cuda/2d/sart.cu
+++ b/cuda/2d/sart.cu
@@ -71,6 +71,8 @@ SART::SART() : ReconAlgo()
 	projectionCount = 0;
 	iteration = 0;
 	customOrder = false;
+
+	fRelaxation = 1.0f;
 }
 
 
@@ -98,6 +100,7 @@ void SART::reset()
 	projectionCount = 0;
 	iteration = 0;
 	customOrder = false;
+	fRelaxation = 1.0f;
 
 	ReconAlgo::reset();
 }
@@ -200,10 +203,10 @@ bool SART::iterate(unsigned int iterations)
 			// BP, mask, and add back
 			// TODO: Try putting the masking directly in the BP
 			zeroVolumeData(D_tmpData, tmpPitch, dims);
-			callBP_SART(D_tmpData, tmpPitch, D_projData, projPitch, angle, 1.0f);
+			callBP_SART(D_tmpData, tmpPitch, D_projData, projPitch, angle, fRelaxation);
 			processVol<opAddMul>(D_volumeData, D_maskData, D_tmpData, volumePitch, dims);
 		} else {
-			callBP_SART(D_volumeData, volumePitch, D_projData, projPitch, angle, 1.0f);
+			callBP_SART(D_volumeData, volumePitch, D_projData, projPitch, angle, fRelaxation);
 		}
 
 		if (useMinConstraint)
diff --git a/cuda/2d/sart.h b/cuda/2d/sart.h
index 7dcd641..eff9ecf 100644
--- a/cuda/2d/sart.h
+++ b/cuda/2d/sart.h
@@ -50,6 +50,8 @@ public:
 
 	virtual float computeDiffNorm();
 
+	void setRelaxation(float r) { fRelaxation = r; }
+
 protected:
 	void reset();
 	bool precomputeWeights();
@@ -78,6 +80,8 @@ protected:
 	// Geometry-specific precomputed data
 	float* D_lineWeight;
 	unsigned int linePitch;
+
+	float fRelaxation;
 };
 
 }
diff --git a/cuda/2d/sirt.cu b/cuda/2d/sirt.cu
index 162ee77..4baaccb 100644
--- a/cuda/2d/sirt.cu
+++ b/cuda/2d/sirt.cu
@@ -50,6 +50,8 @@ SIRT::SIRT() : ReconAlgo()
 	D_minMaskData = 0;
 	D_maxMaskData = 0;
 
+	fRelaxation = 1.0f;
+
 	freeMinMaxMasks = false;
 }
 
@@ -83,6 +85,8 @@ void SIRT::reset()
 	useVolumeMask = false;
 	useSinogramMask = false;
 
+	fRelaxation = 1.0f;
+
 	ReconAlgo::reset();
 }
 
@@ -139,6 +143,9 @@ bool SIRT::precomputeWeights()
 		processVol<opMul>(D_pixelWeight, D_maskData, pixelPitch, dims);
 	}
 
+	// Also fold the relaxation factor into pixel weights
+	processVol<opMul>(D_pixelWeight, fRelaxation, pixelPitch, dims);
+
 	return true;
 }
 
@@ -253,6 +260,7 @@ bool SIRT::iterate(unsigned int iterations)
 
 		callBP(D_tmpData, tmpPitch, D_projData, projPitch, 1.0f);
 
+		// pixel weights also contain the volume mask and relaxation factor
 		processVol<opAddMul>(D_volumeData, D_pixelWeight, D_tmpData, volumePitch, dims);
 
 		if (useMinConstraint)
diff --git a/cuda/2d/sirt.h b/cuda/2d/sirt.h
index 21094a1..bc913f4 100644
--- a/cuda/2d/sirt.h
+++ b/cuda/2d/sirt.h
@@ -53,6 +53,8 @@ public:
 	bool uploadMinMaxMasks(const float* minMaskData, const float* maxMaskData,
 	                       unsigned int pitch);
 
+	void setRelaxation(float r) { fRelaxation = r; }
+
 	virtual bool iterate(unsigned int iterations);
 
 	virtual float computeDiffNorm();
@@ -81,6 +83,8 @@ protected:
 	unsigned int minMaskPitch;
 	float* D_maxMaskData;
 	unsigned int maxMaskPitch;
+
+	float fRelaxation;
 };
 
 bool doSIRT(float* D_volumeData, unsigned int volumePitch,