From 5d5f27a910d1fa63974cf4f4dc2485c60906441c Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Mon, 6 Oct 2014 11:46:53 +0200 Subject: Start supporting flexible 2D volume geometry in Cuda FBP too --- cuda/2d/astra.cu | 36 ++++++----------- cuda/2d/astra.h | 3 +- src/CudaFilteredBackProjectionAlgorithm.cpp | 61 ++++++++++++++++++++++------- 3 files changed, 61 insertions(+), 39 deletions(-) diff --git a/cuda/2d/astra.cu b/cuda/2d/astra.cu index 087905d..be273eb 100644 --- a/cuda/2d/astra.cu +++ b/cuda/2d/astra.cu @@ -68,6 +68,7 @@ public: SDimensions dims; float* angles; float* TOffsets; + astraCUDA::SFanProjection* fanProjections; float fOriginSourceDistance; float fOriginDetectorDistance; @@ -94,6 +95,8 @@ AstraFBP::AstraFBP() pData = new AstraFBP_internal(); pData->angles = 0; + pData->fanProjections = 0; + pData->TOffsets = 0; pData->D_sinoData = 0; pData->D_volumeData = 0; @@ -117,6 +120,9 @@ AstraFBP::~AstraFBP() delete[] pData->TOffsets; pData->TOffsets = 0; + delete[] pData->fanProjections; + pData->fanProjections = 0; + cudaFree(pData->D_sinoData); pData->D_sinoData = 0; @@ -173,6 +179,7 @@ bool AstraFBP::setProjectionGeometry(unsigned int iProjAngles, bool AstraFBP::setFanGeometry(unsigned int iProjAngles, unsigned int iProjDets, + const astraCUDA::SFanProjection *fanProjs, const float* pfAngles, float fOriginSourceDistance, float fOriginDetectorDistance, @@ -186,6 +193,9 @@ bool AstraFBP::setFanGeometry(unsigned int iProjAngles, pData->fOriginSourceDistance = fOriginSourceDistance; pData->fOriginDetectorDistance = fOriginDetectorDistance; + pData->fanProjections = new astraCUDA::SFanProjection[iProjAngles]; + memcpy(pData->fanProjections, fanProjs, iProjAngles * sizeof(fanProjs[0])); + pData->bFanBeam = true; pData->bShortScan = bShortScan; @@ -314,7 +324,7 @@ bool AstraFBP::run() // Call FDK_PreWeight to handle fan beam geometry. We treat // this as a cone beam setup of a single slice: - // TODO: TOffsets... + // TODO: TOffsets affects this preweighting... // We create a fake cudaPitchedPtr cudaPitchedPtr tmp; @@ -358,29 +368,7 @@ bool AstraFBP::run() } 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; + ok = FanBP_FBPWeighted(pData->D_volumeData, pData->volumePitch, pData->D_sinoData, pData->sinoPitch, pData->dims, pData->fanProjections); } else { ok = BP(pData->D_volumeData, pData->volumePitch, pData->D_sinoData, pData->sinoPitch, pData->dims, pData->angles, pData->TOffsets); diff --git a/cuda/2d/astra.h b/cuda/2d/astra.h index d2d7059..efbac1d 100644 --- a/cuda/2d/astra.h +++ b/cuda/2d/astra.h @@ -78,9 +78,10 @@ public: // fDetSize indicates the size of a detector pixel compared to a // volume pixel edge. // - // pfAngles will only be read from during this call. + // pfAngles, fanProjs will only be read from during this call. bool setFanGeometry(unsigned int iProjAngles, unsigned int iProjDets, + const astraCUDA::SFanProjection *fanProjs, const float *pfAngles, float fOriginSourceDistance, float fOriginDetectorDistance, diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp index 39740c8..63f4f40 100644 --- a/src/CudaFilteredBackProjectionAlgorithm.cpp +++ b/src/CudaFilteredBackProjectionAlgorithm.cpp @@ -32,6 +32,7 @@ $Id$ #include #include "astra/AstraObjectManager.h" +#include "../cuda/2d/astra.h" #include @@ -241,28 +242,65 @@ void CCudaFilteredBackProjectionAlgorithm::run(int _iNrIterations /* = 0 */) bool ok = true; - // TODO: off-center geometry, non-square pixels + // TODO: non-square pixels? ok &= m_pFBP->setReconstructionGeometry(volgeom.getGridColCount(), volgeom.getGridRowCount(), volgeom.getPixelLengthX()); - // TODO: off-center geometry int iDetectorCount; if (parprojgeom) { + + float *offsets, *angles, detSize, outputScale; + + ok = convertAstraGeometry(&volgeom, parprojgeom, offsets, angles, detSize, outputScale); + + ok &= m_pFBP->setProjectionGeometry(parprojgeom->getProjectionAngleCount(), parprojgeom->getDetectorCount(), - parprojgeom->getProjectionAngles(), + angles, parprojgeom->getDetectorWidth()); iDetectorCount = parprojgeom->getDetectorCount(); + + // TODO: Are detSize and outputScale handled correctly? + + if (offsets) + ok &= m_pFBP->setTOffsets(offsets); + ASTRA_ASSERT(ok); + + delete[] offsets; + delete[] angles; + } else if (fanprojgeom) { + + astraCUDA::SFanProjection* projs; + float outputScale; + + // FIXME: Implement this, and clean up the interface to AstraFBP. + if (abs(volgeom.getWindowMinX() + volgeom.getWindowMaxX()) > 0.00001 * volgeom.getPixelLengthX()) { + // Off-center volume geometry isn't supported yet + ASTRA_ASSERT(false); + } + if (abs(volgeom.getWindowMinY() + volgeom.getWindowMaxY()) > 0.00001 * volgeom.getPixelLengthY()) { + // Off-center volume geometry isn't supported yet + ASTRA_ASSERT(false); + } + + ok = convertAstraGeometry(&volgeom, fanprojgeom, projs, outputScale); + + // CHECKME: outputScale? + ok &= m_pFBP->setFanGeometry(fanprojgeom->getProjectionAngleCount(), - fanprojgeom->getDetectorCount(), - fanprojgeom->getProjectionAngles(), - fanprojgeom->getOriginSourceDistance(), - fanprojgeom->getOriginDetectorDistance(), - fanprojgeom->getDetectorWidth(), - m_bShortScan); + fanprojgeom->getDetectorCount(), + projs, + fanprojgeom->getProjectionAngles(), + fanprojgeom->getOriginSourceDistance(), + fanprojgeom->getOriginDetectorDistance(), + + fanprojgeom->getDetectorWidth(), + m_bShortScan); iDetectorCount = fanprojgeom->getDetectorCount(); + + delete[] projs; } else { assert(false); } @@ -271,11 +309,6 @@ void CCudaFilteredBackProjectionAlgorithm::run(int _iNrIterations /* = 0 */) ASTRA_ASSERT(ok); - const float *pfTOffsets = m_pSinogram->getGeometry()->getExtraDetectorOffset(); - if (pfTOffsets) - ok &= m_pFBP->setTOffsets(pfTOffsets); - ASTRA_ASSERT(ok); - ok &= m_pFBP->init(m_iGPUIndex); ASTRA_ASSERT(ok); -- cgit v1.2.3