From 248639b4fee8659a4106dcc44d721149a1885018 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 5 Mar 2015 17:01:47 +0100 Subject: Add 3d geometry normalization functions --- cuda/3d/astra3d.cu | 150 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 150 insertions(+) (limited to 'cuda/3d/astra3d.cu') diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu index 0b9c70b..f672d6c 100644 --- a/cuda/3d/astra3d.cu +++ b/cuda/3d/astra3d.cu @@ -40,6 +40,12 @@ $Id$ #include "arith3d.h" #include "astra3d.h" +#include "astra/ParallelProjectionGeometry3D.h" +#include "astra/ParallelVecProjectionGeometry3D.h" +#include "astra/ConeProjectionGeometry3D.h" +#include "astra/ConeVecProjectionGeometry3D.h" +#include "astra/VolumeGeometry3D.h" + #include using namespace astraCUDA3d; @@ -137,6 +143,150 @@ static SPar3DProjection* genPar3DProjections(unsigned int iProjAngles, + +// adjust pProjs to normalize volume geometry +template +static bool convertAstraGeometry_internal(const CVolumeGeometry3D* pVolGeom, + unsigned int iProjectionAngleCount, + ProjectionT*& pProjs, + float& fOutputScale) +{ + assert(pVolGeom); + assert(pProjs); + + // TODO: Relative instead of absolute + const float EPS = 0.00001f; + if (abs(pVolGeom->getPixelLengthX() - pVolGeom->getPixelLengthY()) > EPS) + return false; + if (abs(pVolGeom->getPixelLengthX() - pVolGeom->getPixelLengthZ()) > EPS) + return false; + + + // Translate + float dx = -(pVolGeom->getWindowMinX() + pVolGeom->getWindowMaxX()) / 2; + float dy = -(pVolGeom->getWindowMinY() + pVolGeom->getWindowMaxY()) / 2; + float dz = -(pVolGeom->getWindowMinZ() + pVolGeom->getWindowMaxZ()) / 2; + + float factor = 1.0f / pVolGeom->getPixelLengthX(); + + for (int i = 0; i < iProjectionAngleCount; ++i) { + // CHECKME: Order of scaling and translation + pProjs[i].translate(dx, dy, dz); + pProjs[i].scale(factor); + } + + // CHECKME: Check factor + fOutputScale *= pVolGeom->getPixelLengthX(); + + return true; +} + + + +bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, + const CParallelProjectionGeometry3D* pProjGeom, + SPar3DProjection*& pProjs, float& fOutputScale) +{ + assert(pVolGeom); + assert(pProjGeom); + assert(pProjGeom->getProjectionAngles()); + + int nth = pProjGeom->getProjectionCount(); + + pProjs = genPar3DProjections(nth, + pProjGeom->getDetectorColCount(), + pProjGeom->getDetectorRowCount(), + pProjGeom->getDetectorSpacingX(), + pProjGeom->getDetectorSpacingY(), + pProjGeom->getProjectionAngles()); + + bool ok; + + fOutputScale = 1.0f; + + ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, fOutputScale); + + return ok; +} + +bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, + const CParallelVecProjectionGeometry3D* pProjGeom, + SPar3DProjection*& pProjs, float& fOutputScale) +{ + assert(pVolGeom); + assert(pProjGeom); + assert(pProjGeom->getProjectionVectors()); + + int nth = pProjGeom->getProjectionCount(); + + pProjs = new SPar3DProjection[nth]; + for (int i = 0; i < nth; ++i) + pProjs[i] = pProjGeom->getProjectionVectors()[i]; + + bool ok; + + fOutputScale = 1.0f; + + ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, fOutputScale); + + return ok; +} + +bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, + const CConeProjectionGeometry3D* pProjGeom, + SConeProjection*& pProjs, float& fOutputScale) +{ + assert(pVolGeom); + assert(pProjGeom); + assert(pProjGeom->getProjectionAngles()); + + int nth = pProjGeom->getProjectionCount(); + + pProjs = genConeProjections(nth, + pProjGeom->getDetectorColCount(), + pProjGeom->getDetectorRowCount(), + pProjGeom->getOriginSourceDistance(), + pProjGeom->getOriginDetectorDistance(), + pProjGeom->getDetectorSpacingX(), + pProjGeom->getDetectorSpacingY(), + pProjGeom->getProjectionAngles()); + + bool ok; + + fOutputScale = 1.0f; + + ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, fOutputScale); + + return ok; +} + +bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, + const CConeVecProjectionGeometry3D* pProjGeom, + SConeProjection*& pProjs, float& fOutputScale) +{ + assert(pVolGeom); + assert(pProjGeom); + assert(pProjGeom->getProjectionAngles()); + + int nth = pProjGeom->getProjectionCount(); + + pProjs = new SConeProjection[nth]; + for (int i = 0; i < nth; ++i) + pProjs[i] = pProjGeom->getProjectionVectors()[i]; + + bool ok; + + fOutputScale = 1.0f; + + ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, fOutputScale); + + return ok; +} + + + + + class AstraSIRT3d_internal { public: SDimensions3D dims; -- cgit v1.2.3 From 5304d08cd1ab7b8d778c367912934376eb92370f Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Mon, 9 Mar 2015 15:43:56 +0100 Subject: Allow non-centered volume geometry in SIRT3D and CGLS3D --- cuda/3d/astra3d.cu | 284 ++++++++++++++--------------------------------------- 1 file changed, 76 insertions(+), 208 deletions(-) (limited to 'cuda/3d/astra3d.cu') diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu index f672d6c..426f3a0 100644 --- a/cuda/3d/astra3d.cu +++ b/cuda/3d/astra3d.cu @@ -182,6 +182,20 @@ static bool convertAstraGeometry_internal(const CVolumeGeometry3D* pVolGeom, } +void convertAstraGeometry_dims(const CVolumeGeometry3D* pVolGeom, + const CProjectionGeometry3D* pProjGeom, + SDimensions3D& dims) +{ + dims.iVolX = pVolGeom->getGridColCount(); + dims.iVolY = pVolGeom->getGridRowCount(); + dims.iVolZ = pVolGeom->getGridSliceCount(); + dims.iProjAngles = pProjGeom->getProjectionCount(); + dims.iProjU = pProjGeom->getDetectorColCount(), + dims.iProjV = pProjGeom->getDetectorRowCount(), + dims.iRaysPerDetDim = 1; + dims.iRaysPerVoxelDim = 1; +} + bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, const CParallelProjectionGeometry3D* pProjGeom, @@ -370,127 +384,55 @@ AstraSIRT3d::~AstraSIRT3d() pData = 0; } -bool AstraSIRT3d::setReconstructionGeometry(unsigned int iVolX, - unsigned int iVolY, - unsigned int iVolZ/*, - float fPixelSize = 1.0f*/) +bool AstraSIRT3d::setGeometry(const CVolumeGeometry3D* pVolGeom, + const CProjectionGeometry3D* pProjGeom) { if (pData->initialized) return false; - pData->dims.iVolX = iVolX; - pData->dims.iVolY = iVolY; - pData->dims.iVolZ = iVolZ; + convertAstraGeometry_dims(pVolGeom, pProjGeom, pData->dims); - return (iVolX > 0 && iVolY > 0 && iVolZ > 0); -} - - -bool AstraSIRT3d::setPar3DGeometry(unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - const SPar3DProjection* projs) -{ - if (pData->initialized) + if (pData->dims.iVolX <= 0 || pData->dims.iVolX <= 0 || pData->dims.iVolX <= 0) return false; - - pData->dims.iProjAngles = iProjAngles; - pData->dims.iProjU = iProjU; - pData->dims.iProjV = iProjV; - - if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || projs == 0) + if (pData->dims.iProjAngles <= 0 || pData->dims.iProjU <= 0 || pData->dims.iProjV <= 0) return false; - pData->parprojs = new SPar3DProjection[iProjAngles]; - memcpy(pData->parprojs, projs, iProjAngles * sizeof(projs[0])); - - pData->projType = PROJ_PARALLEL; - - return true; -} - -bool AstraSIRT3d::setPar3DGeometry(unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - float fDetUSize, - float fDetVSize, - const float *pfAngles) -{ - if (pData->initialized) - return false; - - if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0) - return false; - - SPar3DProjection* p = genPar3DProjections(iProjAngles, - iProjU, iProjV, - fDetUSize, fDetVSize, - pfAngles); - pData->dims.iProjAngles = iProjAngles; - pData->dims.iProjU = iProjU; - pData->dims.iProjV = iProjV; - - pData->parprojs = p; - pData->projType = PROJ_PARALLEL; - - return true; -} - - + const CConeProjectionGeometry3D* conegeom = dynamic_cast(pProjGeom); + const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast(pProjGeom); + const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast(pProjGeom); + const CConeVecProjectionGeometry3D* conevec3dgeom = dynamic_cast(pProjGeom); -bool AstraSIRT3d::setConeGeometry(unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - const SConeProjection* projs) -{ - if (pData->initialized) - return false; + float outputScale; + bool ok; - pData->dims.iProjAngles = iProjAngles; - pData->dims.iProjU = iProjU; - pData->dims.iProjV = iProjV; + pData->projs = 0; + pData->parprojs = 0; + + if (conegeom) { + ok = convertAstraGeometry(pVolGeom, conegeom, pData->projs, outputScale); + pData->projType = PROJ_PARALLEL; + } else if (conevec3dgeom) { + ok = convertAstraGeometry(pVolGeom, conevec3dgeom, pData->projs, outputScale); + pData->projType = PROJ_PARALLEL; + } else if (par3dgeom) { + ok = convertAstraGeometry(pVolGeom, par3dgeom, pData->parprojs, outputScale); + pData->projType = PROJ_CONE; + } else if (parvec3dgeom) { + ok = convertAstraGeometry(pVolGeom, parvec3dgeom, pData->parprojs, outputScale); + pData->projType = PROJ_CONE; + } else { + ok = false; + } - if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || projs == 0) + if (!ok) return false; - pData->projs = new SConeProjection[iProjAngles]; - memcpy(pData->projs, projs, iProjAngles * sizeof(projs[0])); - pData->projType = PROJ_CONE; + // TODO: Handle outputScale return true; } -bool AstraSIRT3d::setConeGeometry(unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - float fOriginSourceDistance, - float fOriginDetectorDistance, - float fDetUSize, - float fDetVSize, - const float *pfAngles) -{ - if (pData->initialized) - return false; - - if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0) - return false; - - SConeProjection* p = genConeProjections(iProjAngles, - iProjU, iProjV, - fOriginSourceDistance, - fOriginDetectorDistance, - fDetUSize, fDetVSize, - pfAngles); - pData->dims.iProjAngles = iProjAngles; - pData->dims.iProjU = iProjU; - pData->dims.iProjV = iProjV; - - pData->projs = p; - pData->projType = PROJ_CONE; - - return true; -} bool AstraSIRT3d::enableSuperSampling(unsigned int iVoxelSuperSampling, unsigned int iDetectorSuperSampling) @@ -837,125 +779,51 @@ AstraCGLS3d::~AstraCGLS3d() pData = 0; } -bool AstraCGLS3d::setReconstructionGeometry(unsigned int iVolX, - unsigned int iVolY, - unsigned int iVolZ/*, - float fPixelSize = 1.0f*/) +bool AstraCGLS3d::setGeometry(const CVolumeGeometry3D* pVolGeom, + const CProjectionGeometry3D* pProjGeom) { if (pData->initialized) return false; - pData->dims.iVolX = iVolX; - pData->dims.iVolY = iVolY; - pData->dims.iVolZ = iVolZ; + convertAstraGeometry_dims(pVolGeom, pProjGeom, pData->dims); - return (iVolX > 0 && iVolY > 0 && iVolZ > 0); -} - - -bool AstraCGLS3d::setPar3DGeometry(unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - const SPar3DProjection* projs) -{ - if (pData->initialized) + if (pData->dims.iVolX <= 0 || pData->dims.iVolX <= 0 || pData->dims.iVolX <= 0) return false; - - pData->dims.iProjAngles = iProjAngles; - pData->dims.iProjU = iProjU; - pData->dims.iProjV = iProjV; - - if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || projs == 0) + if (pData->dims.iProjAngles <= 0 || pData->dims.iProjU <= 0 || pData->dims.iProjV <= 0) return false; - pData->parprojs = new SPar3DProjection[iProjAngles]; - memcpy(pData->parprojs, projs, iProjAngles * sizeof(projs[0])); + const CConeProjectionGeometry3D* conegeom = dynamic_cast(pProjGeom); + const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast(pProjGeom); + const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast(pProjGeom); + const CConeVecProjectionGeometry3D* conevec3dgeom = dynamic_cast(pProjGeom); - pData->projType = PROJ_PARALLEL; - - return true; -} - -bool AstraCGLS3d::setPar3DGeometry(unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - float fDetUSize, - float fDetVSize, - const float *pfAngles) -{ - if (pData->initialized) - return false; - - if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0) - return false; - - SPar3DProjection* p = genPar3DProjections(iProjAngles, - iProjU, iProjV, - fDetUSize, fDetVSize, - pfAngles); - pData->dims.iProjAngles = iProjAngles; - pData->dims.iProjU = iProjU; - pData->dims.iProjV = iProjV; - - pData->parprojs = p; - pData->projType = PROJ_PARALLEL; - - return true; -} - - - -bool AstraCGLS3d::setConeGeometry(unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - const SConeProjection* projs) -{ - if (pData->initialized) - return false; - - pData->dims.iProjAngles = iProjAngles; - pData->dims.iProjU = iProjU; - pData->dims.iProjV = iProjV; - - if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || projs == 0) - return false; - - pData->projs = new SConeProjection[iProjAngles]; - memcpy(pData->projs, projs, iProjAngles * sizeof(projs[0])); - - pData->projType = PROJ_CONE; - - return true; -} + float outputScale; + bool ok; -bool AstraCGLS3d::setConeGeometry(unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - float fOriginSourceDistance, - float fOriginDetectorDistance, - float fDetUSize, - float fDetVSize, - const float *pfAngles) -{ - if (pData->initialized) - return false; + pData->projs = 0; + pData->parprojs = 0; + + if (conegeom) { + ok = convertAstraGeometry(pVolGeom, conegeom, pData->projs, outputScale); + pData->projType = PROJ_PARALLEL; + } else if (conevec3dgeom) { + ok = convertAstraGeometry(pVolGeom, conevec3dgeom, pData->projs, outputScale); + pData->projType = PROJ_PARALLEL; + } else if (par3dgeom) { + ok = convertAstraGeometry(pVolGeom, par3dgeom, pData->parprojs, outputScale); + pData->projType = PROJ_CONE; + } else if (parvec3dgeom) { + ok = convertAstraGeometry(pVolGeom, parvec3dgeom, pData->parprojs, outputScale); + pData->projType = PROJ_CONE; + } else { + ok = false; + } - if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0) + if (!ok) return false; - SConeProjection* p = genConeProjections(iProjAngles, - iProjU, iProjV, - fOriginSourceDistance, - fOriginDetectorDistance, - fDetUSize, fDetVSize, - pfAngles); - - pData->dims.iProjAngles = iProjAngles; - pData->dims.iProjU = iProjU; - pData->dims.iProjV = iProjV; - pData->projs = p; - pData->projType = PROJ_CONE; + // TODO: Handle outputScale return true; } -- cgit v1.2.3 From 9eb68c39c62a8e674e3dbe50252528226c6593ff Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 11 Mar 2015 12:06:11 +0100 Subject: Adjust interface slightly --- cuda/3d/astra3d.cu | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) (limited to 'cuda/3d/astra3d.cu') diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu index 426f3a0..5b1f363 100644 --- a/cuda/3d/astra3d.cu +++ b/cuda/3d/astra3d.cu @@ -182,7 +182,7 @@ static bool convertAstraGeometry_internal(const CVolumeGeometry3D* pVolGeom, } -void convertAstraGeometry_dims(const CVolumeGeometry3D* pVolGeom, +bool convertAstraGeometry_dims(const CVolumeGeometry3D* pVolGeom, const CProjectionGeometry3D* pProjGeom, SDimensions3D& dims) { @@ -194,6 +194,13 @@ void convertAstraGeometry_dims(const CVolumeGeometry3D* pVolGeom, dims.iProjV = pProjGeom->getDetectorRowCount(), dims.iRaysPerDetDim = 1; dims.iRaysPerVoxelDim = 1; + + if (dims.iVolX <= 0 || dims.iVolX <= 0 || dims.iVolX <= 0) + return false; + if (dims.iProjAngles <= 0 || dims.iProjU <= 0 || dims.iProjV <= 0) + return false; + + return true; } @@ -390,11 +397,9 @@ bool AstraSIRT3d::setGeometry(const CVolumeGeometry3D* pVolGeom, if (pData->initialized) return false; - convertAstraGeometry_dims(pVolGeom, pProjGeom, pData->dims); + bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, pData->dims); - if (pData->dims.iVolX <= 0 || pData->dims.iVolX <= 0 || pData->dims.iVolX <= 0) - return false; - if (pData->dims.iProjAngles <= 0 || pData->dims.iProjU <= 0 || pData->dims.iProjV <= 0) + if (!ok) return false; const CConeProjectionGeometry3D* conegeom = dynamic_cast(pProjGeom); @@ -403,7 +408,6 @@ bool AstraSIRT3d::setGeometry(const CVolumeGeometry3D* pVolGeom, const CConeVecProjectionGeometry3D* conevec3dgeom = dynamic_cast(pProjGeom); float outputScale; - bool ok; pData->projs = 0; pData->parprojs = 0; @@ -785,11 +789,9 @@ bool AstraCGLS3d::setGeometry(const CVolumeGeometry3D* pVolGeom, if (pData->initialized) return false; - convertAstraGeometry_dims(pVolGeom, pProjGeom, pData->dims); + bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, pData->dims); - if (pData->dims.iVolX <= 0 || pData->dims.iVolX <= 0 || pData->dims.iVolX <= 0) - return false; - if (pData->dims.iProjAngles <= 0 || pData->dims.iProjU <= 0 || pData->dims.iProjV <= 0) + if (!ok) return false; const CConeProjectionGeometry3D* conegeom = dynamic_cast(pProjGeom); @@ -798,7 +800,6 @@ bool AstraCGLS3d::setGeometry(const CVolumeGeometry3D* pVolGeom, const CConeVecProjectionGeometry3D* conevec3dgeom = dynamic_cast(pProjGeom); float outputScale; - bool ok; pData->projs = 0; pData->parprojs = 0; -- cgit v1.2.3 From 140f64028a6c06895ba7dad8997e14b7a05aadab Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 11 Mar 2015 12:07:48 +0100 Subject: Let astraCudaFDK use utility functions --- cuda/3d/astra3d.cu | 36 ++++++++++++++---------------------- 1 file changed, 14 insertions(+), 22 deletions(-) (limited to 'cuda/3d/astra3d.cu') diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu index 5b1f363..0e94fb8 100644 --- a/cuda/3d/astra3d.cu +++ b/cuda/3d/astra3d.cu @@ -1679,33 +1679,19 @@ bool astraCudaPar3DBP_SIRTWeighted(float* pfVolume, bool astraCudaFDK(float* pfVolume, const float* pfProjections, - unsigned int iVolX, - unsigned int iVolY, - unsigned int iVolZ, - unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - float fOriginSourceDistance, - float fOriginDetectorDistance, - float fDetUSize, - float fDetVSize, - const float *pfAngles, + const CVolumeGeometry3D* pVolGeom, + const CConeProjectionGeometry3D* pProjGeom, bool bShortScan, int iGPUIndex, int iVoxelSuperSampling) { SDimensions3D dims; - dims.iVolX = iVolX; - dims.iVolY = iVolY; - dims.iVolZ = iVolZ; - if (iVolX == 0 || iVolY == 0 || iVolZ == 0) - return false; + bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims); - dims.iProjAngles = iProjAngles; - dims.iProjU = iProjU; - dims.iProjV = iProjV; + // TODO: Check that pVolGeom is normalized, since we don't support + // other volume geometries yet - if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0) + if (!ok) return false; dims.iRaysPerVoxelDim = iVoxelSuperSampling; @@ -1722,9 +1708,8 @@ bool astraCudaFDK(float* pfVolume, const float* pfProjections, return false; } - cudaPitchedPtr D_volumeData = allocateVolumeData(dims); - bool ok = D_volumeData.ptr; + ok = D_volumeData.ptr; if (!ok) return false; @@ -1745,6 +1730,13 @@ bool astraCudaFDK(float* pfVolume, const float* pfProjections, return false; } + float fOriginSourceDistance = pProjGeom->getOriginSourceDistance(); + float fOriginDetectorDistance = pProjGeom->getOriginDetectorDistance(); + float fDetUSize = pProjGeom->getDetectorSpacingX(); + float fDetVSize = pProjGeom->getDetectorSpacingY(); + const float *pfAngles = pProjGeom->getProjectionAngles(); + + // TODO: Offer interface for SrcZ, DetZ ok &= FDK(D_volumeData, D_projData, fOriginSourceDistance, fOriginDetectorDistance, 0, 0, fDetUSize, fDetVSize, -- cgit v1.2.3 From e188bcdaaffee075adf5fa4371453d91bcb71225 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 11 Mar 2015 14:43:49 +0100 Subject: Add another utility function --- cuda/3d/astra3d.cu | 99 +++++++++++++++++++++++++++++------------------------- 1 file changed, 53 insertions(+), 46 deletions(-) (limited to 'cuda/3d/astra3d.cu') diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu index 0e94fb8..eff928d 100644 --- a/cuda/3d/astra3d.cu +++ b/cuda/3d/astra3d.cu @@ -305,6 +305,37 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, } +bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, + const CConeVecProjectionGeometry3D* pProjGeom, + SPar3DProjection*& pParProjs, + SConeProjection*& pConeProjs, + float& fOutputScale) +{ + const CConeProjectionGeometry3D* conegeom = dynamic_cast(pProjGeom); + const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast(pProjGeom); + const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast(pProjGeom); + const CConeVecProjectionGeometry3D* conevec3dgeom = dynamic_cast(pProjGeom); + + pConeProjs = 0; + pParProjs = 0; + + bool ok; + + if (conegeom) { + ok = convertAstraGeometry(pVolGeom, conegeom, pConeProjs, outputScale); + } else if (conevec3dgeom) { + ok = convertAstraGeometry(pVolGeom, conevec3dgeom, pConeProjs, outputScale); + } else if (par3dgeom) { + ok = convertAstraGeometry(pVolGeom, par3dgeom, pParProjs, outputScale); + } else if (parvec3dgeom) { + ok = convertAstraGeometry(pVolGeom, parvec3dgeom, pParProjs, outputScale); + } else { + ok = false; + } + + return ok; +} + @@ -402,35 +433,23 @@ bool AstraSIRT3d::setGeometry(const CVolumeGeometry3D* pVolGeom, if (!ok) return false; - const CConeProjectionGeometry3D* conegeom = dynamic_cast(pProjGeom); - const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast(pProjGeom); - const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast(pProjGeom); - const CConeVecProjectionGeometry3D* conevec3dgeom = dynamic_cast(pProjGeom); - - float outputScale; - pData->projs = 0; pData->parprojs = 0; + float outputScale; - if (conegeom) { - ok = convertAstraGeometry(pVolGeom, conegeom, pData->projs, outputScale); - pData->projType = PROJ_PARALLEL; - } else if (conevec3dgeom) { - ok = convertAstraGeometry(pVolGeom, conevec3dgeom, pData->projs, outputScale); - pData->projType = PROJ_PARALLEL; - } else if (par3dgeom) { - ok = convertAstraGeometry(pVolGeom, par3dgeom, pData->parprojs, outputScale); - pData->projType = PROJ_CONE; - } else if (parvec3dgeom) { - ok = convertAstraGeometry(pVolGeom, parvec3dgeom, pData->parprojs, outputScale); - pData->projType = PROJ_CONE; - } else { - ok = false; - } - + ok = convertAstraGeometry(pVolGeom, pProjGeom, + pData->parprojs, pData->projs, + outputScale); if (!ok) return false; + if (pData->projs) { + assert(pData->parprojs == 0); + pData->projType = PROJ_CONE; + } else { + assert(pData->parprojs != 0); + pData->projType = PROJ_PARALLEL; + } // TODO: Handle outputScale @@ -794,35 +813,23 @@ bool AstraCGLS3d::setGeometry(const CVolumeGeometry3D* pVolGeom, if (!ok) return false; - const CConeProjectionGeometry3D* conegeom = dynamic_cast(pProjGeom); - const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast(pProjGeom); - const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast(pProjGeom); - const CConeVecProjectionGeometry3D* conevec3dgeom = dynamic_cast(pProjGeom); - - float outputScale; - pData->projs = 0; pData->parprojs = 0; + float outputScale; - if (conegeom) { - ok = convertAstraGeometry(pVolGeom, conegeom, pData->projs, outputScale); - pData->projType = PROJ_PARALLEL; - } else if (conevec3dgeom) { - ok = convertAstraGeometry(pVolGeom, conevec3dgeom, pData->projs, outputScale); - pData->projType = PROJ_PARALLEL; - } else if (par3dgeom) { - ok = convertAstraGeometry(pVolGeom, par3dgeom, pData->parprojs, outputScale); - pData->projType = PROJ_CONE; - } else if (parvec3dgeom) { - ok = convertAstraGeometry(pVolGeom, parvec3dgeom, pData->parprojs, outputScale); - pData->projType = PROJ_CONE; - } else { - ok = false; - } - + ok = convertAstraGeometry(pVolGeom, pProjGeom, + pData->parprojs, pData->projs, + outputScale); if (!ok) return false; + if (pData->projs) { + assert(pData->parprojs == 0); + pData->projType = PROJ_CONE; + } else { + assert(pData->parprojs != 0); + pData->projType = PROJ_PARALLEL; + } // TODO: Handle outputScale -- cgit v1.2.3 From 18d12242207d1113c3015b451f522531168e626a Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 11 Mar 2015 17:27:44 +0100 Subject: Add flexible volgeom3d support to astraCudaBP_SIRTWeighted --- cuda/3d/astra3d.cu | 95 ++++++++++++++++++++---------------------------------- 1 file changed, 35 insertions(+), 60 deletions(-) (limited to 'cuda/3d/astra3d.cu') diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu index eff928d..2f7ea99 100644 --- a/cuda/3d/astra3d.cu +++ b/cuda/3d/astra3d.cu @@ -306,7 +306,7 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, - const CConeVecProjectionGeometry3D* pProjGeom, + const CProjectionGeometry3D* pProjGeom, SPar3DProjection*& pParProjs, SConeProjection*& pConeProjs, float& fOutputScale) @@ -322,13 +322,13 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, bool ok; if (conegeom) { - ok = convertAstraGeometry(pVolGeom, conegeom, pConeProjs, outputScale); + ok = convertAstraGeometry(pVolGeom, conegeom, pConeProjs, fOutputScale); } else if (conevec3dgeom) { - ok = convertAstraGeometry(pVolGeom, conevec3dgeom, pConeProjs, outputScale); + ok = convertAstraGeometry(pVolGeom, conevec3dgeom, pConeProjs, fOutputScale); } else if (par3dgeom) { - ok = convertAstraGeometry(pVolGeom, par3dgeom, pParProjs, outputScale); + ok = convertAstraGeometry(pVolGeom, par3dgeom, pParProjs, fOutputScale); } else if (parvec3dgeom) { - ok = convertAstraGeometry(pVolGeom, parvec3dgeom, pParProjs, outputScale); + ok = convertAstraGeometry(pVolGeom, parvec3dgeom, pParProjs, fOutputScale); } else { ok = false; } @@ -1471,40 +1471,6 @@ bool astraCudaPar3DBP(float* pfVolume, const float* pfProjections, return ok; } -// This computes the column weights, divides by them, and adds the -// result to the current volume. This is both more expensive and more -// GPU memory intensive than the regular BP, but allows saving system RAM. -bool astraCudaPar3DBP_SIRTWeighted(float* pfVolume, const float* pfProjections, - unsigned int iVolX, - unsigned int iVolY, - unsigned int iVolZ, - unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - float fDetUSize, - float fDetVSize, - const float *pfAngles, - int iGPUIndex, int iVoxelSuperSampling) -{ - if (iVolX == 0 || iVolY == 0 || iVolZ == 0) - return false; - if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0) - return false; - - SPar3DProjection* p = genPar3DProjections(iProjAngles, - iProjU, iProjV, - fDetUSize, fDetVSize, - pfAngles); - - bool ok; - ok = astraCudaPar3DBP_SIRTWeighted(pfVolume, pfProjections, iVolX, iVolY, iVolZ, - iProjAngles, iProjU, iProjV, p, iGPUIndex, iVoxelSuperSampling); - - delete[] p; - - return ok; -} - bool astraCudaPar3DBP(float* pfVolume, const float* pfProjections, unsigned int iVolX, @@ -1582,33 +1548,30 @@ bool astraCudaPar3DBP(float* pfVolume, const float* pfProjections, // This computes the column weights, divides by them, and adds the // result to the current volume. This is both more expensive and more // GPU memory intensive than the regular BP, but allows saving system RAM. -bool astraCudaPar3DBP_SIRTWeighted(float* pfVolume, +bool astraCudaBP_SIRTWeighted(float* pfVolume, const float* pfProjections, - unsigned int iVolX, - unsigned int iVolY, - unsigned int iVolZ, - unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - const SPar3DProjection *pfAngles, + const CVolumeGeometry3D* pVolGeom, + const CProjectionGeometry3D* pProjGeom, int iGPUIndex, int iVoxelSuperSampling) { SDimensions3D dims; - dims.iVolX = iVolX; - dims.iVolY = iVolY; - dims.iVolZ = iVolZ; - if (iVolX == 0 || iVolY == 0 || iVolZ == 0) + bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims); + if (!ok) return false; - dims.iProjAngles = iProjAngles; - dims.iProjU = iProjU; - dims.iProjV = iProjV; + dims.iRaysPerVoxelDim = iVoxelSuperSampling; - if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0) - return false; + SPar3DProjection* pParProjs; + SConeProjection* pConeProjs; - dims.iRaysPerVoxelDim = iVoxelSuperSampling; + float outputScale; + + ok = convertAstraGeometry(pVolGeom, pProjGeom, + pParProjs, pConeProjs, + outputScale); + + // TODO: OutputScale if (iGPUIndex != -1) { cudaSetDevice(iGPUIndex); @@ -1621,7 +1584,7 @@ bool astraCudaPar3DBP_SIRTWeighted(float* pfVolume, cudaPitchedPtr D_pixelWeight = allocateVolumeData(dims); - bool ok = D_pixelWeight.ptr; + ok = D_pixelWeight.ptr; if (!ok) return false; @@ -1643,7 +1606,12 @@ bool astraCudaPar3DBP_SIRTWeighted(float* pfVolume, // Compute weights ok &= zeroVolumeData(D_pixelWeight, dims); processSino3D(D_projData, 1.0f, dims); - ok &= Par3DBP(D_pixelWeight, D_projData, dims, pfAngles); + + if (pParProjs) + ok &= Par3DBP(D_pixelWeight, D_projData, dims, pParProjs); + else + ok &= ConeBP(D_pixelWeight, D_projData, dims, pConeProjs); + processVol3D(D_pixelWeight, dims); if (!ok) { cudaFree(D_pixelWeight.ptr); @@ -1656,7 +1624,11 @@ bool astraCudaPar3DBP_SIRTWeighted(float* pfVolume, dims, dims.iProjU); ok &= zeroVolumeData(D_volumeData, dims); // Do BP into D_volumeData - ok &= Par3DBP(D_volumeData, D_projData, dims, pfAngles); + if (pParProjs) + ok &= Par3DBP(D_volumeData, D_projData, dims, pParProjs); + else + ok &= ConeBP(D_volumeData, D_projData, dims, pConeProjs); + // Multiply with weights processVol3D(D_volumeData, D_pixelWeight, dims); @@ -1679,6 +1651,9 @@ bool astraCudaPar3DBP_SIRTWeighted(float* pfVolume, cudaFree(D_volumeData.ptr); cudaFree(D_projData.ptr); + delete[] pParProjs; + delete[] pConeProjs; + return ok; } -- cgit v1.2.3 From 6909836555afe155ffc3897ef2189ed0562bb045 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 11 Mar 2015 18:44:53 +0100 Subject: Add flexible volgeom3d support to astraCudaBP --- cuda/3d/astra3d.cu | 176 ++++++----------------------------------------------- 1 file changed, 18 insertions(+), 158 deletions(-) (limited to 'cuda/3d/astra3d.cu') diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu index 2f7ea99..97bebf4 100644 --- a/cuda/3d/astra3d.cu +++ b/cuda/3d/astra3d.cu @@ -1331,173 +1331,30 @@ bool astraCudaPar3DFP(const float* pfVolume, float* pfProjections, } -bool astraCudaConeBP(float* pfVolume, const float* pfProjections, - unsigned int iVolX, - unsigned int iVolY, - unsigned int iVolZ, - unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - float fOriginSourceDistance, - float fOriginDetectorDistance, - float fDetUSize, - float fDetVSize, - const float *pfAngles, - int iGPUIndex, int iVoxelSuperSampling) -{ - if (iVolX == 0 || iVolY == 0 || iVolZ == 0) - return false; - if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0) - return false; - - SConeProjection* p = genConeProjections(iProjAngles, - iProjU, iProjV, - fOriginSourceDistance, - fOriginDetectorDistance, - fDetUSize, fDetVSize, - pfAngles); - - bool ok; - ok = astraCudaConeBP(pfVolume, pfProjections, iVolX, iVolY, iVolZ, - iProjAngles, iProjU, iProjV, p, iGPUIndex, iVoxelSuperSampling); - - delete[] p; - - return ok; -} -bool astraCudaConeBP(float* pfVolume, const float* pfProjections, - unsigned int iVolX, - unsigned int iVolY, - unsigned int iVolZ, - unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - const SConeProjection *pfAngles, - int iGPUIndex, int iVoxelSuperSampling) +bool astraCudaBP(float* pfVolume, const float* pfProjections, + const CVolumeGeometry3D* pVolGeom, + const CProjectionGeometry3D* pProjGeom, + int iGPUIndex, int iVoxelSuperSampling) { SDimensions3D dims; - dims.iVolX = iVolX; - dims.iVolY = iVolY; - dims.iVolZ = iVolZ; - if (iVolX == 0 || iVolY == 0 || iVolZ == 0) - return false; - - dims.iProjAngles = iProjAngles; - dims.iProjU = iProjU; - dims.iProjV = iProjV; - - if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0) - return false; - - dims.iRaysPerVoxelDim = iVoxelSuperSampling; - - if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0) - return false; - - if (iGPUIndex != -1) { - cudaSetDevice(iGPUIndex); - cudaError_t err = cudaGetLastError(); - - // Ignore errors caused by calling cudaSetDevice multiple times - if (err != cudaSuccess && err != cudaErrorSetOnActiveProcess) - return false; - } - - cudaPitchedPtr D_volumeData = allocateVolumeData(dims); - bool ok = D_volumeData.ptr; + bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims); if (!ok) return false; - cudaPitchedPtr D_projData = allocateProjectionData(dims); - ok = D_projData.ptr; - if (!ok) { - cudaFree(D_volumeData.ptr); - return false; - } - - ok &= copyProjectionsToDevice(pfProjections, D_projData, - dims, dims.iProjU); - - ok &= zeroVolumeData(D_volumeData, dims); - - if (!ok) { - cudaFree(D_volumeData.ptr); - cudaFree(D_projData.ptr); - return false; - } - - ok &= ConeBP(D_volumeData, D_projData, dims, pfAngles); - - ok &= copyVolumeFromDevice(pfVolume, D_volumeData, dims, dims.iVolX); - - - cudaFree(D_volumeData.ptr); - cudaFree(D_projData.ptr); - - return ok; - -} - -bool astraCudaPar3DBP(float* pfVolume, const float* pfProjections, - unsigned int iVolX, - unsigned int iVolY, - unsigned int iVolZ, - unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - float fDetUSize, - float fDetVSize, - const float *pfAngles, - int iGPUIndex, int iVoxelSuperSampling) -{ - if (iVolX == 0 || iVolY == 0 || iVolZ == 0) - return false; - if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0) - return false; - - SPar3DProjection* p = genPar3DProjections(iProjAngles, - iProjU, iProjV, - fDetUSize, fDetVSize, - pfAngles); - - bool ok; - ok = astraCudaPar3DBP(pfVolume, pfProjections, iVolX, iVolY, iVolZ, - iProjAngles, iProjU, iProjV, p, iGPUIndex, iVoxelSuperSampling); - - delete[] p; - - return ok; -} - - -bool astraCudaPar3DBP(float* pfVolume, const float* pfProjections, - unsigned int iVolX, - unsigned int iVolY, - unsigned int iVolZ, - unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - const SPar3DProjection *pfAngles, - int iGPUIndex, int iVoxelSuperSampling) -{ - SDimensions3D dims; + dims.iRaysPerVoxelDim = iVoxelSuperSampling; - dims.iVolX = iVolX; - dims.iVolY = iVolY; - dims.iVolZ = iVolZ; - if (iVolX == 0 || iVolY == 0 || iVolZ == 0) - return false; + SPar3DProjection* pParProjs; + SConeProjection* pConeProjs; - dims.iProjAngles = iProjAngles; - dims.iProjU = iProjU; - dims.iProjV = iProjV; + float outputScale; - if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0) - return false; + ok = convertAstraGeometry(pVolGeom, pProjGeom, + pParProjs, pConeProjs, + outputScale); - dims.iRaysPerVoxelDim = iVoxelSuperSampling; + // TODO: OutputScale if (iGPUIndex != -1) { cudaSetDevice(iGPUIndex); @@ -1510,7 +1367,7 @@ bool astraCudaPar3DBP(float* pfVolume, const float* pfProjections, cudaPitchedPtr D_volumeData = allocateVolumeData(dims); - bool ok = D_volumeData.ptr; + ok = D_volumeData.ptr; if (!ok) return false; @@ -1532,7 +1389,10 @@ bool astraCudaPar3DBP(float* pfVolume, const float* pfProjections, return false; } - ok &= Par3DBP(D_volumeData, D_projData, dims, pfAngles); + if (pParProjs) + ok &= Par3DBP(D_volumeData, D_projData, dims, pParProjs); + else + ok &= ConeBP(D_volumeData, D_projData, dims, pConeProjs); ok &= copyVolumeFromDevice(pfVolume, D_volumeData, dims, dims.iVolX); -- cgit v1.2.3 From 57ee6b85884b8226b26b7415ef151b4a6e63337c Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 12 Mar 2015 11:53:40 +0100 Subject: Add flexible volgeom3d support to astraCudaFP --- cuda/3d/astra3d.cu | 204 +++++++++-------------------------------------------- 1 file changed, 33 insertions(+), 171 deletions(-) (limited to 'cuda/3d/astra3d.cu') diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu index 97bebf4..b2375f3 100644 --- a/cuda/3d/astra3d.cu +++ b/cuda/3d/astra3d.cu @@ -1103,179 +1103,31 @@ float AstraCGLS3d::computeDiffNorm() -bool astraCudaConeFP(const float* pfVolume, float* pfProjections, - unsigned int iVolX, - unsigned int iVolY, - unsigned int iVolZ, - unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - float fOriginSourceDistance, - float fOriginDetectorDistance, - float fDetUSize, - float fDetVSize, - const float *pfAngles, - int iGPUIndex, int iDetectorSuperSampling) -{ - if (iVolX == 0 || iVolY == 0 || iVolZ == 0) - return false; - if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0) - return false; - - SConeProjection* p = genConeProjections(iProjAngles, - iProjU, iProjV, - fOriginSourceDistance, - fOriginDetectorDistance, - fDetUSize, fDetVSize, - pfAngles); - - bool ok; - ok = astraCudaConeFP(pfVolume, pfProjections, iVolX, iVolY, iVolZ, - iProjAngles, iProjU, iProjV, p, iGPUIndex, iDetectorSuperSampling); - - delete[] p; - - return ok; -} - -bool astraCudaConeFP(const float* pfVolume, float* pfProjections, - unsigned int iVolX, - unsigned int iVolY, - unsigned int iVolZ, - unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - const SConeProjection *pfAngles, - int iGPUIndex, int iDetectorSuperSampling) +bool astraCudaFP(const float* pfVolume, float* pfProjections, + const CVolumeGeometry3D* pVolGeom, + const CProjectionGeometry3D* pProjGeom, + int iGPUIndex, int iDetectorSuperSampling, + Cuda3DProjectionKernel projKernel) { SDimensions3D dims; - dims.iVolX = iVolX; - dims.iVolY = iVolY; - dims.iVolZ = iVolZ; - if (iVolX == 0 || iVolY == 0 || iVolZ == 0) - return false; - - dims.iProjAngles = iProjAngles; - dims.iProjU = iProjU; - dims.iProjV = iProjV; - - if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0) + bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims); + if (!ok) return false; dims.iRaysPerDetDim = iDetectorSuperSampling; - if (iDetectorSuperSampling == 0) return false; - if (iGPUIndex != -1) { - cudaSetDevice(iGPUIndex); - cudaError_t err = cudaGetLastError(); - - // Ignore errors caused by calling cudaSetDevice multiple times - if (err != cudaSuccess && err != cudaErrorSetOnActiveProcess) - return false; - } - - cudaPitchedPtr D_volumeData = allocateVolumeData(dims); - bool ok = D_volumeData.ptr; - if (!ok) - return false; - - cudaPitchedPtr D_projData = allocateProjectionData(dims); - ok = D_projData.ptr; - if (!ok) { - cudaFree(D_volumeData.ptr); - return false; - } - - ok &= copyVolumeToDevice(pfVolume, D_volumeData, dims, dims.iVolX); - - ok &= zeroProjectionData(D_projData, dims); - - if (!ok) { - cudaFree(D_volumeData.ptr); - cudaFree(D_projData.ptr); - return false; - } - - ok &= ConeFP(D_volumeData, D_projData, dims, pfAngles, 1.0f); - - ok &= copyProjectionsFromDevice(pfProjections, D_projData, - dims, dims.iProjU); - - - cudaFree(D_volumeData.ptr); - cudaFree(D_projData.ptr); - - return ok; - -} - -bool astraCudaPar3DFP(const float* pfVolume, float* pfProjections, - unsigned int iVolX, - unsigned int iVolY, - unsigned int iVolZ, - unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - float fDetUSize, - float fDetVSize, - const float *pfAngles, - int iGPUIndex, int iDetectorSuperSampling, - Cuda3DProjectionKernel projKernel) -{ - if (iVolX == 0 || iVolY == 0 || iVolZ == 0) - return false; - if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0) - return false; - - SPar3DProjection* p = genPar3DProjections(iProjAngles, - iProjU, iProjV, - fDetUSize, fDetVSize, - pfAngles); - - bool ok; - ok = astraCudaPar3DFP(pfVolume, pfProjections, iVolX, iVolY, iVolZ, - iProjAngles, iProjU, iProjV, p, iGPUIndex, iDetectorSuperSampling, - projKernel); - - delete[] p; - - return ok; -} - - -bool astraCudaPar3DFP(const float* pfVolume, float* pfProjections, - unsigned int iVolX, - unsigned int iVolY, - unsigned int iVolZ, - unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - const SPar3DProjection *pfAngles, - int iGPUIndex, int iDetectorSuperSampling, - Cuda3DProjectionKernel projKernel) -{ - SDimensions3D dims; - - dims.iVolX = iVolX; - dims.iVolY = iVolY; - dims.iVolZ = iVolZ; - if (iVolX == 0 || iVolY == 0 || iVolZ == 0) - return false; - - dims.iProjAngles = iProjAngles; - dims.iProjU = iProjU; - dims.iProjV = iProjV; + SPar3DProjection* pParProjs; + SConeProjection* pConeProjs; - if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0) - return false; + float outputScale; - dims.iRaysPerDetDim = iDetectorSuperSampling; + ok = convertAstraGeometry(pVolGeom, pProjGeom, + pParProjs, pConeProjs, + outputScale); - if (iDetectorSuperSampling == 0) - return false; if (iGPUIndex != -1) { cudaSetDevice(iGPUIndex); @@ -1288,7 +1140,7 @@ bool astraCudaPar3DFP(const float* pfVolume, float* pfProjections, cudaPitchedPtr D_volumeData = allocateVolumeData(dims); - bool ok = D_volumeData.ptr; + ok = D_volumeData.ptr; if (!ok) return false; @@ -1309,15 +1161,25 @@ bool astraCudaPar3DFP(const float* pfVolume, float* pfProjections, return false; } - switch (projKernel) { - case ker3d_default: - ok &= Par3DFP(D_volumeData, D_projData, dims, pfAngles, 1.0f); - break; - case ker3d_sum_square_weights: - ok &= Par3DFP_SumSqW(D_volumeData, D_projData, dims, pfAngles, 1.0f); - break; - default: - assert(false); + if (pParProjs) { + switch (projKernel) { + case ker3d_default: + ok &= Par3DFP(D_volumeData, D_projData, dims, pParProjs, 1.0f); + break; + case ker3d_sum_square_weights: + ok &= Par3DFP_SumSqW(D_volumeData, D_projData, dims, pParProjs, 1.0f); + break; + default: + assert(false); + } + } else { + switch (projKernel) { + case ker3d_default: + ok &= ConeFP(D_volumeData, D_projData, dims, pConeProjs, 1.0f); + break; + default: + assert(false); + } } ok &= copyProjectionsFromDevice(pfProjections, D_projData, -- cgit v1.2.3 From a70ad8df8fc2a3da63fc91dd18bbfd55be7a89dd Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 12 Mar 2015 12:30:47 +0100 Subject: Add outputScale argument to 3D CUDA BP --- cuda/3d/astra3d.cu | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) (limited to 'cuda/3d/astra3d.cu') diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu index b2375f3..7589416 100644 --- a/cuda/3d/astra3d.cu +++ b/cuda/3d/astra3d.cu @@ -1252,9 +1252,9 @@ bool astraCudaBP(float* pfVolume, const float* pfProjections, } if (pParProjs) - ok &= Par3DBP(D_volumeData, D_projData, dims, pParProjs); + ok &= Par3DBP(D_volumeData, D_projData, dims, pParProjs, 1.0f); else - ok &= ConeBP(D_volumeData, D_projData, dims, pConeProjs); + ok &= ConeBP(D_volumeData, D_projData, dims, pConeProjs, 1.0f); ok &= copyVolumeFromDevice(pfVolume, D_volumeData, dims, dims.iVolX); @@ -1330,9 +1330,9 @@ bool astraCudaBP_SIRTWeighted(float* pfVolume, processSino3D(D_projData, 1.0f, dims); if (pParProjs) - ok &= Par3DBP(D_pixelWeight, D_projData, dims, pParProjs); + ok &= Par3DBP(D_pixelWeight, D_projData, dims, pParProjs, 1.0f); else - ok &= ConeBP(D_pixelWeight, D_projData, dims, pConeProjs); + ok &= ConeBP(D_pixelWeight, D_projData, dims, pConeProjs, 1.0f); processVol3D(D_pixelWeight, dims); if (!ok) { @@ -1347,9 +1347,9 @@ bool astraCudaBP_SIRTWeighted(float* pfVolume, ok &= zeroVolumeData(D_volumeData, dims); // Do BP into D_volumeData if (pParProjs) - ok &= Par3DBP(D_volumeData, D_projData, dims, pParProjs); + ok &= Par3DBP(D_volumeData, D_projData, dims, pParProjs, 1.0f); else - ok &= ConeBP(D_volumeData, D_projData, dims, pConeProjs); + ok &= ConeBP(D_volumeData, D_projData, dims, pConeProjs, 1.0f); // Multiply with weights processVol3D(D_volumeData, D_pixelWeight, dims); -- cgit v1.2.3 From 42db106d9b66639312d874e4f35e4e9ff7a407d0 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 12 Mar 2015 14:03:17 +0100 Subject: Scale CUDA 3D FP/BP output with volume pixel size --- cuda/3d/astra3d.cu | 48 +++++++++++++++++++++--------------------------- 1 file changed, 21 insertions(+), 27 deletions(-) (limited to 'cuda/3d/astra3d.cu') diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu index 7589416..ae79efb 100644 --- a/cuda/3d/astra3d.cu +++ b/cuda/3d/astra3d.cu @@ -353,7 +353,7 @@ public: SConeProjection* projs; SPar3DProjection* parprojs; - float fPixelSize; + float fOutputScale; bool initialized; bool setStartReconstruction; @@ -390,6 +390,8 @@ AstraSIRT3d::AstraSIRT3d() pData->dims.iRaysPerVoxelDim = 1; pData->projs = 0; + pData->parprojs = 0; + pData->fOutputScale = 1.0f; pData->initialized = false; pData->setStartReconstruction = false; @@ -435,11 +437,10 @@ bool AstraSIRT3d::setGeometry(const CVolumeGeometry3D* pVolGeom, pData->projs = 0; pData->parprojs = 0; - float outputScale; ok = convertAstraGeometry(pVolGeom, pProjGeom, pData->parprojs, pData->projs, - outputScale); + pData->fOutputScale); if (!ok) return false; @@ -451,8 +452,6 @@ bool AstraSIRT3d::setGeometry(const CVolumeGeometry3D* pVolGeom, pData->projType = PROJ_PARALLEL; } - // TODO: Handle outputScale - return true; } @@ -519,9 +518,9 @@ bool AstraSIRT3d::init() bool ok; if (pData->projType == PROJ_PARALLEL) { - ok = pData->sirt.setPar3DGeometry(pData->dims, pData->parprojs); + ok = pData->sirt.setPar3DGeometry(pData->dims, pData->parprojs, pData->fOutputScale); } else { - ok = pData->sirt.setConeGeometry(pData->dims, pData->projs); + ok = pData->sirt.setConeGeometry(pData->dims, pData->projs, pData->fOutputScale); } if (!ok) @@ -733,7 +732,7 @@ public: SConeProjection* projs; SPar3DProjection* parprojs; - float fPixelSize; + float fOutputScale; bool initialized; bool setStartReconstruction; @@ -770,6 +769,8 @@ AstraCGLS3d::AstraCGLS3d() pData->dims.iRaysPerVoxelDim = 1; pData->projs = 0; + pData->parprojs = 0; + pData->fOutputScale = 1.0f; pData->initialized = false; pData->setStartReconstruction = false; @@ -815,11 +816,10 @@ bool AstraCGLS3d::setGeometry(const CVolumeGeometry3D* pVolGeom, pData->projs = 0; pData->parprojs = 0; - float outputScale; ok = convertAstraGeometry(pVolGeom, pProjGeom, pData->parprojs, pData->projs, - outputScale); + pData->fOutputScale); if (!ok) return false; @@ -831,8 +831,6 @@ bool AstraCGLS3d::setGeometry(const CVolumeGeometry3D* pVolGeom, pData->projType = PROJ_PARALLEL; } - // TODO: Handle outputScale - return true; } @@ -900,9 +898,9 @@ bool AstraCGLS3d::init() bool ok; if (pData->projType == PROJ_PARALLEL) { - ok = pData->cgls.setPar3DGeometry(pData->dims, pData->parprojs); + ok = pData->cgls.setPar3DGeometry(pData->dims, pData->parprojs, pData->fOutputScale); } else { - ok = pData->cgls.setConeGeometry(pData->dims, pData->projs); + ok = pData->cgls.setConeGeometry(pData->dims, pData->projs, pData->fOutputScale); } if (!ok) @@ -1164,10 +1162,10 @@ bool astraCudaFP(const float* pfVolume, float* pfProjections, if (pParProjs) { switch (projKernel) { case ker3d_default: - ok &= Par3DFP(D_volumeData, D_projData, dims, pParProjs, 1.0f); + ok &= Par3DFP(D_volumeData, D_projData, dims, pParProjs, outputScale); break; case ker3d_sum_square_weights: - ok &= Par3DFP_SumSqW(D_volumeData, D_projData, dims, pParProjs, 1.0f); + ok &= Par3DFP_SumSqW(D_volumeData, D_projData, dims, pParProjs, outputScale*outputScale); break; default: assert(false); @@ -1175,7 +1173,7 @@ bool astraCudaFP(const float* pfVolume, float* pfProjections, } else { switch (projKernel) { case ker3d_default: - ok &= ConeFP(D_volumeData, D_projData, dims, pConeProjs, 1.0f); + ok &= ConeFP(D_volumeData, D_projData, dims, pConeProjs, outputScale); break; default: assert(false); @@ -1216,8 +1214,6 @@ bool astraCudaBP(float* pfVolume, const float* pfProjections, pParProjs, pConeProjs, outputScale); - // TODO: OutputScale - if (iGPUIndex != -1) { cudaSetDevice(iGPUIndex); cudaError_t err = cudaGetLastError(); @@ -1252,9 +1248,9 @@ bool astraCudaBP(float* pfVolume, const float* pfProjections, } if (pParProjs) - ok &= Par3DBP(D_volumeData, D_projData, dims, pParProjs, 1.0f); + ok &= Par3DBP(D_volumeData, D_projData, dims, pParProjs, outputScale); else - ok &= ConeBP(D_volumeData, D_projData, dims, pConeProjs, 1.0f); + ok &= ConeBP(D_volumeData, D_projData, dims, pConeProjs, outputScale); ok &= copyVolumeFromDevice(pfVolume, D_volumeData, dims, dims.iVolX); @@ -1293,8 +1289,6 @@ bool astraCudaBP_SIRTWeighted(float* pfVolume, pParProjs, pConeProjs, outputScale); - // TODO: OutputScale - if (iGPUIndex != -1) { cudaSetDevice(iGPUIndex); cudaError_t err = cudaGetLastError(); @@ -1330,9 +1324,9 @@ bool astraCudaBP_SIRTWeighted(float* pfVolume, processSino3D(D_projData, 1.0f, dims); if (pParProjs) - ok &= Par3DBP(D_pixelWeight, D_projData, dims, pParProjs, 1.0f); + ok &= Par3DBP(D_pixelWeight, D_projData, dims, pParProjs, outputScale); else - ok &= ConeBP(D_pixelWeight, D_projData, dims, pConeProjs, 1.0f); + ok &= ConeBP(D_pixelWeight, D_projData, dims, pConeProjs, outputScale); processVol3D(D_pixelWeight, dims); if (!ok) { @@ -1347,9 +1341,9 @@ bool astraCudaBP_SIRTWeighted(float* pfVolume, ok &= zeroVolumeData(D_volumeData, dims); // Do BP into D_volumeData if (pParProjs) - ok &= Par3DBP(D_volumeData, D_projData, dims, pParProjs, 1.0f); + ok &= Par3DBP(D_volumeData, D_projData, dims, pParProjs, outputScale); else - ok &= ConeBP(D_volumeData, D_projData, dims, pConeProjs, 1.0f); + ok &= ConeBP(D_volumeData, D_projData, dims, pConeProjs, outputScale); // Multiply with weights processVol3D(D_volumeData, D_pixelWeight, dims); -- cgit v1.2.3 From 7584ffbd6748bcca8c3f7ed2dc961be01f2fcfdc Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 16 Jul 2015 18:09:51 +0200 Subject: Fix assert --- cuda/3d/astra3d.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'cuda/3d/astra3d.cu') diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu index ae79efb..3815a1a 100644 --- a/cuda/3d/astra3d.cu +++ b/cuda/3d/astra3d.cu @@ -287,7 +287,7 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, { assert(pVolGeom); assert(pProjGeom); - assert(pProjGeom->getProjectionAngles()); + assert(pProjGeom->getProjectionVectors()); int nth = pProjGeom->getProjectionCount(); -- cgit v1.2.3