summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorWillem Jan Palenstijn <wjp@usecode.org>2019-03-30 20:49:42 +0100
committerWillem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>2019-09-25 14:10:08 +0200
commit48f4e7b165404a0375db300b9fe59da92edf1cce (patch)
tree493a6f578d6b51baf0db79f898a5ed3e9a2074fc
parent3cf63d335ebe392a8c77f0c90395c18150647aeb (diff)
downloadastra-48f4e7b165404a0375db300b9fe59da92edf1cce.tar.gz
astra-48f4e7b165404a0375db300b9fe59da92edf1cce.tar.bz2
astra-48f4e7b165404a0375db300b9fe59da92edf1cce.tar.xz
astra-48f4e7b165404a0375db300b9fe59da92edf1cce.zip
Adjust FBP to line integral scaling
-rw-r--r--cuda/2d/algo.cu13
-rw-r--r--cuda/2d/cgls.cu4
-rw-r--r--cuda/2d/fbp.cu5
-rw-r--r--include/astra/cuda/2d/algo.h7
-rw-r--r--include/astra/cuda/2d/cgls.h2
-rw-r--r--src/CudaFilteredBackProjectionAlgorithm.cpp7
-rw-r--r--src/CudaReconstructionAlgorithm2D.cpp4
-rw-r--r--src/FilteredBackProjectionAlgorithm.cpp8
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<opMul>(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();
}