From 9c869d834ddc681299df9999787eb92bd6010dac Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 29 Mar 2019 19:08:11 +0100 Subject: Adjust strip kernels to line integral scaling --- include/astra/ParallelBeamStripKernelProjector2D.inl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) (limited to 'include/astra/ParallelBeamStripKernelProjector2D.inl') 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); } } -- cgit v1.2.3 From 027dc56ffe325652314b4f9e5d804b3b20fd0db3 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Sat, 30 Mar 2019 22:43:26 +0100 Subject: Work around some warnings --- .../astra/ParallelBeamStripKernelProjector2D.inl | 25 +++++++++++----------- 1 file changed, 12 insertions(+), 13 deletions(-) (limited to 'include/astra/ParallelBeamStripKernelProjector2D.inl') diff --git a/include/astra/ParallelBeamStripKernelProjector2D.inl b/include/astra/ParallelBeamStripKernelProjector2D.inl index bed02ab..2031560 100644 --- a/include/astra/ParallelBeamStripKernelProjector2D.inl +++ b/include/astra/ParallelBeamStripKernelProjector2D.inl @@ -147,19 +147,6 @@ void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, const float32 relPixelArea = pixelArea / rayWidth; bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY); - if (vertical) { - RxOverRy = proj->fRayX/proj->fRayY; - deltac = -m_pVolumeGeometry->getPixelLengthY() * RxOverRy * inv_pixelLengthX; - S = 0.5f - 0.5f*fabs(RxOverRy); - T = 0.5f + 0.5f*fabs(RxOverRy); - invTminS = 1.0f / (T-S); - } else { - RyOverRx = proj->fRayY/proj->fRayX; - deltar = -m_pVolumeGeometry->getPixelLengthX() * RyOverRx * inv_pixelLengthY; - S = 0.5f - 0.5f*fabs(RyOverRx); - T = 0.5f + 0.5f*fabs(RyOverRx); - invTminS = 1.0f / (T-S); - } Ex = m_pVolumeGeometry->getWindowMinX() + pixelLengthX*0.5f; Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f; @@ -180,6 +167,12 @@ void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, // vertically if (vertical) { + RxOverRy = proj->fRayX/proj->fRayY; + deltac = -m_pVolumeGeometry->getPixelLengthY() * RxOverRy * inv_pixelLengthX; + S = 0.5f - 0.5f*fabs(RxOverRy); + T = 0.5f + 0.5f*fabs(RxOverRy); + invTminS = 1.0f / (T-S); + // calculate cL and cR for row 0 cL = (DLx + (Ey - DLy)*RxOverRy - Ex) * inv_pixelLengthX; cR = (DRx + (Ey - DRy)*RxOverRy - Ex) * inv_pixelLengthX; @@ -233,6 +226,12 @@ void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, // horizontally else { + RyOverRx = proj->fRayY/proj->fRayX; + deltar = -m_pVolumeGeometry->getPixelLengthX() * RyOverRx * inv_pixelLengthY; + S = 0.5f - 0.5f*fabs(RyOverRx); + T = 0.5f + 0.5f*fabs(RyOverRx); + invTminS = 1.0f / (T-S); + // calculate rL and rR for row 0 rL = -(DLy + (Ex - DLx)*RyOverRx - Ey) * inv_pixelLengthY; rR = -(DRy + (Ex - DRx)*RyOverRx - Ey) * inv_pixelLengthY; -- cgit v1.2.3