summaryrefslogtreecommitdiffstats
path: root/cuda/2d/cgls.cu
diff options
context:
space:
mode:
Diffstat (limited to 'cuda/2d/cgls.cu')
-rw-r--r--cuda/2d/cgls.cu304
1 files changed, 304 insertions, 0 deletions
diff --git a/cuda/2d/cgls.cu b/cuda/2d/cgls.cu
new file mode 100644
index 0000000..5b1cf46
--- /dev/null
+++ b/cuda/2d/cgls.cu
@@ -0,0 +1,304 @@
+/*
+-----------------------------------------------------------------------
+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 "cgls.h"
+#include "util.h"
+#include "arith.h"
+
+#ifdef STANDALONE
+#include "testutil.h"
+#endif
+
+namespace astraCUDA {
+
+CGLS::CGLS() : ReconAlgo()
+{
+ D_z = 0;
+ D_p = 0;
+ D_r = 0;
+ D_w = 0;
+
+ sliceInitialized = false;
+}
+
+
+CGLS::~CGLS()
+{
+ reset();
+}
+
+void CGLS::reset()
+{
+ cudaFree(D_z);
+ cudaFree(D_p);
+ cudaFree(D_r);
+ cudaFree(D_w);
+
+ D_z = 0;
+ D_p = 0;
+ D_r = 0;
+ D_w = 0;
+
+ ReconAlgo::reset();
+}
+
+bool CGLS::init()
+{
+ // Lifetime of z: within an iteration
+ allocateVolume(D_z, dims.iVolWidth+2, dims.iVolHeight+2, zPitch);
+
+ // Lifetime of p: full algorithm
+ allocateVolume(D_p, dims.iVolWidth+2, dims.iVolHeight+2, pPitch);
+
+ // Lifetime of r: full algorithm
+ allocateVolume(D_r, dims.iProjDets+2, dims.iProjAngles, rPitch);
+
+ // Lifetime of w: within an iteration
+ allocateVolume(D_w, dims.iProjDets+2, dims.iProjAngles, wPitch);
+
+ // TODO: check if allocations succeeded
+ return true;
+}
+
+
+bool CGLS::setBuffers(float* _D_volumeData, unsigned int _volumePitch,
+ float* _D_projData, unsigned int _projPitch)
+{
+ bool ok = ReconAlgo::setBuffers(_D_volumeData, _volumePitch,
+ _D_projData, _projPitch);
+
+ if (!ok)
+ return false;
+
+ sliceInitialized = false;
+
+ return true;
+}
+
+bool CGLS::copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, float fSinogramScale,
+ const float* pfReconstruction, unsigned int iReconstructionPitch,
+ const float* pfVolMask, unsigned int iVolMaskPitch,
+ const float* pfSinoMask, unsigned int iSinoMaskPitch)
+{
+ sliceInitialized = false;
+
+ return ReconAlgo::copyDataToGPU(pfSinogram, iSinogramPitch, fSinogramScale, pfReconstruction, iReconstructionPitch, pfVolMask, iVolMaskPitch, pfSinoMask, iSinoMaskPitch);
+}
+
+bool CGLS::iterate(unsigned int iterations)
+{
+ shouldAbort = false;
+
+ if (!sliceInitialized) {
+
+ // copy sinogram
+ cudaMemcpy2D(D_r, sizeof(float)*rPitch, D_sinoData, sizeof(float)*sinoPitch, sizeof(float)*(dims.iProjDets+2), dims.iProjAngles, cudaMemcpyDeviceToDevice);
+
+ // r = sino - A*x
+ if (useVolumeMask) {
+ // Use z as temporary storage here since it is unused
+ cudaMemcpy2D(D_z, sizeof(float)*zPitch, D_volumeData, sizeof(float)*volumePitch, sizeof(float)*(dims.iVolWidth+2), dims.iVolHeight+2, cudaMemcpyDeviceToDevice);
+ processVol<opMul, VOL>(D_z, D_maskData, zPitch, dims.iVolWidth, dims.iVolHeight);
+ callFP(D_z, zPitch, D_r, rPitch, -1.0f);
+ } else {
+ callFP(D_volumeData, volumePitch, D_r, rPitch, -1.0f);
+ }
+
+
+ // p = A'*r
+ zeroVolume(D_p, pPitch, dims.iVolWidth+2, dims.iVolHeight+2);
+ callBP(D_p, pPitch, D_r, rPitch);
+ if (useVolumeMask)
+ processVol<opMul, VOL>(D_p, D_maskData, pPitch, dims.iVolWidth, dims.iVolHeight);
+
+
+ gamma = dotProduct2D(D_p, pPitch, dims.iVolWidth, dims.iVolHeight, 1, 1);
+
+ sliceInitialized = true;
+ }
+
+
+ // iteration
+ for (unsigned int iter = 0; iter < iterations && !shouldAbort; ++iter) {
+
+ // w = A*p
+ zeroVolume(D_w, wPitch, dims.iProjDets+2, dims.iProjAngles);
+ callFP(D_p, pPitch, D_w, wPitch, 1.0f);
+
+ // alpha = gamma / <w,w>
+ float ww = dotProduct2D(D_w, wPitch, dims.iProjDets, dims.iProjAngles, 1, 0);
+ float alpha = gamma / ww;
+
+ // x += alpha*p
+ processVol<opAddScaled, VOL>(D_volumeData, D_p, alpha, volumePitch, dims.iVolWidth, dims.iVolHeight);
+
+ // r -= alpha*w
+ processVol<opAddScaled, SINO>(D_r, D_w, -alpha, rPitch, dims.iProjDets, dims.iProjAngles);
+
+
+ // z = A'*r
+ zeroVolume(D_z, zPitch, dims.iVolWidth+2, dims.iVolHeight+2);
+ callBP(D_z, zPitch, D_r, rPitch);
+ if (useVolumeMask)
+ processVol<opMul, VOL>(D_z, D_maskData, zPitch, dims.iVolWidth, dims.iVolHeight);
+
+ float beta = 1.0f / gamma;
+ gamma = dotProduct2D(D_z, zPitch, dims.iVolWidth, dims.iVolHeight, 1, 1);
+ beta *= gamma;
+
+ // p = z + beta*p
+ processVol<opScaleAndAdd, VOL>(D_p, D_z, beta, pPitch, dims.iVolWidth, dims.iVolHeight);
+
+ }
+
+ return true;
+}
+
+
+float CGLS::computeDiffNorm()
+{
+ // We can use w and z as temporary storage here since they're not
+ // used outside of iterations.
+
+ // copy sinogram to w
+ cudaMemcpy2D(D_w, sizeof(float)*wPitch, D_sinoData, sizeof(float)*sinoPitch, sizeof(float)*(dims.iProjDets+2), dims.iProjAngles, cudaMemcpyDeviceToDevice);
+
+ // do FP, subtracting projection from sinogram
+ if (useVolumeMask) {
+ cudaMemcpy2D(D_z, sizeof(float)*zPitch, D_volumeData, sizeof(float)*volumePitch, sizeof(float)*(dims.iVolWidth+2), dims.iVolHeight+2, cudaMemcpyDeviceToDevice);
+ processVol<opMul, VOL>(D_z, D_maskData, zPitch, dims.iVolWidth, dims.iVolHeight);
+ callFP(D_z, zPitch, D_w, wPitch, -1.0f);
+ } else {
+ callFP(D_volumeData, volumePitch, D_w, wPitch, -1.0f);
+ }
+
+ // compute norm of D_w
+
+ float s = dotProduct2D(D_w, wPitch, dims.iProjDets, dims.iProjAngles, 1, 0);
+
+ return sqrt(s);
+}
+
+bool doCGLS(float* D_volumeData, unsigned int volumePitch,
+ float* D_sinoData, unsigned int sinoPitch,
+ const SDimensions& dims, /*const SAugmentedData& augs,*/
+ const float* angles, const float* TOffsets, unsigned int iterations)
+{
+ CGLS cgls;
+ bool ok = true;
+
+ ok &= cgls.setGeometry(dims, angles);
+#if 0
+ if (D_maskData)
+ ok &= cgls.enableVolumeMask();
+#endif
+ if (TOffsets)
+ ok &= cgls.setTOffsets(TOffsets);
+
+ if (!ok)
+ return false;
+
+ ok = cgls.init();
+ if (!ok)
+ return false;
+
+#if 0
+ if (D_maskData)
+ ok &= cgls.setVolumeMask(D_maskData, maskPitch);
+#endif
+
+ ok &= cgls.setBuffers(D_volumeData, volumePitch, D_sinoData, sinoPitch);
+ if (!ok)
+ return false;
+
+ ok = cgls.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);
+
+ CGLS cgls;
+
+ cgls.setGeometry(dims, angle);
+ cgls.init();
+
+ cgls.setBuffers(D_volumeData, volumePitch, D_sinoData, sinoPitch);
+
+ cgls.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