Update version to 1.3

1 files changed, 342 insertions, 0 deletions

diff --git a/cuda/2d/sirt.cu b/cuda/2d/sirt.cu
new file mode 100644
index 0000000..31954e4
--- /dev/null
+++ b/cuda/2d/sirt.cu
@@ -0,0 +1,342 @@
+/*
+-----------------------------------------------------------------------
+Copyright 2012 iMinds-Vision Lab, University of Antwerp
+
+Contact: astra@ua.ac.be
+Website: http://astra.ua.ac.be
+
+
+This file is part of the
+All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox").
+
+The ASTRA Toolbox is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+The ASTRA Toolbox is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
+
+-----------------------------------------------------------------------
+$Id$
+*/
+
+#include <cstdio>
+#include <cassert>
+
+#include "sirt.h"
+#include "util.h"
+#include "arith.h"
+
+#ifdef STANDALONE
+#include "testutil.h"
+#endif
+
+namespace astraCUDA {
+
+SIRT::SIRT() : ReconAlgo()
+{
+	D_projData = 0;
+	D_tmpData = 0;
+
+	D_lineWeight = 0;
+	D_pixelWeight = 0;
+
+	D_minMaskData = 0;
+	D_maxMaskData = 0;
+
+	freeMinMaxMasks = false;
+}
+
+
+SIRT::~SIRT()
+{
+	reset();
+}
+
+void SIRT::reset()
+{
+	cudaFree(D_projData);
+	cudaFree(D_tmpData);
+	cudaFree(D_lineWeight);
+	cudaFree(D_pixelWeight);
+	if (freeMinMaxMasks) {
+		cudaFree(D_minMaskData);
+		cudaFree(D_maxMaskData);
+	}
+
+	D_projData = 0;
+	D_tmpData = 0;
+
+	D_lineWeight = 0;
+	D_pixelWeight = 0;
+
+	freeMinMaxMasks = false;
+	D_minMaskData = 0;
+	D_maxMaskData = 0;
+
+	useVolumeMask = false;
+	useSinogramMask = false;
+
+	ReconAlgo::reset();
+}
+
+bool SIRT::init()
+{
+	allocateVolume(D_pixelWeight, dims.iVolWidth+2, dims.iVolHeight+2, pixelPitch);
+	zeroVolume(D_pixelWeight, pixelPitch, dims.iVolWidth+2, dims.iVolHeight+2);
+
+	allocateVolume(D_tmpData, dims.iVolWidth+2, dims.iVolHeight+2, tmpPitch);
+	zeroVolume(D_tmpData, tmpPitch, dims.iVolWidth+2, dims.iVolHeight+2);
+
+	allocateVolume(D_projData, dims.iProjDets+2, dims.iProjAngles, projPitch);
+	zeroVolume(D_projData, projPitch, dims.iProjDets+2, dims.iProjAngles);
+	
+	allocateVolume(D_lineWeight, dims.iProjDets+2, dims.iProjAngles, linePitch);
+	zeroVolume(D_lineWeight, linePitch, dims.iProjDets+2, dims.iProjAngles);
+
+	// We can't precompute lineWeights and pixelWeights when using a mask
+	if (!useVolumeMask && !useSinogramMask)
+		precomputeWeights();
+
+	// TODO: check if allocations succeeded
+	return true;
+}
+
+bool SIRT::precomputeWeights()
+{
+	zeroVolume(D_lineWeight, linePitch, dims.iProjDets+2, dims.iProjAngles);
+	if (useVolumeMask) {
+		callFP(D_maskData, maskPitch, D_lineWeight, linePitch, 1.0f);
+	} else {
+		processVol<opSet, VOL>(D_tmpData, 1.0f, tmpPitch, dims.iVolWidth, dims.iVolHeight);
+		callFP(D_tmpData, tmpPitch, D_lineWeight, linePitch, 1.0f);
+	}
+	processVol<opInvert, SINO>(D_lineWeight, linePitch, dims.iProjDets, dims.iProjAngles);
+
+	if (useSinogramMask) {
+		// scale line weights with sinogram mask to zero out masked sinogram pixels
+		processVol<opMul, SINO>(D_lineWeight, D_smaskData, linePitch, dims.iProjDets, dims.iProjAngles);
+	}
+
+
+	zeroVolume(D_pixelWeight, pixelPitch, dims.iVolWidth+2, dims.iVolHeight+2);
+	if (useSinogramMask) {
+		callBP(D_pixelWeight, pixelPitch, D_smaskData, smaskPitch);
+	} else {
+		processVol<opSet, SINO>(D_projData, 1.0f, projPitch, dims.iProjDets, dims.iProjAngles);
+		callBP(D_pixelWeight, pixelPitch, D_projData, projPitch);
+	}
+	processVol<opInvert, VOL>(D_pixelWeight, pixelPitch, dims.iVolWidth, dims.iVolHeight);
+
+	if (useVolumeMask) {
+		// scale pixel weights with mask to zero out masked pixels
+		processVol<opMul, VOL>(D_pixelWeight, D_maskData, pixelPitch, dims.iVolWidth, dims.iVolHeight);
+	}
+
+	return true;
+}
+
+bool SIRT::setMinMaxMasks(float* D_minMaskData_, float* D_maxMaskData_,
+	                      unsigned int iPitch)
+{
+	D_minMaskData = D_minMaskData_;
+	D_maxMaskData = D_maxMaskData_;
+	minMaskPitch = iPitch;
+	maxMaskPitch = iPitch;
+
+	freeMinMaxMasks = false;
+	return true;
+}
+
+bool SIRT::uploadMinMaxMasks(const float* pfMinMaskData, const float* pfMaxMaskData,
+	                         unsigned int iPitch)
+{
+	freeMinMaxMasks = true;
+	bool ok = true;
+	if (pfMinMaskData) {
+		allocateVolume(D_minMaskData, dims.iVolWidth+2, dims.iVolHeight+2, minMaskPitch);
+		ok = copyVolumeToDevice(pfMinMaskData, iPitch,
+		                        dims.iVolWidth, dims.iVolHeight,
+		                        D_minMaskData, minMaskPitch);
+	}
+	if (!ok)
+		return false;
+
+	if (pfMaxMaskData) {
+		allocateVolume(D_maxMaskData, dims.iVolWidth+2, dims.iVolHeight+2, maxMaskPitch);
+		ok = copyVolumeToDevice(pfMaxMaskData, iPitch,
+		                        dims.iVolWidth, dims.iVolHeight,
+		                        D_maxMaskData, maxMaskPitch);
+	}
+	if (!ok)
+		return false;
+
+	return true;
+}
+
+bool SIRT::iterate(unsigned int iterations)
+{
+	shouldAbort = false;
+
+	if (useVolumeMask || useSinogramMask)
+		precomputeWeights();
+
+	// iteration
+	for (unsigned int iter = 0; iter < iterations && !shouldAbort; ++iter) {
+
+		// copy sinogram to projection data
+		cudaMemcpy2D(D_projData, sizeof(float)*projPitch, D_sinoData, sizeof(float)*sinoPitch, sizeof(float)*(dims.iProjDets+2), dims.iProjAngles, cudaMemcpyDeviceToDevice);
+
+		// do FP, subtracting projection from sinogram
+		if (useVolumeMask) {
+				cudaMemcpy2D(D_tmpData, sizeof(float)*tmpPitch, D_volumeData, sizeof(float)*volumePitch, sizeof(float)*(dims.iVolWidth+2), dims.iVolHeight+2, cudaMemcpyDeviceToDevice);
+				processVol<opMul, VOL>(D_tmpData, D_maskData, tmpPitch, dims.iVolWidth, dims.iVolHeight);
+				callFP(D_tmpData, tmpPitch, D_projData, projPitch, -1.0f);
+		} else {
+				callFP(D_volumeData, volumePitch, D_projData, projPitch, -1.0f);
+		}
+
+		processVol<opMul, SINO>(D_projData, D_lineWeight, projPitch, dims.iProjDets, dims.iProjAngles);
+
+		zeroVolume(D_tmpData, tmpPitch, dims.iVolWidth+2, dims.iVolHeight+2);
+
+		callBP(D_tmpData, tmpPitch, D_projData, projPitch);
+
+		processVol<opAddMul, VOL>(D_volumeData, D_pixelWeight, D_tmpData, volumePitch, dims.iVolWidth, dims.iVolHeight);
+
+		if (useMinConstraint)
+			processVol<opClampMin, VOL>(D_volumeData, fMinConstraint, volumePitch, dims.iVolWidth, dims.iVolHeight);
+		if (useMaxConstraint)
+			processVol<opClampMax, VOL>(D_volumeData, fMaxConstraint, volumePitch, dims.iVolWidth, dims.iVolHeight);
+		if (D_minMaskData)
+			processVol<opClampMinMask, VOL>(D_volumeData, D_minMaskData, volumePitch, dims.iVolWidth, dims.iVolHeight);
+		if (D_maxMaskData)
+			processVol<opClampMaxMask, VOL>(D_volumeData, D_maxMaskData, volumePitch, dims.iVolWidth, dims.iVolHeight);
+	}
+
+	return true;
+}
+
+float SIRT::computeDiffNorm()
+{
+	// copy sinogram to projection data
+	cudaMemcpy2D(D_projData, sizeof(float)*projPitch, D_sinoData, sizeof(float)*sinoPitch, sizeof(float)*(dims.iProjDets+2), dims.iProjAngles, cudaMemcpyDeviceToDevice);
+
+	// do FP, subtracting projection from sinogram
+	if (useVolumeMask) {
+			cudaMemcpy2D(D_tmpData, sizeof(float)*tmpPitch, D_volumeData, sizeof(float)*volumePitch, sizeof(float)*(dims.iVolWidth+2), dims.iVolHeight+2, cudaMemcpyDeviceToDevice);
+			processVol<opMul, VOL>(D_tmpData, D_maskData, tmpPitch, dims.iVolWidth, dims.iVolHeight);
+			callFP(D_tmpData, tmpPitch, D_projData, projPitch, -1.0f);
+	} else {
+			callFP(D_volumeData, volumePitch, D_projData, projPitch, -1.0f);
+	}
+
+
+	// compute norm of D_projData
+
+	float s = dotProduct2D(D_projData, projPitch, dims.iProjDets, dims.iProjAngles, 1, 0);
+
+	return sqrt(s);
+}
+
+
+bool doSIRT(float* D_volumeData, unsigned int volumePitch,
+            float* D_sinoData, unsigned int sinoPitch,
+            float* D_maskData, unsigned int maskPitch,
+            const SDimensions& dims, const float* angles,
+            const float* TOffsets, unsigned int iterations)
+{
+	SIRT sirt;
+	bool ok = true;
+
+	ok &= sirt.setGeometry(dims, angles);
+	if (D_maskData)
+		ok &= sirt.enableVolumeMask();
+	if (TOffsets)
+		ok &= sirt.setTOffsets(TOffsets);
+
+	if (!ok)
+		return false;
+
+	ok = sirt.init();
+	if (!ok)
+		return false;
+
+	if (D_maskData)
+		ok &= sirt.setVolumeMask(D_maskData, maskPitch);
+
+	ok &= sirt.setBuffers(D_volumeData, volumePitch, D_sinoData, sinoPitch);
+	if (!ok)
+		return false;
+
+	ok = sirt.iterate(iterations);
+
+	return ok;
+}
+
+}
+
+#ifdef STANDALONE
+
+using namespace astraCUDA;
+
+int main()
+{
+	float* D_volumeData;
+	float* D_sinoData;
+
+	SDimensions dims;
+	dims.iVolWidth = 1024;
+	dims.iVolHeight = 1024;
+	dims.iProjAngles = 512;
+	dims.iProjDets = 1536;
+	dims.fDetScale = 1.0f;
+	dims.iRaysPerDet = 1;
+	unsigned int volumePitch, sinoPitch;
+
+	allocateVolume(D_volumeData, dims.iVolWidth+2, dims.iVolHeight+2, volumePitch);
+	zeroVolume(D_volumeData, volumePitch, dims.iVolWidth+2, dims.iVolHeight+2);
+	printf("pitch: %u\n", volumePitch);
+
+	allocateVolume(D_sinoData, dims.iProjDets+2, dims.iProjAngles, sinoPitch);
+	zeroVolume(D_sinoData, sinoPitch, dims.iProjDets+2, dims.iProjAngles);
+	printf("pitch: %u\n", sinoPitch);
+	
+	unsigned int y, x;
+	float* sino = loadImage("sino.png", y, x);
+
+	float* img = new float[dims.iVolWidth*dims.iVolHeight];
+
+	copySinogramToDevice(sino, dims.iProjDets, dims.iProjDets, dims.iProjAngles, D_sinoData, sinoPitch);
+
+	float* angle = new float[dims.iProjAngles];
+
+	for (unsigned int i = 0; i < dims.iProjAngles; ++i)
+		angle[i] = i*(M_PI/dims.iProjAngles);
+
+	SIRT sirt;
+
+	sirt.setGeometry(dims, angle);
+	sirt.init();
+
+	sirt.setBuffers(D_volumeData, volumePitch, D_sinoData, sinoPitch);
+
+	sirt.iterate(25);
+
+
+	delete[] angle;
+
+	copyVolumeFromDevice(img, dims.iVolWidth, dims.iVolWidth, dims.iVolHeight, D_volumeData, volumePitch);
+
+	saveImage("vol.png",dims.iVolHeight,dims.iVolWidth,img);
+
+	return 0;
+}
+#endif
+