From 0c77eee16e2f4161c1ebc110b15ce6563d4a96c2 Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <WillemJan.Palenstijn@uantwerpen.be>
Date: Wed, 16 Apr 2014 11:13:22 +0000
Subject: Add fan beam support to FBP_CUDA
---
cuda/2d/astra.cu | 91 ++++++++++++++++++++++++++++-
cuda/2d/astra.h | 13 +++++
cuda/3d/fdk.cu | 23 +++++++-
cuda/3d/fdk.h | 8 +++
src/CudaFilteredBackProjectionAlgorithm.cpp | 30 ++++++++--
5 files changed, 155 insertions(+), 10 deletions(-)
diff --git a/cuda/2d/astra.cu b/cuda/2d/astra.cu
index 1c2e623..e317f6c 100644
--- a/cuda/2d/astra.cu
+++ b/cuda/2d/astra.cu
@@ -33,6 +33,7 @@ $Id$
#include "par_fp.h"
#include "fan_fp.h"
#include "par_bp.h"
+#include "fan_bp.h"
#include "arith.h"
#include "astra.h"
@@ -43,6 +44,9 @@ $Id$
#include "../../include/astra/Logger.h"
+// For fan beam FBP weighting
+#include "../3d/fdk.h"
+
using namespace astraCUDA;
using namespace std;
@@ -61,8 +65,14 @@ public:
float* angles;
float* TOffsets;
+ float fOriginSourceDistance;
+ float fOriginDetectorDistance;
+
float fPixelSize;
+ bool bFanBeam;
+ bool bShortScan;
+
bool initialized;
bool setStartReconstruction;
@@ -152,9 +162,33 @@ bool AstraFBP::setProjectionGeometry(unsigned int iProjAngles,
pData->angles = new float[iProjAngles];
memcpy(pData->angles, pfAngles, iProjAngles * sizeof(pfAngles[0]));
+ pData->bFanBeam = false;
+
+ return true;
+}
+
+bool AstraFBP::setFanGeometry(unsigned int iProjAngles,
+ unsigned int iProjDets,
+ const float* pfAngles,
+ float fOriginSourceDistance,
+ float fOriginDetectorDistance,
+ float fDetSize,
+ bool bShortScan)
+{
+ // Slightly abusing setProjectionGeometry for this...
+ if (!setProjectionGeometry(iProjAngles, iProjDets, pfAngles, fDetSize))
+ return false;
+
+ pData->fOriginSourceDistance = fOriginSourceDistance;
+ pData->fOriginDetectorDistance = fOriginDetectorDistance;
+
+ pData->bFanBeam = true;
+ pData->bShortScan = bShortScan;
+
return true;
}
+
bool AstraFBP::setPixelSuperSampling(unsigned int iPixelSuperSampling)
{
if (pData->initialized)
@@ -274,6 +308,34 @@ bool AstraFBP::run()
bool ok = false;
+ if (pData->bFanBeam) {
+ // Call FDK_PreWeight to handle fan beam geometry. We treat
+ // this as a cone beam setup of a single slice:
+
+ // TODO: TOffsets...
+
+ // We create a fake cudaPitchedPtr
+ cudaPitchedPtr tmp;
+ tmp.ptr = pData->D_sinoData;
+ tmp.pitch = pData->sinoPitch * sizeof(float);
+ tmp.xsize = pData->dims.iProjDets;
+ tmp.ysize = pData->dims.iProjAngles;
+ // and a fake Dimensions3D
+ astraCUDA3d::SDimensions3D dims3d;
+ dims3d.iVolX = pData->dims.iVolWidth;
+ dims3d.iVolY = pData->dims.iVolHeight;
+ dims3d.iVolZ = 1;
+ dims3d.iProjAngles = pData->dims.iProjAngles;
+ dims3d.iProjU = pData->dims.iProjDets;
+ dims3d.iProjV = 1;
+ dims3d.iRaysPerDetDim = dims3d.iRaysPerVoxelDim = 1;
+
+ astraCUDA3d::FDK_PreWeight(tmp, pData->fOriginSourceDistance,
+ pData->fOriginDetectorDistance, 0.0f, 0.0f,
+ pData->dims.fDetScale, 1.0f, // TODO: Are these correct?
+ pData->bShortScan, dims3d, pData->angles);
+ }
+
if (pData->m_pDevFilter) {
int iFFTRealDetCount = calcNextPowerOfTwo(2 * pData->dims.iProjDets);
@@ -293,7 +355,34 @@ bool AstraFBP::run()
}
- ok = BP(pData->D_volumeData, pData->volumePitch, pData->D_sinoData, pData->sinoPitch, pData->dims, pData->angles, pData->TOffsets);
+ if (pData->bFanBeam) {
+ // TODO: TOffsets?
+ // TODO: Remove this code duplication with CudaReconstructionAlgorithm
+ SFanProjection* projs;
+ projs = new SFanProjection[pData->dims.iProjAngles];
+
+ float fSrcX0 = 0.0f;
+ float fSrcY0 = -pData->fOriginSourceDistance / pData->fPixelSize;
+ float fDetUX0 = pData->dims.fDetScale;
+ float fDetUY0 = 0.0f;
+ float fDetSX0 = pData->dims.iProjDets * fDetUX0 / -2.0f;
+ float fDetSY0 = pData->fOriginDetectorDistance / pData->fPixelSize;
+
+#define ROTATE0(name,i,alpha) do { projs[i].f##name##X = f##name##X0 * cos(alpha) - f##name##Y0 * sin(alpha); projs[i].f##name##Y = f##name##X0 * sin(alpha) + f##name##Y0 * cos(alpha); } while(0)
+ for (unsigned int i = 0; i < pData->dims.iProjAngles; ++i) {
+ ROTATE0(Src, i, pData->angles[i]);
+ ROTATE0(DetS, i, pData->angles[i]);
+ ROTATE0(DetU, i, pData->angles[i]);
+ }
+
+#undef ROTATE0
+ ok = FanBP_FBPWeighted(pData->D_volumeData, pData->volumePitch, pData->D_sinoData, pData->sinoPitch, pData->dims, projs);
+
+ delete[] projs;
+
+ } else {
+ ok = BP(pData->D_volumeData, pData->volumePitch, pData->D_sinoData, pData->sinoPitch, pData->dims, pData->angles, pData->TOffsets);
+ }
if(!ok)
{
return false;
diff --git a/cuda/2d/astra.h b/cuda/2d/astra.h
index 9e58301..42b15e8 100644
--- a/cuda/2d/astra.h
+++ b/cuda/2d/astra.h
@@ -68,6 +68,19 @@ public:
unsigned int iProjDets,
const float *pfAngles,
float fDetSize = 1.0f);
+ // Set the projection angles and number of detector pixels per angle.
+ // pfAngles must be a float array of length iProjAngles.
+ // fDetSize indicates the size of a detector pixel compared to a
+ // volume pixel edge.
+ //
+ // pfAngles will only be read from during this call.
+ bool setFanGeometry(unsigned int iProjAngles,
+ unsigned int iProjDets,
+ const float *pfAngles,
+ float fOriginSourceDistance,
+ float fOriginDetectorDistance,
+ float fDetSize = 1.0f,
+ bool bShortScan = false);
// Set linear supersampling factor for the BP.
// (The number of rays is the square of this)
diff --git a/cuda/3d/fdk.cu b/cuda/3d/fdk.cu
index 883187d..fd5a8a2 100644
--- a/cuda/3d/fdk.cu
+++ b/cuda/3d/fdk.cu
@@ -298,8 +298,7 @@ __global__ void devFDK_ParkerWeight(void* D_projData, unsigned int projPitch, un
// Perform the FDK pre-weighting and filtering
-bool FDK_Filter(cudaPitchedPtr D_projData,
- cufftComplex * D_filter,
+bool FDK_PreWeight(cudaPitchedPtr D_projData,
float fSrcOrigin, float fDetOrigin,
float fSrcZ, float fDetZ,
float fDetUSize, float fDetVSize, bool bShortScan,
@@ -335,13 +334,22 @@ bool FDK_Filter(cudaPitchedPtr D_projData,
}
cudaTextForceKernelsCompletion();
+ return true;
+}
+bool FDK_Filter(cudaPitchedPtr D_projData,
+ cufftComplex * D_filter,
+ float fSrcOrigin, float fDetOrigin,
+ float fSrcZ, float fDetZ,
+ float fDetUSize, float fDetVSize, bool bShortScan,
+ const SDimensions3D& dims, const float* angles)
+{
// The filtering is a regular ramp filter per detector line.
int iPaddedDetCount = calcNextPowerOfTwo(2 * dims.iProjU);
int iHalfFFTSize = calcFFTFourSize(iPaddedDetCount);
-
+ int projPitch = D_projData.pitch/sizeof(float);
// We process one sinogram at a time.
@@ -390,11 +398,18 @@ bool FDK(cudaPitchedPtr D_volumeData,
int iPaddedDetCount = calcNextPowerOfTwo(2 * dims.iProjU);
int iHalfFFTSize = calcFFTFourSize(iPaddedDetCount);
+ ok = FDK_PreWeight(D_projData, fSrcOrigin, fDetOrigin,
+ fSrcZ, fDetZ, fDetUSize, fDetVSize,
+ bShortScan, dims, angles);
+ if (!ok)
+ return false;
+
cufftComplex *pHostFilter = new cufftComplex[dims.iProjAngles * iHalfFFTSize];
memset(pHostFilter, 0, sizeof(cufftComplex) * dims.iProjAngles * iHalfFFTSize);
genFilter(FILTER_RAMLAK, 1.0f, dims.iProjAngles, pHostFilter, iPaddedDetCount, iHalfFFTSize);
+
allocateComplexOnDevice(dims.iProjAngles, iHalfFFTSize, &D_filter);
uploadComplexArrayToDevice(dims.iProjAngles, iHalfFFTSize, pHostFilter, D_filter);
@@ -403,6 +418,8 @@ bool FDK(cudaPitchedPtr D_volumeData,
// Perform filtering
+
+
ok = FDK_Filter(D_projData, D_filter, fSrcOrigin, fDetOrigin,
fSrcZ, fDetZ, fDetUSize, fDetVSize,
bShortScan, dims, angles);
diff --git a/cuda/3d/fdk.h b/cuda/3d/fdk.h
index 5443b19..c9123b9 100644
--- a/cuda/3d/fdk.h
+++ b/cuda/3d/fdk.h
@@ -29,8 +29,16 @@ $Id$
#ifndef _CUDA_FDK_H
#define _CUDA_FDK_H
+#include "dims3d.h"
+
namespace astraCUDA3d {
+bool FDK_PreWeight(cudaPitchedPtr D_projData,
+ float fSrcOrigin, float fDetOrigin,
+ float fSrcZ, float fDetZ,
+ float fDetUSize, float fDetVSize, bool bShortScan,
+ const SDimensions3D& dims, const float* angles);
+
bool FDK(cudaPitchedPtr D_volumeData,
cudaPitchedPtr D_projData,
float fSrcOrigin, float fDetOrigin,
diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp
index 3656f3f..c53ac2d 100644
--- a/src/CudaFilteredBackProjectionAlgorithm.cpp
+++ b/src/CudaFilteredBackProjectionAlgorithm.cpp
@@ -27,6 +27,7 @@ $Id$
*/
#include <astra/CudaFilteredBackProjectionAlgorithm.h>
+#include <astra/FanFlatProjectionGeometry2D.h>
#include <boost/lexical_cast.hpp>
#include <cstring>
@@ -226,7 +227,8 @@ void CCudaFilteredBackProjectionAlgorithm::run(int _iNrIterations /* = 0 */)
if (!m_bAstraFBPInit) {
const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry();
- const CParallelProjectionGeometry2D& projgeom = *dynamic_cast<CParallelProjectionGeometry2D*>(m_pSinogram->getGeometry());
+ const CParallelProjectionGeometry2D* parprojgeom = dynamic_cast<CParallelProjectionGeometry2D*>(m_pSinogram->getGeometry());
+ const CFanFlatProjectionGeometry2D* fanprojgeom = dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pSinogram->getGeometry());
bool ok = true;
@@ -235,10 +237,26 @@ void CCudaFilteredBackProjectionAlgorithm::run(int _iNrIterations /* = 0 */)
volgeom.getGridRowCount(),
volgeom.getPixelLengthX());
// TODO: off-center geometry
- ok &= m_pFBP->setProjectionGeometry(projgeom.getProjectionAngleCount(),
- projgeom.getDetectorCount(),
- projgeom.getProjectionAngles(),
- projgeom.getDetectorWidth());
+ int iDetectorCount;
+ if (parprojgeom) {
+ ok &= m_pFBP->setProjectionGeometry(parprojgeom->getProjectionAngleCount(),
+ parprojgeom->getDetectorCount(),
+ parprojgeom->getProjectionAngles(),
+ parprojgeom->getDetectorWidth());
+ iDetectorCount = parprojgeom->getDetectorCount();
+ } else if (fanprojgeom) {
+ ok &= m_pFBP->setFanGeometry(fanprojgeom->getProjectionAngleCount(),
+ fanprojgeom->getDetectorCount(),
+ fanprojgeom->getProjectionAngles(),
+ fanprojgeom->getOriginSourceDistance(),
+ fanprojgeom->getOriginDetectorDistance(),
+ fanprojgeom->getDetectorWidth(),
+ false); // TODO: Support short-scan
+
+ iDetectorCount = fanprojgeom->getDetectorCount();
+ } else {
+ assert(false);
+ }
ok &= m_pFBP->setPixelSuperSampling(m_iPixelSuperSampling);
@@ -252,7 +270,7 @@ void CCudaFilteredBackProjectionAlgorithm::run(int _iNrIterations /* = 0 */)
ok &= m_pFBP->init(m_iGPUIndex);
ASTRA_ASSERT(ok);
- ok &= m_pFBP->setSinogram(m_pSinogram->getDataConst(), projgeom.getDetectorCount());
+ ok &= m_pFBP->setSinogram(m_pSinogram->getDataConst(), iDetectorCount);
ASTRA_ASSERT(ok);
ok &= m_pFBP->setFilter(m_eFilter, m_pfFilter, m_iFilterWidth, m_fFilterD, m_fFilterParameter);
--
cgit v1.2.3