From ddc188532c306d464ee9ce9989f3d17b7be97e34 Mon Sep 17 00:00:00 2001 From: Simon Rit Date: Mon, 27 Nov 2023 08:11:12 +0100 Subject: [PATCH] Fix FFD integration direction following a change on RTK projectors RTK integrates X-rays from pixel to source since RTKConsortium/RTK@a4ba257d39c2c1f1b7eef79f59631cc7bb2005ae for implementing attenuated projectors. --- .../GateFixedForcedDetectionFunctors.hh | 36 +++++++++---------- 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/source/digits_hits/include/GateFixedForcedDetectionFunctors.hh b/source/digits_hits/include/GateFixedForcedDetectionFunctors.hh index f77ac84d7..b76fac84f 100644 --- a/source/digits_hits/include/GateFixedForcedDetectionFunctors.hh +++ b/source/digits_hits/include/GateFixedForcedDetectionFunctors.hh @@ -543,7 +543,7 @@ namespace GateFixedForcedDetectionFunctor const double & itkNotUsed(rayCastValue), const VectorType & stepInMM, const VectorType & itkNotUsed(source), - const VectorType & sourceToPixel, + const VectorType & pixelToSource, const VectorType & nearestPoint, const VectorType & farthestPoint) { @@ -559,7 +559,7 @@ namespace GateFixedForcedDetectionFunctor /* The last material is the world material. One must fill the weight with the length from source to nearest point and farthest point to pixel point. */ - VectorType worldVector = sourceToPixel + nearestPoint - farthestPoint; + VectorType worldVector = farthestPoint - nearestPoint - pixelToSource; for (int i = 0; i < 3; i++) { worldVector[i] *= m_VolumeSpacing[i]; @@ -651,14 +651,14 @@ namespace GateFixedForcedDetectionFunctor const double &itkNotUsed(rayCastValue), const VectorType &stepInMM, const VectorType &itkNotUsed(source), - const VectorType &sourceToPixel, + const VectorType &pixelToSource, const VectorType &nearestPoint, const VectorType &farthestPoint) { /* Compute ray length in world material This is used to compute the length in world as well as the direction of the ray in mm. */ - VectorType worldVector = sourceToPixel + nearestPoint - farthestPoint; + VectorType worldVector = farthestPoint - nearestPoint - pixelToSource; for (int i = 0; i < 3; i++) { worldVector[i] *= m_VolumeSpacing[i]; @@ -710,7 +710,7 @@ namespace GateFixedForcedDetectionFunctor /* Final computation */ // double weight = std::exp(-rayIntegral) * DCScompton * GetSolidAngle(sourceToPixel) * (*m_ResponseDetector)(energy); - double weight = std::exp(-rayIntegral) * DCScompton * GetSolidAngle(sourceToPixel); + double weight = std::exp(-rayIntegral) * DCScompton * GetSolidAngle(pixelToSource); if (m_generatePhotons) { //Accumulate(threadId, output, weight, m_Energy); @@ -730,7 +730,7 @@ namespace GateFixedForcedDetectionFunctor for (int i = 0; i < 3; i++) { photonDirection[i] = worldVector[i] / worldVectorNorm; - photonPosition[i] = sourceToPixel[i] * m_VolumeSpacing[i]; + photonPosition[i] = -pixelToSource[i] * m_VolumeSpacing[i]; } SavePhotonsparameters(threadId, photonPosition, photonDirection, weight, energy); } @@ -844,12 +844,12 @@ namespace GateFixedForcedDetectionFunctor const double & itkNotUsed(rayCastValue), const VectorType & stepInMM, const VectorType & itkNotUsed(source), - const VectorType & sourceToPixel, + const VectorType & pixelToSource, const VectorType & nearestPoint, const VectorType & farthestPoint) { /* Compute ray length in world material. This is used to compute the length in world as well as the direction of the ray in mm. */ - VectorType worldVector = sourceToPixel + nearestPoint - farthestPoint; + VectorType worldVector = farthestPoint - nearestPoint - pixelToSource; for (int i = 0; i < 3; i++) { worldVector[i] *= m_VolumeSpacing[i]; @@ -886,7 +886,7 @@ namespace GateFixedForcedDetectionFunctor } /* Final computation */ - double weight = std::exp(-rayIntegral) * DCSrayleigh * GetSolidAngle(sourceToPixel); + double weight = std::exp(-rayIntegral) * DCSrayleigh * GetSolidAngle(pixelToSource); if (m_generatePhotons) { VectorType photonDirection; @@ -905,7 +905,7 @@ namespace GateFixedForcedDetectionFunctor for (int i = 0; i < 3; i++) { photonDirection[i] = worldVector[i] / worldVectorNorm; - photonPosition[i] = sourceToPixel[i] * m_VolumeSpacing[i]; + photonPosition[i] = -pixelToSource[i] * m_VolumeSpacing[i]; } SavePhotonsparameters(threadId, photonPosition, photonDirection, weight, m_Energy); } @@ -989,14 +989,14 @@ namespace GateFixedForcedDetectionFunctor const double &itkNotUsed(rayCastValue), const VectorType & stepInMM, const VectorType & itkNotUsed(source), - const VectorType & sourceToPixel, + const VectorType & pixelToSource, const VectorType & nearestPoint, const VectorType & farthestPoint) { /* Compute ray length in world material This is used to compute the length in world as well as the direction of the ray in mm. */ - VectorType worldVector = sourceToPixel + nearestPoint - farthestPoint; + VectorType worldVector = farthestPoint - nearestPoint - pixelToSource; for (int i = 0; i < 3; i++) { worldVector[i] *= m_VolumeSpacing[i]; @@ -1023,7 +1023,7 @@ namespace GateFixedForcedDetectionFunctor } /* Final computation */ - double weight = m_Weight * std::exp(-rayIntegral)*GetSolidAngle(sourceToPixel)/(4*itk::Math::pi); + double weight = m_Weight * std::exp(-rayIntegral)*GetSolidAngle(pixelToSource)/(4*itk::Math::pi); if (m_generatePhotons) { VectorType photonDirection; @@ -1042,7 +1042,7 @@ namespace GateFixedForcedDetectionFunctor for (int i = 0; i < 3; i++) { photonDirection[i] = worldVector[i] / worldVectorNorm; - photonPosition[i] = sourceToPixel[i] * m_VolumeSpacing[i]; + photonPosition[i] = -pixelToSource[i] * m_VolumeSpacing[i]; } SavePhotonsparameters(threadId, photonPosition, photonDirection, weight, m_Energy); } @@ -1095,14 +1095,14 @@ namespace GateFixedForcedDetectionFunctor const double & itkNotUsed(rayCastValue), const VectorType & stepInMM, const VectorType & itkNotUsed(source), - const VectorType & sourceToPixel, + const VectorType & pixelToSource, const VectorType & nearestPoint, const VectorType & farthestPoint) { /* Compute ray length in world material This is used to compute the length in world as well as the direction of the ray in mm. */ - VectorType worldVector = sourceToPixel + nearestPoint - farthestPoint; + VectorType worldVector = farthestPoint - nearestPoint - pixelToSource; for (int i = 0; i < 3; i++) { worldVector[i] *= m_VolumeSpacing[i]; @@ -1136,7 +1136,7 @@ namespace GateFixedForcedDetectionFunctor } /* Final computation */ - double weight = m_Weight * std::exp(-rayIntegral)*GetSolidAngle(sourceToPixel)/(4*itk::Math::pi); + double weight = m_Weight * std::exp(-rayIntegral)*GetSolidAngle(pixelToSource)/(4*itk::Math::pi); if (m_generatePhotons) { VectorType photonDirection; @@ -1155,7 +1155,7 @@ namespace GateFixedForcedDetectionFunctor for (int i = 0; i < 3; i++) { photonDirection[i] = worldVector[i] / worldVectorNorm; - photonPosition[i] = sourceToPixel[i] * m_VolumeSpacing[i]; + photonPosition[i] = -pixelToSource[i] * m_VolumeSpacing[i]; } SavePhotonsparameters(threadId, photonPosition, photonDirection, weight, m_Energy); }