summaryrefslogtreecommitdiffstats
path: root/include/astra
diff options
context:
space:
mode:
authorWillem Jan Palenstijn <wjp@usecode.org>2019-03-29 19:08:11 +0100
committerWillem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>2019-09-25 14:10:08 +0200
commit9c869d834ddc681299df9999787eb92bd6010dac (patch)
treec21e3efd97ca910388db92062b44dc9dfdd76955 /include/astra
parent15f53527e2f38e8eeb8bece79376b5e4686f3dd4 (diff)
downloadastra-9c869d834ddc681299df9999787eb92bd6010dac.tar.gz
astra-9c869d834ddc681299df9999787eb92bd6010dac.tar.bz2
astra-9c869d834ddc681299df9999787eb92bd6010dac.tar.xz
astra-9c869d834ddc681299df9999787eb92bd6010dac.zip
Adjust strip kernels to line integral scaling
Diffstat (limited to 'include/astra')
-rw-r--r--include/astra/FanFlatBeamStripKernelProjector2D.inl12
-rw-r--r--include/astra/ParallelBeamStripKernelProjector2D.inl8
2 files changed, 14 insertions, 6 deletions
diff --git a/include/astra/FanFlatBeamStripKernelProjector2D.inl b/include/astra/FanFlatBeamStripKernelProjector2D.inl
index 889d0a3..f5a688c 100644
--- a/include/astra/FanFlatBeamStripKernelProjector2D.inl
+++ b/include/astra/FanFlatBeamStripKernelProjector2D.inl
@@ -109,8 +109,12 @@ 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 dist_srcDetPixSquared = projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance() + (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * DW * DW;
+ dist_srcDetPixSquared = dist_srcDetPixSquared * dist_srcDetPixSquared / (projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance() * DW * DW);
+
+
+
+ //float32 InvRayWidthSquared = (projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance()) / dist_srcDetPixSquared;
float32 sin_theta_left, cos_theta_left;
float32 sin_theta_right, cos_theta_right;
@@ -257,8 +261,8 @@ 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 dist_srcDetPixSquared = projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance() + (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * DW * DW;
+ dist_srcDetPixSquared = dist_srcDetPixSquared * dist_srcDetPixSquared / (projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance() * DW * DW);
// get theta_l = alpha_left + theta and theta_r = alpha_right + theta
float32 sin_theta_left, cos_theta_left;
diff --git a/include/astra/ParallelBeamStripKernelProjector2D.inl b/include/astra/ParallelBeamStripKernelProjector2D.inl
index 0d775b3..bed02ab 100644
--- a/include/astra/ParallelBeamStripKernelProjector2D.inl
+++ b/include/astra/ParallelBeamStripKernelProjector2D.inl
@@ -142,6 +142,10 @@ void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom,
const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle];
+ const float32 rayWidth = fabs(proj->fDetUX * proj->fRayY - proj->fDetUY * proj->fRayX) /
+ sqrt(proj->fRayX * proj->fRayX + proj->fRayY * proj->fRayY);
+ const float32 relPixelArea = pixelArea / rayWidth;
+
bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY);
if (vertical) {
RxOverRy = proj->fRayX/proj->fRayY;
@@ -219,7 +223,7 @@ void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom,
else if (-S < offsetL) res -= 0.5f + offsetL;
else if (-T < offsetL) res -= 0.5f*(offsetL+T)*(offsetL+T)*invTminS;
- p.addWeight(iRayIndex, iVolumeIndex, pixelArea*res);
+ p.addWeight(iRayIndex, iVolumeIndex, relPixelArea*res);
p.pixelPosterior(iVolumeIndex);
}
}
@@ -272,7 +276,7 @@ void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom,
else if (-S < offsetL) res -= 0.5f + offsetL;
else if (-T < offsetL) res -= 0.5f*(offsetL+T)*(offsetL+T)*invTminS;
- p.addWeight(iRayIndex, iVolumeIndex, pixelArea*res);
+ p.addWeight(iRayIndex, iVolumeIndex, relPixelArea*res);
p.pixelPosterior(iVolumeIndex);
}
}