summaryrefslogtreecommitdiffstats
path: root/cuda
diff options
context:
space:
mode:
Diffstat (limited to 'cuda')
-rw-r--r--cuda/2d/fbp.cu16
-rw-r--r--cuda/3d/cone_bp.cu21
-rw-r--r--cuda/3d/fdk.cu26
-rw-r--r--cuda/3d/fdk.h3
4 files changed, 50 insertions, 16 deletions
diff --git a/cuda/2d/fbp.cu b/cuda/2d/fbp.cu
index e7ed277..ecaf544 100644
--- a/cuda/2d/fbp.cu
+++ b/cuda/2d/fbp.cu
@@ -260,6 +260,7 @@ bool FBP::iterate(unsigned int iterations)
bool ok = false;
+ float fFanDetSize = 0.0f;
if (fanProjs) {
// Call FDK_PreWeight to handle fan beam geometry. We treat
// this as a cone beam setup of a single slice:
@@ -271,12 +272,12 @@ bool FBP::iterate(unsigned int iterations)
float *pfAngles = new float[dims.iProjAngles];
- float fOriginSource, fOriginDetector, fDetSize, fOffset;
+ float fOriginSource, fOriginDetector, fOffset;
for (unsigned int i = 0; i < dims.iProjAngles; ++i) {
bool ok = astra::getFanParameters(fanProjs[i], dims.iProjDets,
pfAngles[i],
fOriginSource, fOriginDetector,
- fDetSize, fOffset);
+ fFanDetSize, fOffset);
if (!ok) {
ASTRA_ERROR("FBP_CUDA: Failed to extract circular fan beam parameters from fan beam geometry");
return false;
@@ -300,7 +301,7 @@ bool FBP::iterate(unsigned int iterations)
astraCUDA3d::FDK_PreWeight(tmp, fOriginSource,
fOriginDetector, 0.0f,
- fDetSize, 1.0f,
+ fFanDetSize, 1.0f, /* fPixelSize */ 1.0f,
m_bShortScan, dims3d, pfAngles);
} else {
// TODO: How should different detector pixel size in different
@@ -326,12 +327,17 @@ bool FBP::iterate(unsigned int iterations)
}
- float fOutputScale = (M_PI / 2.0f) / (float)dims.iProjAngles;
-
if (fanProjs) {
+ float fOutputScale = 1.0 / (/*fPixelSize * fPixelSize * fPixelSize * */ fFanDetSize * fFanDetSize);
+
ok = FanBP_FBPWeighted(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, fanProjs, fOutputScale);
} else {
+ // scale by number of angles. For the fan-beam case, this is already
+ // handled by FDK_PreWeight
+ float fOutputScale = (M_PI / 2.0f) / (float)dims.iProjAngles;
+ //fOutputScale /= fDetSize * fDetSize;
+
ok = BP(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, parProjs, fOutputScale);
}
if(!ok)
diff --git a/cuda/3d/cone_bp.cu b/cuda/3d/cone_bp.cu
index 3dd0c97..2d12d00 100644
--- a/cuda/3d/cone_bp.cu
+++ b/cuda/3d/cone_bp.cu
@@ -126,6 +126,9 @@ __global__ void dev_cone_BP(void* D_volData, unsigned int volPitch, int startAng
float fVNum = fCv.w + fX * fCv.x + fY * fCv.y + fZ * fCv.z;
float fDen = fCd.w + fX * fCd.x + fY * fCd.y + fZ * fCd.z;
+ // fCd.w = -|| u v s || (determinant of 3x3 matrix with cols u,v,s)
+ // fDen = || u v (x-s) ||
+
float fU,fV, fr;
for (int idx = 0; idx < ZSIZE; idx++)
@@ -134,9 +137,17 @@ __global__ void dev_cone_BP(void* D_volData, unsigned int volPitch, int startAng
fU = fUNum * fr;
fV = fVNum * fr;
float fVal = tex3D(gT_coneProjTexture, fU, fAngle, fV);
- if (FDKWEIGHT)
+ if (FDKWEIGHT) {
+ // The correct factor here is this one:
+ // Z[idx] += (fr*fCd.w)*(fr*fCd.w)*fVal;
+ // This is the square of the inverse magnification factor
+ // from fX,fY,fZ to the detector.
+
+ // Since we are assuming we have a circular cone
+ // beam trajectory, fCd.w is constant, and we instead
+ // multiply by fCd.w*fCd.w in the FDK preweighting step.
Z[idx] += fr*fr*fVal;
- else
+ } else
Z[idx] += fVal;
fUNum += fCu.z;
@@ -255,7 +266,11 @@ bool ConeBP_Array(cudaPitchedPtr D_volumeData,
{
bindProjDataTexture(D_projArray);
- float fOutputScale = params.fOutputScale * params.fVolScaleX * params.fVolScaleY * params.fVolScaleZ;
+ float fOutputScale;
+ if (params.bFDKWeighting)
+ fOutputScale = params.fOutputScale / (params.fVolScaleX * params.fVolScaleY * params.fVolScaleZ);
+ else
+ fOutputScale = params.fOutputScale * (params.fVolScaleX * params.fVolScaleY * params.fVolScaleZ);
for (unsigned int th = 0; th < dims.iProjAngles; th += g_MaxAngles) {
unsigned int angleCount = g_MaxAngles;
diff --git a/cuda/3d/fdk.cu b/cuda/3d/fdk.cu
index 4e926f2..48194c4 100644
--- a/cuda/3d/fdk.cu
+++ b/cuda/3d/fdk.cu
@@ -59,7 +59,7 @@ __constant__ float gC_angle[g_MaxAngles];
// per-detector u/v shifts?
-__global__ void devFDK_preweight(void* D_projData, unsigned int projPitch, unsigned int startAngle, unsigned int endAngle, float fSrcOrigin, float fDetOrigin, float fZShift, float fDetUSize, float fDetVSize, const SDimensions3D dims)
+__global__ void devFDK_preweight(void* D_projData, unsigned int projPitch, unsigned int startAngle, unsigned int endAngle, float fSrcOrigin, float fDetOrigin, float fZShift, float fDetUSize, float fDetVSize, float fVoxSize, const SDimensions3D dims)
{
float* projData = (float*)D_projData;
int angle = startAngle + blockIdx.y * g_anglesPerWeightBlock + threadIdx.y;
@@ -83,10 +83,18 @@ __global__ void devFDK_preweight(void* D_projData, unsigned int projPitch, unsig
float fV = (startDetectorV - 0.5f*dims.iProjV + 0.5f) * fDetVSize + fZShift;
- //const float fW = fCentralRayLength;
- //const float fW = fCentralRayLength * (M_PI / 2.0f) / (float)dims.iProjAngles;
+ // Contributions to the weighting factors:
+ // fCentralRayLength / fRayLength : the main FDK preweighting factor
+ // fSrcOrigin / (fDetUSize * fCentralRayLength)
+ // : to adjust the filter to the det width
+ // || u v s || ^ 2 : see cone_bp.cu, FDKWEIGHT
+ // pi / (2 * iProjAngles) : scaling of the integral over angles
+ // fVoxSize ^ 2 : ...
+
const float fW1 = fSrcOrigin * fDetUSize * fDetVSize;
- const float fW = fCentralRayLength * fW1 * fW1 * (M_PI / 2.0f) / (float)dims.iProjAngles;
+ const float fW2 = fCentralRayLength / (fDetUSize * fSrcOrigin);
+ const float fW3 = fVoxSize * fVoxSize;
+ const float fW = fCentralRayLength * fW1 * fW1 * fW2 * fW3 * (M_PI / 2.0f) / (float)dims.iProjAngles;
for (int detectorV = startDetectorV; detectorV < endDetectorV; ++detectorV)
{
@@ -142,6 +150,8 @@ __global__ void devFDK_ParkerWeight(void* D_projData, unsigned int projPitch, un
fWeight = 0.0f;
}
+ fWeight *= 2; // adjust to effectively halved angular range
+
for (int detectorV = startDetectorV; detectorV < endDetectorV; ++detectorV)
{
@@ -156,7 +166,8 @@ __global__ void devFDK_ParkerWeight(void* D_projData, unsigned int projPitch, un
bool FDK_PreWeight(cudaPitchedPtr D_projData,
float fSrcOrigin, float fDetOrigin,
float fZShift,
- float fDetUSize, float fDetVSize, bool bShortScan,
+ float fDetUSize, float fDetVSize, float fVoxSize,
+ bool bShortScan,
const SDimensions3D& dims, const float* angles)
{
// The pre-weighting factor for a ray is the cosine of the angle between
@@ -168,7 +179,7 @@ bool FDK_PreWeight(cudaPitchedPtr D_projData,
int projPitch = D_projData.pitch/sizeof(float);
- devFDK_preweight<<<dimGrid, dimBlock>>>(D_projData.ptr, projPitch, 0, dims.iProjAngles, fSrcOrigin, fDetOrigin, fZShift, fDetUSize, fDetVSize, dims);
+ devFDK_preweight<<<dimGrid, dimBlock>>>(D_projData.ptr, projPitch, 0, dims.iProjAngles, fSrcOrigin, fDetOrigin, fZShift, fDetUSize, fDetVSize, fVoxSize, dims);
cudaTextForceKernelsCompletion();
@@ -310,8 +321,9 @@ bool FDK(cudaPitchedPtr D_volumeData,
#if 1
+ // NB: assuming cube voxels (params.fVolScaleX)
ok = FDK_PreWeight(D_projData, fSrcOrigin, fDetOrigin,
- fZShift, fDetUSize, fDetVSize,
+ fZShift, fDetUSize, fDetVSize, params.fVolScaleX,
bShortScan, dims, pfAngles);
#else
ok = true;
diff --git a/cuda/3d/fdk.h b/cuda/3d/fdk.h
index f4049e6..6f6e73b 100644
--- a/cuda/3d/fdk.h
+++ b/cuda/3d/fdk.h
@@ -35,7 +35,8 @@ namespace astraCUDA3d {
bool FDK_PreWeight(cudaPitchedPtr D_projData,
float fSrcOrigin, float fDetOrigin,
float fZShift,
- float fDetUSize, float fDetVSize, bool bShortScan,
+ float fDetUSize, float fDetVSize, float fVoxSize,
+ bool bShortScan,
const SDimensions3D& dims, const float* angles);
bool FDK(cudaPitchedPtr D_volumeData,
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-01 06:35:39 +0000