summaryrefslogtreecommitdiffstats
path: root/include
diff options
context:
space:
mode:
Diffstat (limited to 'include')
-rw-r--r--include/astra/FanFlatBeamStripKernelProjector2D.inl54
1 files changed, 46 insertions, 8 deletions
diff --git a/include/astra/FanFlatBeamStripKernelProjector2D.inl b/include/astra/FanFlatBeamStripKernelProjector2D.inl
index 732c302..889d0a3 100644
--- a/include/astra/FanFlatBeamStripKernelProjector2D.inl
+++ b/include/astra/FanFlatBeamStripKernelProjector2D.inl
@@ -109,6 +109,9 @@ void CFanFlatBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, i
// POLICY: RAY PRIOR
if (!p.rayPrior(iRayIndex)) continue;
+ float32 dist_srcDetPixSquared = projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance() +
+ (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * DW * DW;
+
float32 sin_theta_left, cos_theta_left;
float32 sin_theta_right, cos_theta_right;
@@ -138,10 +141,11 @@ void CFanFlatBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, i
float32 inv_cos_theta_left = 1.0f / cos_theta_left;
float32 inv_cos_theta_right = 1.0f / cos_theta_right;
- float32 updateX_left = sin_theta_left * inv_cos_theta_left;
- float32 updateX_right = sin_theta_right * inv_cos_theta_right;
+ float32 updateX_left = sin_theta_left * inv_cos_theta_left; // tan(theta_left)
+ float32 updateX_right = sin_theta_right * inv_cos_theta_right; // tan(theta_right)
// Precalculate kernel limits
+ // BUG: If updateX_left > 1 (which can happen if theta_left >= pi/4 > theta), then T_l > U_l, and the expressions for res are no longer correct
float32 S_l = -0.5f * updateX_left;
if (S_l > 0) {S_l = -S_l;}
float32 T_l = -S_l;
@@ -185,6 +189,13 @@ void CFanFlatBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, i
xL += updateX_left;
xR += updateX_right;
+ float32 diffSrcYSquared;
+ if (switch_t)
+ diffSrcYSquared = m_pVolumeGeometry->pixelRowToCenterY(row) + cos_theta * projgeom->getOriginSourceDistance();
+ else
+ diffSrcYSquared = m_pVolumeGeometry->pixelRowToCenterY(row) - cos_theta * projgeom->getOriginSourceDistance();
+ diffSrcYSquared = diffSrcYSquared * diffSrcYSquared;
+
// for each affected col
for (col = x1L; col <= x1R; ++col) {
@@ -199,17 +210,25 @@ void CFanFlatBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, i
else if (x2R > U_r) res = x2R - (x2R-U_r)*(x2R-U_r)*inv_4T_r;
else if (x2R >= T_r) res = x2R;
else if (x2R > S_r) res = (x2R-S_r)*(x2R-S_r) * inv_4T_r;
- else { x2L -= 1.0f; x2R -= 1.0f; continue; }
+ else { x2L -= 1.0f; x2R -= 1.0f; p.pixelPosterior(iVolumeIndex); continue; }
// left
if (x2L <= S_l) {}
else if (x2L < T_l) res -= (x2L-S_l)*(x2L-S_l) * inv_4T_l;
else if (x2L <= U_l) res -= x2L;
else if (x2L < V_l) res -= x2L - (x2L-U_l)*(x2L-U_l)*inv_4T_l;
- else { x2L -= 1.0f; x2R -= 1.0f; continue; }
+ else { x2L -= 1.0f; x2R -= 1.0f; p.pixelPosterior(iVolumeIndex); continue; }
+
+ float32 diffSrcX;
+ if (switch_t)
+ diffSrcX = m_pVolumeGeometry->pixelColToCenterX(col) - sin_theta * projgeom->getOriginSourceDistance();
+ else
+ diffSrcX = m_pVolumeGeometry->pixelColToCenterX(col) + sin_theta * projgeom->getOriginSourceDistance();
+
+ float32 scale = sqrt(dist_srcDetPixSquared / (diffSrcYSquared + diffSrcX * diffSrcX));
// POLICY: ADD
- p.addWeight(iRayIndex, iVolumeIndex, PW*PH * res);
+ p.addWeight(iRayIndex, iVolumeIndex, PW*PH * res * scale);
// POLICY: PIXEL POSTERIOR
p.pixelPosterior(iVolumeIndex);
@@ -238,6 +257,9 @@ void CFanFlatBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, i
// POLICY: RAY PRIOR
if (!p.rayPrior(iRayIndex)) continue;
+ float32 dist_srcDetPixSquared = projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance() +
+ (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * DW * DW;
+
// get theta_l = alpha_left + theta and theta_r = alpha_right + theta
float32 sin_theta_left, cos_theta_left;
float32 sin_theta_right, cos_theta_right;
@@ -270,6 +292,7 @@ void CFanFlatBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, i
float32 updateX_right = cos_theta_right * inv_sin_theta_right;
// Precalculate kernel limits
+ // BUG: If updateX_left > 1 (which can happen if theta_left < pi/4 <= theta), then T_l > U_l, and the expressions for res are no longer correct
float32 S_l = -0.5f * updateX_left;
if (S_l > 0) { S_l = -S_l; }
float32 T_l = -S_l;
@@ -313,6 +336,13 @@ void CFanFlatBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, i
xL += updateX_left;
xR += updateX_right;
+ float32 diffSrcXSquared;
+ if (switch_t)
+ diffSrcXSquared = m_pVolumeGeometry->pixelColToCenterX(col) - sin_theta * projgeom->getOriginSourceDistance();
+ else
+ diffSrcXSquared = m_pVolumeGeometry->pixelColToCenterX(col) + sin_theta * projgeom->getOriginSourceDistance();
+ diffSrcXSquared = diffSrcXSquared * diffSrcXSquared;
+
// for each affected row
for (row = x1L; row <= x1R; ++row) {
@@ -328,17 +358,25 @@ void CFanFlatBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, i
else if (x2R > U_r) res = x2R - (x2R-U_r)*(x2R-U_r)*inv_4T_r;
else if (x2R >= T_r) res = x2R;
else if (x2R > S_r) res = (x2R-S_r)*(x2R-S_r) * inv_4T_r;
- else { x2L -= 1.0f; x2R -= 1.0f; continue; }
+ else { x2L -= 1.0f; x2R -= 1.0f; p.pixelPosterior(iVolumeIndex); continue; }
// left
if (x2L <= S_l) {}
else if (x2L < T_l) res -= (x2L-S_l)*(x2L-S_l) * inv_4T_l;
else if (x2L <= U_l) res -= x2L;
else if (x2L < V_l) res -= x2L - (x2L-U_l)*(x2L-U_l)*inv_4T_l;
- else { x2L -= 1.0f; x2R -= 1.0f; continue; }
+ else { x2L -= 1.0f; x2R -= 1.0f; p.pixelPosterior(iVolumeIndex); continue; }
+
+ float32 diffSrcY;
+ if (switch_t)
+ diffSrcY = m_pVolumeGeometry->pixelRowToCenterY(row) + cos_theta * projgeom->getOriginSourceDistance();
+ else
+ diffSrcY = m_pVolumeGeometry->pixelRowToCenterY(row) - cos_theta * projgeom->getOriginSourceDistance();
+
+ float32 scale = sqrt(dist_srcDetPixSquared / (diffSrcXSquared + diffSrcY * diffSrcY));
// POLICY: ADD
- p.addWeight(iRayIndex, iVolumeIndex, PW*PH * res);
+ p.addWeight(iRayIndex, iVolumeIndex, PW*PH * res * scale);
// POLICY: PIXEL POSTERIOR
p.pixelPosterior(iVolumeIndex);