Skip to content

Commit

Permalink
Fix FFD integration direction following a change on RTK projectors
Browse files Browse the repository at this point in the history 
RTK integrates X-rays from pixel to source since
RTKConsortium/RTK@a4ba257 for
implementing attenuated projectors.
  • Loading branch information
Simon Rit committed Nov 27, 2023
1 parent 41a3cff commit ddc1885
Showing 1 changed file with 18 additions and 18 deletions.
36 changes: 18 additions & 18 deletions source/digits_hits/include/GateFixedForcedDetectionFunctors.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
{
Expand All @@ -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];
Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -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);
Expand All @@ -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);
}
Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -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;
Expand All @@ -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);
}
Expand Down Expand Up @@ -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];
Expand All @@ -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;
Expand All @@ -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);
}
Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -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;
Expand All @@ -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);
}
Expand Down

0 comments on commit ddc1885

Please sign in to comment.