summaryrefslogtreecommitdiffstats
path: root/cuda/3d/astra3d.cu
diff options
context:
space:
mode:
Diffstat (limited to 'cuda/3d/astra3d.cu')
-rw-r--r--cuda/3d/astra3d.cu213
1 files changed, 60 insertions, 153 deletions
diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu
index 5670873..91f4530 100644
--- a/cuda/3d/astra3d.cu
+++ b/cuda/3d/astra3d.cu
@@ -67,34 +67,41 @@ template<typename ProjectionT>
static bool convertAstraGeometry_internal(const CVolumeGeometry3D* pVolGeom,
unsigned int iProjectionAngleCount,
ProjectionT*& pProjs,
- float& fOutputScale)
+ SProjectorParams3D& params)
{
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
- fOutputScale *= pVolGeom->getPixelLengthX();
+ //params.fOutputScale *= pVolGeom->getPixelLengthX();
return true;
}
@@ -108,10 +115,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 +129,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 +146,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 +167,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 +193,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 +214,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 +224,7 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
const CProjectionGeometry3D* pProjGeom,
SPar3DProjection*& pParProjs,
SConeProjection*& pConeProjs,
- float& fOutputScale)
+ SProjectorParams3D& params)
{
const CConeProjectionGeometry3D* conegeom = dynamic_cast<const CConeProjectionGeometry3D*>(pProjGeom);
const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast<const CParallelProjectionGeometry3D*>(pProjGeom);
@@ -240,13 +237,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,20 +257,17 @@ 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;
float fRelaxation;
SConeProjection* projs;
SPar3DProjection* parprojs;
- float fOutputScale;
-
bool initialized;
bool setStartReconstruction;
@@ -305,12 +299,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->fRelaxation = 1.0f;
@@ -361,7 +352,7 @@ bool AstraSIRT3d::setGeometry(const CVolumeGeometry3D* pVolGeom,
ok = convertAstraGeometry(pVolGeom, pProjGeom,
pData->parprojs, pData->projs,
- pData->fOutputScale);
+ pData->params);
if (!ok)
return false;
@@ -386,8 +377,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;
}
@@ -447,9 +438,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)
@@ -652,19 +643,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;
@@ -696,12 +684,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;
@@ -750,7 +735,7 @@ bool AstraCGLS3d::setGeometry(const CVolumeGeometry3D* pVolGeom,
ok = convertAstraGeometry(pVolGeom, pProjGeom,
pData->parprojs, pData->projs,
- pData->fOutputScale);
+ pData->params);
if (!ok)
return false;
@@ -774,8 +759,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;
}
@@ -829,9 +814,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)
@@ -1039,23 +1024,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) {
@@ -1093,10 +1078,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);
@@ -1104,7 +1089,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);
@@ -1129,21 +1114,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);
@@ -1179,9 +1163,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);
@@ -1204,21 +1188,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);
@@ -1255,9 +1239,9 @@ bool astraCudaBP_SIRTWeighted(float* pfVolume,
processSino3D<opSet>(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<opInvert>(D_pixelWeight, dims);
if (!ok) {
@@ -1272,9 +1256,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<opMul>(D_volumeData, D_pixelWeight, dims);
@@ -1307,81 +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;
-
- 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;
-
- dims.iRaysPerVoxelDim = iVoxelSuperSampling;
-
- 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
- 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;
-
-}
-
-
-
-
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-09-21 12:49:59 +0000