Skip to content

ENH: Lambda functions instead of functor classes in JosephForwardProjectionFilter #694

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
MichaelColonel wants to merge 3 commits into
base: main
Choose a base branch
from
Open

ENH: Lambda functions instead of functor classes in JosephForwardProjectionFilter #694

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
de18d67
ENH: Use lambda functions instead of functor classes in rtkJosephForw…
Feb 28, 2025
a7eea4b
STYLE: Verbose output in tests
Mar 8, 2025
e80a2d7
STYLE: Fix review requests
Mar 13, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
236 changes: 9 additions & 227 deletions include/rtkJosephForwardAttenuatedProjectionImageFilter.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,213 +29,6 @@

namespace rtk
{
namespace Functor
{
/** \class InterpolationWeightMultiplicationAttenuated
* \brief Function to multiply the interpolation weights with the projected
* volume values and attenuation map.
*
* \author Antoine Robert
*
* \ingroup RTK Functions
*/
template <class TInput, class TCoordinateType, class TOutput = TInput>
class ITK_TEMPLATE_EXPORT InterpolationWeightMultiplicationAttenuated
{
public:
InterpolationWeightMultiplicationAttenuated()
{
for (std::size_t i = 0; i < itk::ITK_MAX_THREADS; i++)
{
m_AttenuationRay[i] = 0;
m_AttenuationPixel[i] = 0;
m_Ex1[i] = 1;
}
}

~InterpolationWeightMultiplicationAttenuated() = default;
bool
operator!=(const InterpolationWeightMultiplicationAttenuated &) const
{
return false;
}

bool
operator==(const InterpolationWeightMultiplicationAttenuated & other) const
{
return !(*this != other);
}

inline TOutput
operator()(const ThreadIdType threadId,
const double stepLengthInVoxel,
const TCoordinateType weight,
const TInput * p,
const int i)
{
const double w = weight * stepLengthInVoxel;

m_AttenuationRay[threadId] += w * (p + m_AttenuationMinusEmissionMapsPtrDiff)[i];
m_AttenuationPixel[threadId] += w * (p + m_AttenuationMinusEmissionMapsPtrDiff)[i];
return weight * p[i];
}

void
SetAttenuationMinusEmissionMapsPtrDiff(std::ptrdiff_t pd)
{
m_AttenuationMinusEmissionMapsPtrDiff = pd;
}
TOutput *
GetAttenuationRay()
{
return m_AttenuationRay;
}
TOutput *
GetAttenuationPixel()
{
return m_AttenuationPixel;
}
TOutput *
GetEx1()
{
return m_Ex1;
}

private:
std::ptrdiff_t m_AttenuationMinusEmissionMapsPtrDiff;
TInput m_AttenuationRay[itk::ITK_MAX_THREADS];
TInput m_AttenuationPixel[itk::ITK_MAX_THREADS];
TInput m_Ex1[itk::ITK_MAX_THREADS];
};

/** \class ComputeAttenuationCorrection
* \brief Function to compute the attenuation correction on the projection.
*
* \author Antoine Robert
*
* \ingroup RTK Functions
*/
template <class TInput, class TOutput>
class ITK_TEMPLATE_EXPORT ComputeAttenuationCorrection
{
public:
using VectorType = itk::Vector<double, 3>;

ComputeAttenuationCorrection() = default;
~ComputeAttenuationCorrection() = default;
bool
operator!=(const ComputeAttenuationCorrection &) const
{
return false;
}

bool
operator==(const ComputeAttenuationCorrection & other) const
{
return !(*this != other);
}

inline void
operator()(const ThreadIdType threadId, TOutput & sumValue, const TInput volumeValue, const VectorType & stepInMM)
{
TInput ex2 = exp(-m_AttenuationRay[threadId] * stepInMM.GetNorm());
TInput wf;

if (m_AttenuationPixel[threadId] > 0)
{
wf = (m_Ex1[threadId] - ex2) / m_AttenuationPixel[threadId];
}
else
{
wf = m_Ex1[threadId] * stepInMM.GetNorm();
}

m_Ex1[threadId] = ex2;
m_AttenuationPixel[threadId] = 0;
sumValue += wf * volumeValue;
}

void
SetAttenuationRayVector(TInput * attenuationRayVector)
{
m_AttenuationRay = attenuationRayVector;
}
void
SetAttenuationPixelVector(TInput * attenuationPixelVector)
{
m_AttenuationPixel = attenuationPixelVector;
}
void
SetEx1(TInput * ex1)
{
m_Ex1 = ex1;
}

private:
TInput * m_AttenuationRay;
TInput * m_AttenuationPixel;
TInput * m_Ex1;
};

/** \class ProjectedValueAccumulationAttenuated
* \brief Function to accumulate the ray casting on the projection.
*
* \author Antoine Robert
*
* \ingroup RTK Functions
*/
template <class TInput, class TOutput>
class ITK_TEMPLATE_EXPORT ProjectedValueAccumulationAttenuated
{
public:
using VectorType = itk::Vector<double, 3>;

ProjectedValueAccumulationAttenuated() = default;
~ProjectedValueAccumulationAttenuated() = default;
bool
operator!=(const ProjectedValueAccumulationAttenuated &) const
{
return false;
}

bool
operator==(const ProjectedValueAccumulationAttenuated & other) const
{
return !(*this != other);
}

inline void
operator()(const ThreadIdType threadId,
const TInput & input,
TOutput & output,
const TOutput & rayCastValue,
const VectorType & /*stepInMM*/,
const VectorType & itkNotUsed(source),
const VectorType & itkNotUsed(sourceToPixel),
const VectorType & itkNotUsed(nearestPoint),
const VectorType & itkNotUsed(farthestPoint))
{
output = input + rayCastValue;
m_Attenuation[threadId] = 0;
m_Ex1[threadId] = 1;
}

void
SetAttenuationVector(TInput * attenuationVector)
{
m_Attenuation = attenuationVector;
}
void
SetEx1(TInput * ex1)
{
m_Ex1 = ex1;
}

private:
TInput * m_Attenuation;
TInput * m_Ex1;
};
} // end namespace Functor

/** \class JosephForwardAttenuatedProjectionImageFilter
* \brief Joseph forward projection.
Expand All @@ -252,40 +45,24 @@ class ITK_TEMPLATE_EXPORT ProjectedValueAccumulationAttenuated
* \ingroup RTK Projector
*/

template <
class TInputImage,
class TOutputImage,
class TInterpolationWeightMultiplication = Functor::InterpolationWeightMultiplicationAttenuated<
typename TInputImage::PixelType,
typename itk::PixelTraits<typename TInputImage::PixelType>::ValueType>,
class TProjectedValueAccumulation =
Functor::ProjectedValueAccumulationAttenuated<typename TInputImage::PixelType, typename TOutputImage::PixelType>,
class TSumAlongRay =
Functor::ComputeAttenuationCorrection<typename TInputImage::PixelType, typename TOutputImage::PixelType>>
template <class TInputImage, class TOutputImage>
class ITK_TEMPLATE_EXPORT JosephForwardAttenuatedProjectionImageFilter
: public JosephForwardProjectionImageFilter<TInputImage,
TOutputImage,
TInterpolationWeightMultiplication,
TProjectedValueAccumulation,
TSumAlongRay>
: public JosephForwardProjectionImageFilter<TInputImage, TOutputImage>
{
public:
ITK_DISALLOW_COPY_AND_MOVE(JosephForwardAttenuatedProjectionImageFilter);

/** Standard class type alias. */
using Self = JosephForwardAttenuatedProjectionImageFilter;
using Superclass = JosephForwardProjectionImageFilter<TInputImage,
TOutputImage,
TInterpolationWeightMultiplication,
TProjectedValueAccumulation,
TSumAlongRay>;
using Superclass = JosephForwardProjectionImageFilter<TInputImage, TOutputImage>;
using Pointer = itk::SmartPointer<Self>;
using ConstPointer = itk::SmartPointer<const Self>;
using InputPixelType = typename TInputImage::PixelType;
using OutputPixelType = typename TOutputImage::PixelType;
using OutputImageRegionType = typename TOutputImage::RegionType;
using CoordinateType = double;
using VectorType = itk::Vector<CoordinateType, TInputImage::ImageDimension>;
using WeightCoordinateType = typename itk::PixelTraits<InputPixelType>::ValueType;

/** ImageDimension constants */
static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension;
Expand All @@ -311,6 +88,11 @@ class ITK_TEMPLATE_EXPORT JosephForwardAttenuatedProjectionImageFilter
* to overwrite the method. */
void
VerifyInputInformation() const override;

std::ptrdiff_t m_AttenuationMinusEmissionMapsPtrDiff;
std::array<InputPixelType, itk::ITK_MAX_THREADS> m_AttenuationRay;
std::array<InputPixelType, itk::ITK_MAX_THREADS> m_AttenuationPixel;
std::array<InputPixelType, itk::ITK_MAX_THREADS> m_Ex1;
};
} // end namespace rtk

Expand Down
Loading
Loading