From 559d3e599b7306e2de64f2a584d72bc5c98b692b Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 4 Feb 2016 13:56:06 +0100 Subject: Refactor CUDA projector params into struct --- cuda/3d/astra3d.cu | 128 ++++++++++++++++++++++------------------------------- 1 file changed, 54 insertions(+), 74 deletions(-) (limited to 'cuda/3d/astra3d.cu') diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu index 8328229..454530e 100644 --- a/cuda/3d/astra3d.cu +++ b/cuda/3d/astra3d.cu @@ -67,7 +67,7 @@ template static bool convertAstraGeometry_internal(const CVolumeGeometry3D* pVolGeom, unsigned int iProjectionAngleCount, ProjectionT*& pProjs, - float& fOutputScale) + SProjectorParams3D& params) { assert(pVolGeom); assert(pProjs); @@ -94,7 +94,7 @@ static bool convertAstraGeometry_internal(const CVolumeGeometry3D* pVolGeom, } // CHECKME: Check factor - fOutputScale *= pVolGeom->getPixelLengthX(); + params.fOutputScale *= pVolGeom->getPixelLengthX(); return true; } @@ -108,10 +108,8 @@ bool convertAstraGeometry_dims(const CVolumeGeometry3D* pVolGeom, 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; + dims.iProjU = pProjGeom->getDetectorColCount(); + dims.iProjV = pProjGeom->getDetectorRowCount(); if (dims.iVolX <= 0 || dims.iVolX <= 0 || dims.iVolX <= 0) return false; @@ -124,7 +122,7 @@ bool convertAstraGeometry_dims(const CVolumeGeometry3D* pVolGeom, bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, const CParallelProjectionGeometry3D* pProjGeom, - SPar3DProjection*& pProjs, float& fOutputScale) + SPar3DProjection*& pProjs, SProjectorParams3D& params) { assert(pVolGeom); assert(pProjGeom); @@ -141,16 +139,14 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, bool ok; - fOutputScale = 1.0f; - - ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, fOutputScale); + ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, params); return ok; } bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, const CParallelVecProjectionGeometry3D* pProjGeom, - SPar3DProjection*& pProjs, float& fOutputScale) + SPar3DProjection*& pProjs, SProjectorParams3D& params) { assert(pVolGeom); assert(pProjGeom); @@ -164,16 +160,14 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, bool ok; - fOutputScale = 1.0f; - - ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, fOutputScale); + ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, params); return ok; } bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, const CConeProjectionGeometry3D* pProjGeom, - SConeProjection*& pProjs, float& fOutputScale) + SConeProjection*& pProjs, SProjectorParams3D& params) { assert(pVolGeom); assert(pProjGeom); @@ -192,16 +186,14 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, bool ok; - fOutputScale = 1.0f; - - ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, fOutputScale); + ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, params); return ok; } bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, const CConeVecProjectionGeometry3D* pProjGeom, - SConeProjection*& pProjs, float& fOutputScale) + SConeProjection*& pProjs, SProjectorParams3D& params) { assert(pVolGeom); assert(pProjGeom); @@ -215,9 +207,7 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, bool ok; - fOutputScale = 1.0f; - - ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, fOutputScale); + ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, params); return ok; } @@ -227,7 +217,7 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, const CProjectionGeometry3D* pProjGeom, SPar3DProjection*& pParProjs, SConeProjection*& pConeProjs, - float& fOutputScale) + SProjectorParams3D& params) { const CConeProjectionGeometry3D* conegeom = dynamic_cast(pProjGeom); const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast(pProjGeom); @@ -240,13 +230,13 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, bool ok; if (conegeom) { - ok = convertAstraGeometry(pVolGeom, conegeom, pConeProjs, fOutputScale); + ok = convertAstraGeometry(pVolGeom, conegeom, pConeProjs, params); } else if (conevec3dgeom) { - ok = convertAstraGeometry(pVolGeom, conevec3dgeom, pConeProjs, fOutputScale); + ok = convertAstraGeometry(pVolGeom, conevec3dgeom, pConeProjs, params); } else if (par3dgeom) { - ok = convertAstraGeometry(pVolGeom, par3dgeom, pParProjs, fOutputScale); + ok = convertAstraGeometry(pVolGeom, par3dgeom, pParProjs, params); } else if (parvec3dgeom) { - ok = convertAstraGeometry(pVolGeom, parvec3dgeom, pParProjs, fOutputScale); + ok = convertAstraGeometry(pVolGeom, parvec3dgeom, pParProjs, params); } else { ok = false; } @@ -260,19 +250,16 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom, class AstraSIRT3d_internal { public: SDimensions3D dims; + SProjectorParams3D params; CUDAProjectionType3d projType; float* angles; float fOriginSourceDistance; float fOriginDetectorDistance; - float fSourceZ; - float fDetSize; SConeProjection* projs; SPar3DProjection* parprojs; - float fOutputScale; - bool initialized; bool setStartReconstruction; @@ -304,12 +291,9 @@ AstraSIRT3d::AstraSIRT3d() pData->dims.iProjAngles = 0; pData->dims.iProjU = 0; pData->dims.iProjV = 0; - pData->dims.iRaysPerDetDim = 1; - pData->dims.iRaysPerVoxelDim = 1; pData->projs = 0; pData->parprojs = 0; - pData->fOutputScale = 1.0f; pData->initialized = false; pData->setStartReconstruction = false; @@ -358,7 +342,7 @@ bool AstraSIRT3d::setGeometry(const CVolumeGeometry3D* pVolGeom, ok = convertAstraGeometry(pVolGeom, pProjGeom, pData->parprojs, pData->projs, - pData->fOutputScale); + pData->params); if (!ok) return false; @@ -383,8 +367,8 @@ bool AstraSIRT3d::enableSuperSampling(unsigned int iVoxelSuperSampling, if (iVoxelSuperSampling == 0 || iDetectorSuperSampling == 0) return false; - pData->dims.iRaysPerVoxelDim = iVoxelSuperSampling; - pData->dims.iRaysPerDetDim = iDetectorSuperSampling; + pData->params.iRaysPerVoxelDim = iVoxelSuperSampling; + pData->params.iRaysPerDetDim = iDetectorSuperSampling; return true; } @@ -436,9 +420,9 @@ bool AstraSIRT3d::init() bool ok; if (pData->projType == PROJ_PARALLEL) { - ok = pData->sirt.setPar3DGeometry(pData->dims, pData->parprojs, pData->fOutputScale); + ok = pData->sirt.setPar3DGeometry(pData->dims, pData->parprojs, pData->params); } else { - ok = pData->sirt.setConeGeometry(pData->dims, pData->projs, pData->fOutputScale); + ok = pData->sirt.setConeGeometry(pData->dims, pData->projs, pData->params); } if (!ok) @@ -639,19 +623,16 @@ float AstraSIRT3d::computeDiffNorm() class AstraCGLS3d_internal { public: SDimensions3D dims; + SProjectorParams3D params; CUDAProjectionType3d projType; float* angles; float fOriginSourceDistance; float fOriginDetectorDistance; - float fSourceZ; - float fDetSize; SConeProjection* projs; SPar3DProjection* parprojs; - float fOutputScale; - bool initialized; bool setStartReconstruction; @@ -683,12 +664,9 @@ AstraCGLS3d::AstraCGLS3d() pData->dims.iProjAngles = 0; pData->dims.iProjU = 0; pData->dims.iProjV = 0; - pData->dims.iRaysPerDetDim = 1; - pData->dims.iRaysPerVoxelDim = 1; pData->projs = 0; pData->parprojs = 0; - pData->fOutputScale = 1.0f; pData->initialized = false; pData->setStartReconstruction = false; @@ -737,7 +715,7 @@ bool AstraCGLS3d::setGeometry(const CVolumeGeometry3D* pVolGeom, ok = convertAstraGeometry(pVolGeom, pProjGeom, pData->parprojs, pData->projs, - pData->fOutputScale); + pData->params); if (!ok) return false; @@ -761,8 +739,8 @@ bool AstraCGLS3d::enableSuperSampling(unsigned int iVoxelSuperSampling, if (iVoxelSuperSampling == 0 || iDetectorSuperSampling == 0) return false; - pData->dims.iRaysPerVoxelDim = iVoxelSuperSampling; - pData->dims.iRaysPerDetDim = iDetectorSuperSampling; + pData->params.iRaysPerVoxelDim = iVoxelSuperSampling; + pData->params.iRaysPerDetDim = iDetectorSuperSampling; return true; } @@ -816,9 +794,9 @@ bool AstraCGLS3d::init() bool ok; if (pData->projType == PROJ_PARALLEL) { - ok = pData->cgls.setPar3DGeometry(pData->dims, pData->parprojs, pData->fOutputScale); + ok = pData->cgls.setPar3DGeometry(pData->dims, pData->parprojs, pData->params); } else { - ok = pData->cgls.setConeGeometry(pData->dims, pData->projs, pData->fOutputScale); + ok = pData->cgls.setConeGeometry(pData->dims, pData->projs, pData->params); } if (!ok) @@ -1026,23 +1004,23 @@ bool astraCudaFP(const float* pfVolume, float* pfProjections, Cuda3DProjectionKernel projKernel) { SDimensions3D dims; + SProjectorParams3D params; + + params.iRaysPerDetDim = iDetectorSuperSampling; bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims); if (!ok) return false; - dims.iRaysPerDetDim = iDetectorSuperSampling; if (iDetectorSuperSampling == 0) return false; SPar3DProjection* pParProjs; SConeProjection* pConeProjs; - float outputScale; - ok = convertAstraGeometry(pVolGeom, pProjGeom, pParProjs, pConeProjs, - outputScale); + params); if (iGPUIndex != -1) { @@ -1080,10 +1058,10 @@ bool astraCudaFP(const float* pfVolume, float* pfProjections, if (pParProjs) { switch (projKernel) { case ker3d_default: - ok &= Par3DFP(D_volumeData, D_projData, dims, pParProjs, outputScale); + ok &= Par3DFP(D_volumeData, D_projData, dims, pParProjs, params); break; case ker3d_sum_square_weights: - ok &= Par3DFP_SumSqW(D_volumeData, D_projData, dims, pParProjs, outputScale*outputScale); + ok &= Par3DFP_SumSqW(D_volumeData, D_projData, dims, pParProjs, params); break; default: assert(false); @@ -1091,7 +1069,7 @@ bool astraCudaFP(const float* pfVolume, float* pfProjections, } else { switch (projKernel) { case ker3d_default: - ok &= ConeFP(D_volumeData, D_projData, dims, pConeProjs, outputScale); + ok &= ConeFP(D_volumeData, D_projData, dims, pConeProjs, params); break; default: assert(false); @@ -1116,21 +1094,20 @@ bool astraCudaBP(float* pfVolume, const float* pfProjections, int iGPUIndex, int iVoxelSuperSampling) { SDimensions3D dims; + SProjectorParams3D params; + + params.iRaysPerVoxelDim = iVoxelSuperSampling; bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims); if (!ok) return false; - dims.iRaysPerVoxelDim = iVoxelSuperSampling; - SPar3DProjection* pParProjs; SConeProjection* pConeProjs; - float outputScale; - ok = convertAstraGeometry(pVolGeom, pProjGeom, pParProjs, pConeProjs, - outputScale); + params); if (iGPUIndex != -1) { cudaSetDevice(iGPUIndex); @@ -1166,9 +1143,9 @@ bool astraCudaBP(float* pfVolume, const float* pfProjections, } if (pParProjs) - ok &= Par3DBP(D_volumeData, D_projData, dims, pParProjs, outputScale); + ok &= Par3DBP(D_volumeData, D_projData, dims, pParProjs, params); else - ok &= ConeBP(D_volumeData, D_projData, dims, pConeProjs, outputScale); + ok &= ConeBP(D_volumeData, D_projData, dims, pConeProjs, params); ok &= copyVolumeFromDevice(pfVolume, D_volumeData, dims, dims.iVolX); @@ -1191,21 +1168,21 @@ bool astraCudaBP_SIRTWeighted(float* pfVolume, int iGPUIndex, int iVoxelSuperSampling) { SDimensions3D dims; + SProjectorParams3D params; + + params.iRaysPerVoxelDim = iVoxelSuperSampling; bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims); if (!ok) return false; - dims.iRaysPerVoxelDim = iVoxelSuperSampling; SPar3DProjection* pParProjs; SConeProjection* pConeProjs; - float outputScale; - ok = convertAstraGeometry(pVolGeom, pProjGeom, pParProjs, pConeProjs, - outputScale); + params); if (iGPUIndex != -1) { cudaSetDevice(iGPUIndex); @@ -1242,9 +1219,9 @@ bool astraCudaBP_SIRTWeighted(float* pfVolume, processSino3D(D_projData, 1.0f, dims); if (pParProjs) - ok &= Par3DBP(D_pixelWeight, D_projData, dims, pParProjs, outputScale); + ok &= Par3DBP(D_pixelWeight, D_projData, dims, pParProjs, params); else - ok &= ConeBP(D_pixelWeight, D_projData, dims, pConeProjs, outputScale); + ok &= ConeBP(D_pixelWeight, D_projData, dims, pConeProjs, params); processVol3D(D_pixelWeight, dims); if (!ok) { @@ -1259,9 +1236,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, outputScale); + ok &= Par3DBP(D_volumeData, D_projData, dims, pParProjs, params); else - ok &= ConeBP(D_volumeData, D_projData, dims, pConeProjs, outputScale); + ok &= ConeBP(D_volumeData, D_projData, dims, pConeProjs, params); // Multiply with weights processVol3D(D_volumeData, D_pixelWeight, dims); @@ -1301,6 +1278,9 @@ bool astraCudaFDK(float* pfVolume, const float* pfProjections, int iGPUIndex, int iVoxelSuperSampling) { SDimensions3D dims; + SProjectorParams3D params; + + params.iRaysPerVoxelDim = iVoxelSuperSampling; bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims); @@ -1310,7 +1290,6 @@ bool astraCudaFDK(float* pfVolume, const float* pfProjections, if (!ok) return false; - dims.iRaysPerVoxelDim = iVoxelSuperSampling; if (iVoxelSuperSampling == 0) return false; @@ -1354,6 +1333,7 @@ bool astraCudaFDK(float* pfVolume, const float* pfProjections, // TODO: Offer interface for SrcZ, DetZ + // TODO: Voxel scaling ok &= FDK(D_volumeData, D_projData, fOriginSourceDistance, fOriginDetectorDistance, 0, 0, fDetUSize, fDetVSize, dims, pfAngles, bShortScan); -- cgit v1.2.3 From ae518f2954d8c6b9f1d5156595ccb1d7dc2ec581 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 5 Feb 2016 14:54:47 +0100 Subject: Process non-cubic-voxel astra geometries --- cuda/3d/astra3d.cu | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) (limited to 'cuda/3d/astra3d.cu') diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu index 454530e..42cece2 100644 --- a/cuda/3d/astra3d.cu +++ b/cuda/3d/astra3d.cu @@ -72,29 +72,36 @@ static bool convertAstraGeometry_internal(const CVolumeGeometry3D* pVolGeom, assert(pVolGeom); assert(pProjs); +#if 0 // 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; - +#endif // 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(); + float fx = 1.0f / pVolGeom->getPixelLengthX(); + float fy = 1.0f / pVolGeom->getPixelLengthY(); + float fz = 1.0f / pVolGeom->getPixelLengthZ(); for (int i = 0; i < iProjectionAngleCount; ++i) { // CHECKME: Order of scaling and translation pProjs[i].translate(dx, dy, dz); - pProjs[i].scale(factor); + pProjs[i].scale(fx, fy, fz); } + params.fVolScaleX = pVolGeom->getPixelLengthX(); + params.fVolScaleY = pVolGeom->getPixelLengthY(); + params.fVolScaleZ = pVolGeom->getPixelLengthZ(); + // CHECKME: Check factor - params.fOutputScale *= pVolGeom->getPixelLengthX(); + //params.fOutputScale *= pVolGeom->getPixelLengthX(); return true; } -- cgit v1.2.3 From 048755bab6b77c1da0050ed091e5007a60564adf Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Tue, 8 Mar 2016 15:41:38 +0100 Subject: Use CompositeGeometryManager for FDK Also fix a number of scaling/weighting issues in FDK, and switch to standard cone_bp with FDKWeighting for the BP step. --- cuda/3d/astra3d.cu | 80 ------------------------------------------------------ 1 file changed, 80 deletions(-) (limited to 'cuda/3d/astra3d.cu') diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu index dd6d551..91f4530 100644 --- a/cuda/3d/astra3d.cu +++ b/cuda/3d/astra3d.cu @@ -1291,84 +1291,4 @@ bool astraCudaBP_SIRTWeighted(float* pfVolume, -bool astraCudaFDK(float* pfVolume, const float* pfProjections, - const CVolumeGeometry3D* pVolGeom, - const CConeProjectionGeometry3D* pProjGeom, - bool bShortScan, - int iGPUIndex, int iVoxelSuperSampling) -{ - SDimensions3D dims; - SProjectorParams3D params; - - params.iRaysPerVoxelDim = iVoxelSuperSampling; - - bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims); - - // TODO: Check that pVolGeom is normalized, since we don't support - // other volume geometries yet - - if (!ok) - return false; - - - if (iVoxelSuperSampling == 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); - 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 &= copyProjectionsToDevice(pfProjections, D_projData, dims, dims.iProjU); - - ok &= zeroVolumeData(D_volumeData, dims); - - if (!ok) { - cudaFree(D_volumeData.ptr); - cudaFree(D_projData.ptr); - 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 - // TODO: Voxel scaling - ok &= FDK(D_volumeData, D_projData, fOriginSourceDistance, - fOriginDetectorDistance, 0, 0, fDetUSize, fDetVSize, - dims, pfAngles, bShortScan); - - ok &= copyVolumeFromDevice(pfVolume, D_volumeData, dims, dims.iVolX); - - - cudaFree(D_volumeData.ptr); - cudaFree(D_projData.ptr); - - return ok; - -} - - - - } -- cgit v1.2.3