From 48f4e7b165404a0375db300b9fe59da92edf1cce Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Sat, 30 Mar 2019 20:49:42 +0100 Subject: Adjust FBP to line integral scaling --- cuda/2d/algo.cu | 13 +++++++------ cuda/2d/cgls.cu | 4 ++-- cuda/2d/fbp.cu | 5 ++--- include/astra/cuda/2d/algo.h | 7 +++++-- include/astra/cuda/2d/cgls.h | 2 +- src/CudaFilteredBackProjectionAlgorithm.cpp | 7 +++++++ src/CudaReconstructionAlgorithm2D.cpp | 4 +--- src/FilteredBackProjectionAlgorithm.cpp | 8 ++++++-- 8 files changed, 31 insertions(+), 19 deletions(-) diff --git a/cuda/2d/algo.cu b/cuda/2d/algo.cu index 11422ff..59ea984 100644 --- a/cuda/2d/algo.cu +++ b/cuda/2d/algo.cu @@ -153,6 +153,12 @@ bool ReconAlgo::setSuperSampling(int raysPerDet, int raysPerPixelDim) return true; } +bool ReconAlgo::setReconstructionScale(float fScale) +{ + fOutputScale *= fScale; + return true; +} + bool ReconAlgo::setVolumeMask(float* _D_maskData, unsigned int _maskPitch) { assert(useVolumeMask); @@ -242,7 +248,7 @@ bool ReconAlgo::allocateBuffers() return true; } -bool ReconAlgo::copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, float fSinogramScale, +bool ReconAlgo::copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, const float* pfReconstruction, unsigned int iReconstructionPitch, const float* pfVolMask, unsigned int iVolMaskPitch, const float* pfSinoMask, unsigned int iSinoMaskPitch) @@ -258,11 +264,6 @@ bool ReconAlgo::copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPit if (!ok) return false; - // rescale sinogram - if (fSinogramScale != 1.0f) - processSino(D_sinoData, fSinogramScale, - sinoPitch, dims); - ok = copyVolumeToDevice(pfReconstruction, iReconstructionPitch, dims, D_volumeData, volumePitch); diff --git a/cuda/2d/cgls.cu b/cuda/2d/cgls.cu index b6a9fae..905b960 100644 --- a/cuda/2d/cgls.cu +++ b/cuda/2d/cgls.cu @@ -102,14 +102,14 @@ bool CGLS::setBuffers(float* _D_volumeData, unsigned int _volumePitch, return true; } -bool CGLS::copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, float fSinogramScale, +bool CGLS::copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, 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); + return ReconAlgo::copyDataToGPU(pfSinogram, iSinogramPitch, pfReconstruction, iReconstructionPitch, pfVolMask, iVolMaskPitch, pfSinoMask, iSinoMaskPitch); } bool CGLS::iterate(unsigned int iterations) diff --git a/cuda/2d/fbp.cu b/cuda/2d/fbp.cu index a5b8a7a..65f1d68 100644 --- a/cuda/2d/fbp.cu +++ b/cuda/2d/fbp.cu @@ -321,15 +321,14 @@ bool FBP::iterate(unsigned int iterations) if (fanProjs) { float fOutputScale = 1.0 / (/*fPixelSize * fPixelSize * fPixelSize * */ fFanDetSize * fFanDetSize); - ok = FanBP_FBPWeighted(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, fanProjs, fOutputScale); + ok = FanBP_FBPWeighted(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, fanProjs, fOutputScale * this->fOutputScale); } else { // scale by number of angles. For the fan-beam case, this is already // handled by FDK_PreWeight float fOutputScale = (M_PI / 2.0f) / (float)dims.iProjAngles; - //fOutputScale /= fDetSize * fDetSize; - ok = BP(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, parProjs, fOutputScale); + ok = BP(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, parProjs, fOutputScale * this->fOutputScale); } if(!ok) { diff --git a/include/astra/cuda/2d/algo.h b/include/astra/cuda/2d/algo.h index 3fccdef..b6ec231 100644 --- a/include/astra/cuda/2d/algo.h +++ b/include/astra/cuda/2d/algo.h @@ -56,6 +56,10 @@ public: bool setSuperSampling(int raysPerDet, int raysPerPixelDim); + // Scale the final reconstruction. + // May be called at any time after setGeometry and before iterate(). Multiple calls stack. + bool setReconstructionScale(float fScale); + virtual bool enableVolumeMask(); virtual bool enableSinogramMask(); @@ -88,8 +92,7 @@ public: // sinogram, reconstruction, volume mask and sinogram mask in system RAM, // respectively. The corresponding pitch variables give the pitches // of these buffers, measured in floats. - // The sinogram is multiplied by fSinogramScale after uploading it. - virtual bool copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, float fSinogramScale, + virtual bool copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, const float* pfReconstruction, unsigned int iReconstructionPitch, const float* pfVolMask, unsigned int iVolMaskPitch, const float* pfSinoMask, unsigned int iSinoMaskPitch); diff --git a/include/astra/cuda/2d/cgls.h b/include/astra/cuda/2d/cgls.h index 375a425..a854a74 100644 --- a/include/astra/cuda/2d/cgls.h +++ b/include/astra/cuda/2d/cgls.h @@ -47,7 +47,7 @@ public: virtual bool setBuffers(float* D_volumeData, unsigned int volumePitch, float* D_projData, unsigned int projPitch); - virtual bool copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, float fSinogramScale, + virtual bool copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, const float* pfReconstruction, unsigned int iReconstructionPitch, const float* pfVolMask, unsigned int iVolMaskPitch, const float* pfSinoMask, unsigned int iSinoMaskPitch); diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp index 88e235b..69140fc 100644 --- a/src/CudaFilteredBackProjectionAlgorithm.cpp +++ b/src/CudaFilteredBackProjectionAlgorithm.cpp @@ -151,6 +151,13 @@ void CCudaFilteredBackProjectionAlgorithm::initCUDAAlgorithm() if (!ok) { ASTRA_ERROR("CCudaFilteredBackProjectionAlgorithm: Failed to set short-scan mode"); } + + const CVolumeGeometry2D& volGeom = *m_pProjector->getVolumeGeometry(); + float fPixelArea = volGeom.getPixelArea(); + ok &= pFBP->setReconstructionScale(1.0f/fPixelArea); + if (!ok) { + ASTRA_ERROR("CCudaFilteredBackProjectionAlgorithm: Failed to set reconstruction scale"); + } } diff --git a/src/CudaReconstructionAlgorithm2D.cpp b/src/CudaReconstructionAlgorithm2D.cpp index 939a026..6730cea 100644 --- a/src/CudaReconstructionAlgorithm2D.cpp +++ b/src/CudaReconstructionAlgorithm2D.cpp @@ -309,9 +309,7 @@ void CCudaReconstructionAlgorithm2D::run(int _iNrIterations) m_bAlgoInit = true; } - float fSinogramScale = 1.0f; - - ok = m_pAlgo->copyDataToGPU(m_pSinogram->getDataConst(), m_pSinogram->getGeometry()->getDetectorCount(), fSinogramScale, + ok = m_pAlgo->copyDataToGPU(m_pSinogram->getDataConst(), m_pSinogram->getGeometry()->getDetectorCount(), m_pReconstruction->getDataConst(), volgeom.getGridColCount(), m_bUseReconstructionMask ? m_pReconstructionMask->getDataConst() : 0, volgeom.getGridColCount(), m_bUseSinogramMask ? m_pSinogramMask->getDataConst() : 0, m_pSinogram->getGeometry()->getDetectorCount()); diff --git a/src/FilteredBackProjectionAlgorithm.cpp b/src/FilteredBackProjectionAlgorithm.cpp index 423dc6c..95bef3c 100644 --- a/src/FilteredBackProjectionAlgorithm.cpp +++ b/src/FilteredBackProjectionAlgorithm.cpp @@ -264,8 +264,12 @@ void CFilteredBackProjectionAlgorithm::run(int _iNrIterations) DefaultBPPolicy(m_pReconstruction, &filteredSinogram)); // Scale data - int iAngleCount = m_pProjector->getProjectionGeometry()->getProjectionAngleCount(); - (*m_pReconstruction) *= (PI/2)/iAngleCount; + const CVolumeGeometry2D& volGeom = *m_pProjector->getVolumeGeometry(); + const CProjectionGeometry2D& projGeom = *m_pProjector->getProjectionGeometry(); + + int iAngleCount = projGeom.getProjectionAngleCount(); + float fPixelArea = volGeom.getPixelArea(); + (*m_pReconstruction) *= PI/(2*iAngleCount*fPixelArea); m_pReconstruction->updateStatistics(); } -- cgit v1.2.3