From cbbaabb81b5cceadea3c03bd0ebe0ade6f458645 Mon Sep 17 00:00:00 2001 From: Mikhail Polkovnikov Date: Fri, 28 Feb 2025 13:06:55 +0300 Subject: [PATCH 1/2] ENH: Use lambda functions instead of functor classes in rtkJosephForwardProjectionFilter Lambdas were added to the classes: JosephForwardProjectionImageFilter, JosephForwardAttenuatedProjectionImageFilter, MaximumIntensityProjectionImageFilter. Python wrappers for classes above have been modified. Similar lambda functions can be added to back projection classes. --- ...phForwardAttenuatedProjectionImageFilter.h | 238 +----------------- ...ForwardAttenuatedProjectionImageFilter.hxx | 129 ++++++---- .../rtkJosephForwardProjectionImageFilter.h | 215 ++++------------ .../rtkJosephForwardProjectionImageFilter.hxx | 154 +++++------- ...rtkMaximumIntensityProjectionImageFilter.h | 122 ++------- ...kMaximumIntensityProjectionImageFilter.hxx | 60 +++++ test/rtkforwardprojectiontest.cxx | 11 + test/rtkmaximumintensityprojectiontest.cxx | 87 ++++++- ...orwardAttenuatedProjectionImageFilter.wrap | 25 -- ...MaximumIntensityProjectionImageFilter.wrap | 27 +- 10 files changed, 369 insertions(+), 699 deletions(-) create mode 100644 include/rtkMaximumIntensityProjectionImageFilter.hxx diff --git a/include/rtkJosephForwardAttenuatedProjectionImageFilter.h b/include/rtkJosephForwardAttenuatedProjectionImageFilter.h index 0282f0967..02012d584 100644 --- a/include/rtkJosephForwardAttenuatedProjectionImageFilter.h +++ b/include/rtkJosephForwardAttenuatedProjectionImageFilter.h @@ -29,214 +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 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); - } - - 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 ITK_TEMPLATE_EXPORT ComputeAttenuationCorrection -{ -public: - using VectorType = itk::Vector; - - ComputeAttenuationCorrection() = default; - ~ComputeAttenuationCorrection() = default; - bool - operator!=(const ComputeAttenuationCorrection &) const - { - return false; - } - - bool - operator==(const ComputeAttenuationCorrection & other) const - { - return !(*this != other); - } - - 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 ITK_TEMPLATE_EXPORT ProjectedValueAccumulationAttenuated -{ -public: - using VectorType = itk::Vector; - using PointType = itk::Point; - - ProjectedValueAccumulationAttenuated() = default; - ~ProjectedValueAccumulationAttenuated() = default; - bool - operator!=(const ProjectedValueAccumulationAttenuated &) const - { - return false; - } - - bool - operator==(const ProjectedValueAccumulationAttenuated & other) const - { - return !(*this != other); - } - - void - operator()(const ThreadIdType threadId, - const TInput & input, - TOutput & output, - const TOutput & rayCastValue, - const VectorType & /*stepInMM*/, - const PointType & itkNotUsed(source), - const VectorType & itkNotUsed(sourceToPixel), - const PointType & itkNotUsed(nearestPoint), - const PointType & 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. @@ -253,40 +45,25 @@ class ITK_TEMPLATE_EXPORT ProjectedValueAccumulationAttenuated * \ingroup RTK Projector */ -template < - class TInputImage, - class TOutputImage, - class TInterpolationWeightMultiplication = Functor::InterpolationWeightMultiplicationAttenuated< - typename TInputImage::PixelType, - typename itk::PixelTraits::ValueType>, - class TProjectedValueAccumulation = - Functor::ProjectedValueAccumulationAttenuated, - class TSumAlongRay = - Functor::ComputeAttenuationCorrection> +template class ITK_TEMPLATE_EXPORT JosephForwardAttenuatedProjectionImageFilter - : public JosephForwardProjectionImageFilter + : public JosephForwardProjectionImageFilter { public: ITK_DISALLOW_COPY_AND_MOVE(JosephForwardAttenuatedProjectionImageFilter); /** Standard class type alias. */ using Self = JosephForwardAttenuatedProjectionImageFilter; - using Superclass = JosephForwardProjectionImageFilter; + using Superclass = JosephForwardProjectionImageFilter; using Pointer = itk::SmartPointer; using ConstPointer = itk::SmartPointer; using InputPixelType = typename TInputImage::PixelType; using OutputPixelType = typename TOutputImage::PixelType; using OutputImageRegionType = typename TOutputImage::RegionType; + using PointType = typename TInputImage::PointType; using CoordinateType = double; using VectorType = itk::Vector; + using WeightCoordinateType = typename itk::PixelTraits::ValueType; /** ImageDimension constants */ static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension; @@ -312,6 +89,11 @@ class ITK_TEMPLATE_EXPORT JosephForwardAttenuatedProjectionImageFilter * to overwrite the method. */ void VerifyInputInformation() const override; + + std::ptrdiff_t m_AttenuationMinusEmissionMapsPtrDiff; + std::array m_AttenuationRay; + std::array m_AttenuationPixel; + std::array m_Ex1; }; } // end namespace rtk diff --git a/include/rtkJosephForwardAttenuatedProjectionImageFilter.hxx b/include/rtkJosephForwardAttenuatedProjectionImageFilter.hxx index 545343554..19df360c6 100644 --- a/include/rtkJosephForwardAttenuatedProjectionImageFilter.hxx +++ b/include/rtkJosephForwardAttenuatedProjectionImageFilter.hxx @@ -31,33 +31,82 @@ namespace rtk { -template -JosephForwardAttenuatedProjectionImageFilter< - TInputImage, - TOutputImage, - TInterpolationWeightMultiplication, - TProjectedValueAccumulation, - TComputeAttenuationCorrection>::JosephForwardAttenuatedProjectionImageFilter() +template +JosephForwardAttenuatedProjectionImageFilter::JosephForwardAttenuatedProjectionImageFilter() + : JosephForwardProjectionImageFilter() { + std::fill(m_AttenuationRay.begin(), m_AttenuationRay.end(), 0.); + std::fill(m_AttenuationPixel.begin(), m_AttenuationPixel.end(), 0.); + std::fill(m_Ex1.begin(), m_Ex1.end(), 1.); + m_AttenuationMinusEmissionMapsPtrDiff = 0; + + /** \brief Function to multiply the interpolation weights with the projected + * volume values and attenuation map. + * + */ + auto interpolationWeightMultiplicationAttenuatedFunc = [this](const ThreadIdType threadId, + double stepLengthInVoxel, + const WeightCoordinateType weight, + const InputPixelType * p, + const int i) -> OutputPixelType { + const double w = weight * stepLengthInVoxel; + + this->m_AttenuationRay[threadId] += w * (p + this->m_AttenuationMinusEmissionMapsPtrDiff)[i]; + this->m_AttenuationPixel[threadId] += w * (p + this->m_AttenuationMinusEmissionMapsPtrDiff)[i]; + return weight * p[i]; + }; + this->SetInterpolationWeightMultiplication(interpolationWeightMultiplicationAttenuatedFunc); + + /** \brief Function to compute the attenuation correction on the projection. + * + */ + auto computeAttenuationCorrectionFunc = [this](const ThreadIdType threadId, + OutputPixelType & sumValue, + const InputPixelType volumeValue, + const VectorType & stepInMM) { + InputPixelType ex2 = exp(-1. * this->m_AttenuationRay[threadId] * stepInMM.GetNorm()); + InputPixelType wf; + + if (this->m_AttenuationPixel[threadId] > 0) + { + wf = (this->m_Ex1[threadId] - ex2) / this->m_AttenuationPixel[threadId]; + } + else + { + wf = this->m_Ex1[threadId] * stepInMM.GetNorm(); + } + + this->m_Ex1[threadId] = ex2; + this->m_AttenuationPixel[threadId] = 0; + sumValue += wf * volumeValue; + }; + this->SetSumAlongRay(computeAttenuationCorrectionFunc); + + /** \brief Function to accumulate the ray casting on the projection. + * + */ + auto projectedValueAccumulationAttenuatedFunc = [this](const ThreadIdType threadId, + const InputPixelType & input, + OutputPixelType & output, + const OutputPixelType & rayCastValue, + const VectorType & itkNotUsed(stepInMM), + const PointType & itkNotUsed(source), + const VectorType & itkNotUsed(sourceToPixel), + const PointType & itkNotUsed(nearestPoint), + const PointType & itkNotUsed(farthestPoint)) { + output = input + rayCastValue; + this->m_AttenuationRay[threadId] = 0; + this->m_Ex1[threadId] = 1; + }; + this->SetProjectedValueAccumulation(projectedValueAccumulationAttenuatedFunc); + this->SetNumberOfRequiredInputs(3); this->DynamicMultiThreadingOff(); } -template +template void -JosephForwardAttenuatedProjectionImageFilter::GenerateInputRequestedRegion() +JosephForwardAttenuatedProjectionImageFilter::GenerateInputRequestedRegion() { Superclass::GenerateInputRequestedRegion(); // Input 2 is the attenuation map relative to the volume @@ -69,17 +118,9 @@ JosephForwardAttenuatedProjectionImageFilterSetRequestedRegion(reqRegion2); } -template +template void -JosephForwardAttenuatedProjectionImageFilter::VerifyInputInformation() const +JosephForwardAttenuatedProjectionImageFilter::VerifyInputInformation() const { Superclass::VerifyInputInformation(); using ImageBaseType = const itk::ImageBase; @@ -160,30 +201,12 @@ JosephForwardAttenuatedProjectionImageFilter +template void -JosephForwardAttenuatedProjectionImageFilter::BeforeThreadedGenerateData() +JosephForwardAttenuatedProjectionImageFilter::BeforeThreadedGenerateData() { - if (!this->GetInput(2)) - { - itkExceptionMacro("Attenuation map (input 2) must be set before running the attenuated forward projector."); - } - this->GetInterpolationWeightMultiplication().SetAttenuationMinusEmissionMapsPtrDiff( - this->GetInput(2)->GetBufferPointer() - this->GetInput(1)->GetBufferPointer()); - this->GetProjectedValueAccumulation().SetAttenuationVector( - this->GetInterpolationWeightMultiplication().GetAttenuationRay()); - this->GetSumAlongRay().SetAttenuationRayVector(this->GetInterpolationWeightMultiplication().GetAttenuationRay()); - this->GetSumAlongRay().SetAttenuationPixelVector(this->GetInterpolationWeightMultiplication().GetAttenuationPixel()); - this->GetProjectedValueAccumulation().SetEx1(this->GetInterpolationWeightMultiplication().GetEx1()); - this->GetSumAlongRay().SetEx1(this->GetInterpolationWeightMultiplication().GetEx1()); + this->m_AttenuationMinusEmissionMapsPtrDiff = + this->GetInput(2)->GetBufferPointer() - this->GetInput(1)->GetBufferPointer(); } } // end namespace rtk diff --git a/include/rtkJosephForwardProjectionImageFilter.h b/include/rtkJosephForwardProjectionImageFilter.h index d8f9798ec..e6c394a9b 100644 --- a/include/rtkJosephForwardProjectionImageFilter.h +++ b/include/rtkJosephForwardProjectionImageFilter.h @@ -29,124 +29,6 @@ #include namespace rtk { -namespace Functor -{ -/** \class InterpolationWeightMultiplication - * \brief Function to multiply the interpolation weights with the projected - * volume values. - * - * \author Simon Rit - * - * \ingroup RTK Functions - */ -template -class ITK_TEMPLATE_EXPORT InterpolationWeightMultiplication -{ -public: - InterpolationWeightMultiplication() = default; - ~InterpolationWeightMultiplication() = default; - bool - operator!=(const InterpolationWeightMultiplication &) const - { - return false; - } - bool - operator==(const InterpolationWeightMultiplication & other) const - { - return !(*this != other); - } - - TOutput - operator()(const ThreadIdType itkNotUsed(threadId), - const double itkNotUsed(stepLengthInVoxel), - const TCoordinateType weight, - const TInput * p, - const int i) const - { - return weight * p[i]; - } -}; - -/** \class SumAlongRay - * \brief Function to compute the attenuation correction on the projection. - * - * \author Antoine Robert - * - * \ingroup RTK Functions - */ -template -class ITK_TEMPLATE_EXPORT SumAlongRay -{ -public: - using VectorType = itk::Vector; - - SumAlongRay() = default; - ~SumAlongRay() = default; - bool - operator!=(const SumAlongRay &) const - { - return false; - } - bool - operator==(const SumAlongRay & other) const - { - return !(*this != other); - } - - void - operator()(const ThreadIdType itkNotUsed(threadId), - TOutput & sumValue, - const TInput volumeValue, - const VectorType & itkNotUsed(stepInMM)) - { - sumValue += static_cast(volumeValue); - } -}; - -/** \class ProjectedValueAccumulation - * \brief Function to accumulate the ray casting on the projection. - * - * \author Simon Rit - * - * \ingroup RTK Functions - */ -template -class ITK_TEMPLATE_EXPORT ProjectedValueAccumulation -{ -public: - using VectorType = itk::Vector; - using PointType = itk::Point; - - ProjectedValueAccumulation() = default; - ~ProjectedValueAccumulation() = default; - bool - operator!=(const ProjectedValueAccumulation &) const - { - return false; - } - bool - operator==(const ProjectedValueAccumulation & other) const - { - return !(*this != other); - } - - void - operator()(const ThreadIdType itkNotUsed(threadId), - const TInput & input, - TOutput & output, - const TOutput & rayCastValue, - const VectorType & stepInMM, - const PointType & itkNotUsed(source), - const VectorType & itkNotUsed(sourceToPixel), - const PointType & itkNotUsed(nearestPoint), - const PointType & itkNotUsed(farthestPoint)) const - { - output = input + rayCastValue * stepInMM.GetNorm(); - } -}; - -} // end namespace Functor - /** \class JosephForwardProjectionImageFilter * \brief Joseph forward projection. @@ -162,15 +44,7 @@ class ITK_TEMPLATE_EXPORT ProjectedValueAccumulation * * \ingroup RTK Projector */ - -template ::ValueType>, - class TProjectedValueAccumulation = - Functor::ProjectedValueAccumulation, - class TSumAlongRay = Functor::SumAlongRay> +template class ITK_TEMPLATE_EXPORT JosephForwardProjectionImageFilter : public ForwardProjectionImageFilter { @@ -185,79 +59,98 @@ class ITK_TEMPLATE_EXPORT JosephForwardProjectionImageFilter using InputPixelType = typename TInputImage::PixelType; using OutputPixelType = typename TOutputImage::PixelType; using OutputImageRegionType = typename TOutputImage::RegionType; + using PointType = typename TInputImage::PointType; using CoordinateType = double; + using WeightCoordinateType = typename itk::PixelTraits::ValueType; using VectorType = itk::Vector; using TClipImageType = itk::Image; using TClipImagePointer = typename TClipImageType::Pointer; + /** \brief Function to multiply the interpolation weights with the projected + * volume values. + * + */ + using InterpolationWeightMultiplicationFunc = + std::function; + + /** \brief Function to compute the attenuation correction on the projection. + * + */ + using SumAlongRayFunc = + std::function; + + /** \brief Function to accumulate the ray casting on the projection. + * + */ + using ProjectedValueAccumulationFunc = std::function; + /** Method for creation through the object factory. */ itkNewMacro(Self); /** Run-time type information (and related methods). */ itkOverrideGetNameOfClassMacro(JosephForwardProjectionImageFilter); - /** Get/Set the functor that is used to multiply each interpolation value with a volume value */ - TInterpolationWeightMultiplication & + /** Get/Set the lambda function that is used to multiply each interpolation value with a volume value */ + InterpolationWeightMultiplicationFunc & GetInterpolationWeightMultiplication() { return m_InterpolationWeightMultiplication; } - const TInterpolationWeightMultiplication & + const InterpolationWeightMultiplicationFunc & GetInterpolationWeightMultiplication() const { return m_InterpolationWeightMultiplication; } void - SetInterpolationWeightMultiplication(const TInterpolationWeightMultiplication & _arg) + SetInterpolationWeightMultiplication(const InterpolationWeightMultiplicationFunc & _arg) { - if (m_InterpolationWeightMultiplication != _arg) - { - m_InterpolationWeightMultiplication = _arg; - this->Modified(); - } + m_InterpolationWeightMultiplication = _arg; + this->Modified(); } - /** Get/Set the functor that is used to accumulate values in the projection image after the ray + /** Get/Set the lambda function that is used to accumulate values in the projection image after the ray * casting has been performed. */ - TProjectedValueAccumulation & + ProjectedValueAccumulationFunc & GetProjectedValueAccumulation() { return m_ProjectedValueAccumulation; } - const TProjectedValueAccumulation & + const ProjectedValueAccumulationFunc & GetProjectedValueAccumulation() const { return m_ProjectedValueAccumulation; } void - SetProjectedValueAccumulation(const TProjectedValueAccumulation & _arg) + SetProjectedValueAccumulation(const ProjectedValueAccumulationFunc & _arg) { - if (m_ProjectedValueAccumulation != _arg) - { - m_ProjectedValueAccumulation = _arg; - this->Modified(); - } + m_ProjectedValueAccumulation = _arg; + this->Modified(); } - /** Get/Set the functor that is used to compute the sum along the ray*/ - TSumAlongRay & + /** Get/Set the lambda function that is used to compute the sum along the ray */ + SumAlongRayFunc & GetSumAlongRay() { return m_SumAlongRay; } - const TSumAlongRay & + const SumAlongRayFunc & GetSumAlongRay() const { return m_SumAlongRay; } void - SetSumAlongRay(const TSumAlongRay & _arg) + SetSumAlongRay(const SumAlongRayFunc & _arg) { - if (m_SumAlongRay != _arg) - { - m_SumAlongRay = _arg; - this->Modified(); - } + m_SumAlongRay = _arg; + this->Modified(); } /** Set/Get the inferior clip image. Each pixel of the image @@ -343,12 +236,14 @@ class ITK_TEMPLATE_EXPORT JosephForwardProjectionImageFilter double maxy); private: - // Functors - TInterpolationWeightMultiplication m_InterpolationWeightMultiplication; - TProjectedValueAccumulation m_ProjectedValueAccumulation; - TSumAlongRay m_SumAlongRay; - double m_InferiorClip{ 0. }; - double m_SuperiorClip{ 1. }; + // lambdas or std::functions + // MUST BE initiated in constructor, or custom functions + // must be set by means of Set methods in application + InterpolationWeightMultiplicationFunc m_InterpolationWeightMultiplication; + SumAlongRayFunc m_SumAlongRay; + ProjectedValueAccumulationFunc m_ProjectedValueAccumulation; + double m_InferiorClip{ 0. }; + double m_SuperiorClip{ 1. }; }; } // end namespace rtk diff --git a/include/rtkJosephForwardProjectionImageFilter.hxx b/include/rtkJosephForwardProjectionImageFilter.hxx index 97f01881e..458927d35 100644 --- a/include/rtkJosephForwardProjectionImageFilter.hxx +++ b/include/rtkJosephForwardProjectionImageFilter.hxx @@ -31,32 +31,37 @@ namespace rtk { -template -JosephForwardProjectionImageFilter::JosephForwardProjectionImageFilter() +template +JosephForwardProjectionImageFilter::JosephForwardProjectionImageFilter() { + this->m_InterpolationWeightMultiplication = [](const ThreadIdType itkNotUsed(threadId), + double itkNotUsed(stepLengthInVoxel), + const WeightCoordinateType weight, + const InputPixelType * p, + int i) -> OutputPixelType { return weight * p[i]; }; + + this->m_SumAlongRay = + [](const ThreadIdType, OutputPixelType & sumValue, const InputPixelType volumeValue, const VectorType &) { + sumValue += static_cast(volumeValue); + }; + + this->m_ProjectedValueAccumulation = [](const ThreadIdType, + const InputPixelType & input, + OutputPixelType & output, + const OutputPixelType & rayCastValue, + const VectorType & stepInMM, + const PointType &, + const VectorType &, + const PointType &, + const PointType &) { output = input + rayCastValue * stepInMM.GetNorm(); }; + this->DynamicMultiThreadingOff(); } -template +template void -JosephForwardProjectionImageFilter::GenerateInputRequestedRegion() +JosephForwardProjectionImageFilter::GenerateInputRequestedRegion() { Superclass::GenerateInputRequestedRegion(); @@ -78,17 +83,9 @@ JosephForwardProjectionImageFilter +template void -JosephForwardProjectionImageFilter::VerifyInputInformation() const +JosephForwardProjectionImageFilter::VerifyInputInformation() const { using ImageBaseType = const itk::ImageBase; @@ -169,19 +166,11 @@ JosephForwardProjectionImageFilter +template void -JosephForwardProjectionImageFilter::ThreadedGenerateData(const OutputImageRegionType & - outputRegionForThread, - ThreadIdType threadId) +JosephForwardProjectionImageFilter::ThreadedGenerateData( + const OutputImageRegionType & outputRegionForThread, + ThreadIdType threadId) { const unsigned int Dimension = TInputImage::ImageDimension; int offsets[3]; @@ -426,30 +415,18 @@ JosephForwardProjectionImageFilter -typename JosephForwardProjectionImageFilter::OutputPixelType -JosephForwardProjectionImageFilter::BilinearInterpolation(const ThreadIdType threadId, - const double stepLengthInVoxel, - const InputPixelType * pxiyi, - const InputPixelType * pxsyi, - const InputPixelType * pxiys, - const InputPixelType * pxsys, - const CoordinateType x, - const CoordinateType y, - const int ox, - const int oy) +template +typename JosephForwardProjectionImageFilter::OutputPixelType +JosephForwardProjectionImageFilter::BilinearInterpolation(const ThreadIdType threadId, + const double stepLengthInVoxel, + const InputPixelType * pxiyi, + const InputPixelType * pxsyi, + const InputPixelType * pxiys, + const InputPixelType * pxsys, + const CoordinateType x, + const CoordinateType y, + const int ox, + const int oy) { int ix = itk::Math::floor(x); int iy = itk::Math::floor(y); @@ -476,34 +453,23 @@ JosephForwardProjectionImageFilter -typename JosephForwardProjectionImageFilter::OutputPixelType -JosephForwardProjectionImageFilter::BilinearInterpolationOnBorders(const ThreadIdType threadId, - const double stepLengthInVoxel, - const InputPixelType * pxiyi, - const InputPixelType * pxsyi, - const InputPixelType * pxiys, - const InputPixelType * pxsys, - const CoordinateType x, - const CoordinateType y, - const int ox, - const int oy, - const CoordinateType minx, - const CoordinateType miny, - const CoordinateType maxx, - const CoordinateType maxy) +template +typename JosephForwardProjectionImageFilter::OutputPixelType +JosephForwardProjectionImageFilter::BilinearInterpolationOnBorders( + const ThreadIdType threadId, + const double stepLengthInVoxel, + const InputPixelType * pxiyi, + const InputPixelType * pxsyi, + const InputPixelType * pxiys, + const InputPixelType * pxsys, + const CoordinateType x, + const CoordinateType y, + const int ox, + const int oy, + const CoordinateType minx, + const CoordinateType miny, + const CoordinateType maxx, + const CoordinateType maxy) { int ix = itk::Math::floor(x); int iy = itk::Math::floor(y); diff --git a/include/rtkMaximumIntensityProjectionImageFilter.h b/include/rtkMaximumIntensityProjectionImageFilter.h index bf46231b5..23e7ad94f 100644 --- a/include/rtkMaximumIntensityProjectionImageFilter.h +++ b/include/rtkMaximumIntensityProjectionImageFilter.h @@ -23,96 +23,6 @@ namespace rtk { -namespace Functor -{ - -/** \class MaximumIntensityAlongRay - * \brief Function to compute the maximum intensity (MIP) value along the ray projection. - * - * \author Mikhail Polkovnikov - * - * \ingroup RTK Functions - */ -template -class ITK_TEMPLATE_EXPORT MaximumIntensityAlongRay -{ -public: - using VectorType = itk::Vector; - - MaximumIntensityAlongRay() = default; - ~MaximumIntensityAlongRay() = default; - bool - operator!=(const MaximumIntensityAlongRay &) const - { - return false; - } - bool - operator==(const MaximumIntensityAlongRay & other) const - { - return !(*this != other); - } - - void - operator()(const ThreadIdType itkNotUsed(threadId), - TOutput & mipValue, - const TInput volumeValue, - const VectorType & itkNotUsed(stepInMM)) - { - auto tmp = static_cast(volumeValue); - if (tmp > mipValue) - { - mipValue = tmp; - } - } -}; - -/** \class MaximumIntensityProjectedValueAccumulation - * \brief Function to calculate maximum intensity step along the ray projection. - * - * \author Mikhail Polkovnikov - * - * \ingroup RTK Functions - */ -template -class ITK_TEMPLATE_EXPORT MaximumIntensityProjectedValueAccumulation -{ -public: - using VectorType = itk::Vector; - using PointType = itk::Point; - - bool - operator!=(const MaximumIntensityProjectedValueAccumulation &) const - { - return false; - } - bool - operator==(const MaximumIntensityProjectedValueAccumulation & other) const - { - return !(*this != other); - } - - void - operator()(const ThreadIdType itkNotUsed(threadId), - const TInput & input, - TOutput & output, - const TOutput & rayCastValue, - const VectorType & stepInMM, - const PointType & itkNotUsed(source), - const VectorType & itkNotUsed(sourceToPixel), - const PointType & itkNotUsed(nearestPoint), - const PointType & itkNotUsed(farthestPoint)) const - { - auto tmp = static_cast(input); - if (tmp < rayCastValue) - { - tmp = rayCastValue; - } - output = tmp * stepInMM.GetNorm(); - } -}; - -} // end namespace Functor - /** \class MaximumIntensityProjectionImageFilter * \brief MIP filter. @@ -120,27 +30,16 @@ class ITK_TEMPLATE_EXPORT MaximumIntensityProjectedValueAccumulation * Performs a MIP forward projection, i.e. calculation of a maximum intensity * step along the x-ray line. * + * \test rtkmaximumintensityprojectiontest.cxx + * * \author Mikhail Polkovnikov * * \ingroup RTK Projector */ -template ::ValueType>, - class TProjectedValueAccumulation = - Functor::MaximumIntensityProjectedValueAccumulation, - class TSumAlongRay = - Functor::MaximumIntensityAlongRay> +template class ITK_TEMPLATE_EXPORT MaximumIntensityProjectionImageFilter - : public JosephForwardProjectionImageFilter + : public JosephForwardProjectionImageFilter { public: ITK_DISALLOW_COPY_AND_MOVE(MaximumIntensityProjectionImageFilter); @@ -148,6 +47,13 @@ class ITK_TEMPLATE_EXPORT MaximumIntensityProjectionImageFilter /** Standard class type alias. */ using Self = MaximumIntensityProjectionImageFilter; using Pointer = itk::SmartPointer; + using Superclass = JosephForwardProjectionImageFilter; + using ConstPointer = itk::SmartPointer; + using InputPixelType = typename TInputImage::PixelType; + using OutputPixelType = typename TOutputImage::PixelType; + using PointType = typename TInputImage::PointType; + using CoordinateType = double; + using VectorType = itk::Vector; /** Method for creation through the object factory. */ itkNewMacro(Self); @@ -156,9 +62,13 @@ class ITK_TEMPLATE_EXPORT MaximumIntensityProjectionImageFilter itkOverrideGetNameOfClassMacro(MaximumIntensityProjectionImageFilter); protected: - MaximumIntensityProjectionImageFilter() = default; + MaximumIntensityProjectionImageFilter(); ~MaximumIntensityProjectionImageFilter() override = default; }; } // end namespace rtk +#ifndef ITK_MANUAL_INSTANTIATION +# include "rtkMaximumIntensityProjectionImageFilter.hxx" +#endif + #endif diff --git a/include/rtkMaximumIntensityProjectionImageFilter.hxx b/include/rtkMaximumIntensityProjectionImageFilter.hxx new file mode 100644 index 000000000..05b64b2e3 --- /dev/null +++ b/include/rtkMaximumIntensityProjectionImageFilter.hxx @@ -0,0 +1,60 @@ +/*========================================================================= + * + * Copyright RTK Consortium + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#ifndef rtkMaximumIntensityProjectionImageFilter_hxx +#define rtkMaximumIntensityProjectionImageFilter_hxx + +namespace rtk +{ + +template +MaximumIntensityProjectionImageFilter::MaximumIntensityProjectionImageFilter() + : JosephForwardProjectionImageFilter() +{ + auto sumAlongRayFunc = + [](const ThreadIdType, OutputPixelType & mipValue, const InputPixelType volumeValue, const VectorType &) { + OutputPixelType tmp = static_cast(volumeValue); + if (tmp > mipValue) + { + mipValue = tmp; + } + }; + this->SetSumAlongRay(sumAlongRayFunc); + + auto projectedValueAccumulationFunc = [](const ThreadIdType, + const InputPixelType & input, + OutputPixelType & output, + const OutputPixelType & rayCastValue, + const VectorType & stepInMM, + const PointType &, + const VectorType &, + const PointType &, + const PointType &) { + OutputPixelType tmp = static_cast(input); + if (tmp < rayCastValue) + { + tmp = rayCastValue; + } + output = tmp * stepInMM.GetNorm(); + }; + this->SetProjectedValueAccumulation(projectedValueAccumulationFunc); +} + +} // end namespace rtk + +#endif diff --git a/test/rtkforwardprojectiontest.cxx b/test/rtkforwardprojectiontest.cxx index 4d56c7b13..00f99fede 100644 --- a/test/rtkforwardprojectiontest.cxx +++ b/test/rtkforwardprojectiontest.cxx @@ -83,6 +83,17 @@ rtkforwardprojectiontest(int, char *[]) jfp->SetInput(projInput->GetOutput()); jfp->SetInput(1, volInput->GetOutput()); +#ifndef USE_CUDA + // Adding custom SumAlongRay lambda function doesn't work with CUDA filter + JFPType::SumAlongRayFunc sumFuncTest = [](const itk::ThreadIdType, + JFPType::OutputPixelType & sum, + const JFPType::InputPixelType in, + const JFPType::VectorType &) -> void { + sum += static_cast(in); + }; + jfp->SetSumAlongRay(sumFuncTest); +#endif + // Ray Box Intersection filter (reference) #ifdef USE_CUDA jfp->SetStepSize(10); diff --git a/test/rtkmaximumintensityprojectiontest.cxx b/test/rtkmaximumintensityprojectiontest.cxx index 2f8229ada..714278d65 100644 --- a/test/rtkmaximumintensityprojectiontest.cxx +++ b/test/rtkmaximumintensityprojectiontest.cxx @@ -1,5 +1,6 @@ #include "rtkTest.h" #include "rtkConstantImageSource.h" +#include "rtkJosephForwardProjectionImageFilter.h" #include "rtkMaximumIntensityProjectionImageFilter.h" #include "rtkRayBoxIntersectionImageFilter.h" #include "rtkThreeDCircularProjectionGeometryXMLFile.h" @@ -54,18 +55,54 @@ rtkmaximumintensityprojectiontest(int, char *[]) projInput->SetConstant(0.); projInput->Update(); - // MIP Forward Projection filter - auto mipfp = rtk::MaximumIntensityProjectionImageFilter::New(); - mipfp->SetInput(projInput->GetOutput()); - mipfp->SetInput(1, volInput->GetOutput()); + // Test-1 + // Joseph Forward Projection filter for MIP projection + std::cout << "\n\n****** Case 1: JosephForwardProjection filter with custom MIP functions ******" << std::endl; + using JFPType = rtk::JosephForwardProjectionImageFilter; + auto jfp = rtk::JosephForwardProjectionImageFilter::New(); + jfp->SetInput(projInput->GetOutput()); + jfp->SetInput(1, volInput->GetOutput()); + + // Function to compute the maximum intensity (MIP) value along the ray projection. + JFPType::SumAlongRayFunc sumAlongFunc = [](const itk::ThreadIdType, + JFPType::OutputPixelType & mipValue, + const JFPType::InputPixelType volumeValue, + const JFPType::VectorType &) -> void { + JFPType::OutputPixelType tmp = static_cast(volumeValue); + if (tmp > mipValue) + { + mipValue = tmp; + } + }; + jfp->SetSumAlongRay(sumAlongFunc); + + // Performs a MIP forward projection, i.e. calculation of a maximum intensity + // step along the x-ray line. + JFPType::ProjectedValueAccumulationFunc projAccumFunc = [](const itk::ThreadIdType itkNotUsed(threadId), + const JFPType::InputPixelType & input, + JFPType::OutputPixelType & output, + const JFPType::OutputPixelType & rayCastValue, + const JFPType::VectorType & stepInMM, + const JFPType::PointType & itkNotUsed(source), + const JFPType::VectorType & itkNotUsed(sourceToPixel), + const JFPType::PointType & itkNotUsed(nearestPoint), + const JFPType::PointType & itkNotUsed(farthestPoint)) { + JFPType::OutputPixelType tmp = static_cast(input); + if (tmp < rayCastValue) + { + tmp = rayCastValue; + } + output = tmp * stepInMM.GetNorm(); + }; + jfp->SetProjectedValueAccumulation(projAccumFunc); auto geometry = rtk::ThreeDCircularProjectionGeometry::New(); geometry->AddProjection(700, 800, 0); - mipfp->SetGeometry(geometry); - mipfp->Update(); + jfp->SetGeometry(geometry); + TRY_AND_EXIT_ON_ITK_EXCEPTION(jfp->Update()); - itk::ImageRegionConstIterator inputIt(mipfp->GetOutput(), mipfp->GetOutput()->GetRequestedRegion()); + itk::ImageRegionConstIterator inputIt(jfp->GetOutput(), jfp->GetOutput()->GetRequestedRegion()); inputIt.GoToBegin(); @@ -84,5 +121,41 @@ rtkmaximumintensityprojectiontest(int, char *[]) return EXIT_FAILURE; } std::cout << "\n\nTest PASSED! " << std::endl; + + // Test-2 + // Maximum Intensity Projection (MIP) filter, derived from Joseph Forward Projection filter + std::cout << "\n\n****** Case 2: MaximumIntensityProjection (MIP) filter ******" << std::endl; + using MIPType = rtk::MaximumIntensityProjectionImageFilter; + MIPType::Pointer mipfp = MIPType::New(); + mipfp->SetInput(projInput->GetOutput()); + mipfp->SetInput(1, volInput->GetOutput()); + + mipfp->SetGeometry(geometry); + TRY_AND_EXIT_ON_ITK_EXCEPTION(mipfp->Update()); + + itk::ImageRegionConstIterator inputIt2(mipfp->GetOutput(), mipfp->GetOutput()->GetRequestedRegion()); + + inputIt2.GoToBegin(); + + res = false; + while (!inputIt2.IsAtEnd()) + { + OutputPixelType pixel = inputIt2.Get(); + if (pixel < 4. || pixel > 4.25) + { + res = true; + } + ++inputIt2; + } + if (res) + { + return EXIT_FAILURE; + } + std::cout << "\n\nTest PASSED! " << std::endl; + + std::cout << "\n\n****** Compare two resulted images ******" << std::endl; + CheckImageQuality(jfp->GetOutput(), mipfp->GetOutput(), 0.001, 0.000001, 1.E+19); + std::cout << "\n\nTest PASSED! " << std::endl; + return EXIT_SUCCESS; } diff --git a/wrapping/rtkJosephForwardAttenuatedProjectionImageFilter.wrap b/wrapping/rtkJosephForwardAttenuatedProjectionImageFilter.wrap index 89457c8d3..df8e6dfec 100644 --- a/wrapping/rtkJosephForwardAttenuatedProjectionImageFilter.wrap +++ b/wrapping/rtkJosephForwardAttenuatedProjectionImageFilter.wrap @@ -1,28 +1,3 @@ -set(WRAPPER_AUTO_INCLUDE_HEADERS OFF) -itk_wrap_named_class("rtk::Functor::InterpolationWeightMultiplicationAttenuated" "rtkFunctorInterpolationWeightMultiplicationAttenuatedBackProjectionAttenuated") - foreach(t ${WRAP_ITK_REAL}) - itk_wrap_template("${ITKM_${t}}${ITKM_${t}}${ITKM_${t}}" "${ITKT_${t}}, ${ITKT_${t}}, ${ITKT_${t}}") - endforeach() -itk_end_wrap_class() -itk_wrap_named_class("rtk::Functor::ProjectedValueAccumulationAttenuated" "rtkProjectedValueAccumulationAttenuated") - foreach(t ${WRAP_ITK_REAL}) - itk_wrap_template("${ITKM_${t}}${ITKM_${t}}" "${ITKT_${t}}, ${ITKT_${t}}") - endforeach() -itk_end_wrap_class() -itk_wrap_named_class("rtk::Functor::ComputeAttenuationCorrection" "rtkComputeAttenuationCorrection") - foreach(t ${WRAP_ITK_REAL}) - itk_wrap_template("${ITKM_${t}}${ITKM_${t}}" "${ITKT_${t}}, ${ITKT_${t}}") - endforeach() -itk_end_wrap_class() -set(WRAPPER_AUTO_INCLUDE_HEADERS ON) - -itk_wrap_class("rtk::JosephForwardProjectionImageFilter" POINTER) - foreach(t ${WRAP_ITK_REAL}) - itk_wrap_template("I${ITKM_${t}}3I${ITKM_${t}}3SWM${ITKM_${t}}D${ITKM_${t}}IPC" - "itk::Image<${ITKT_${t}}, 3>, itk::Image< ${ITKT_${t}}, 3>, rtk::Functor::InterpolationWeightMultiplicationAttenuated<${ITKT_${t}}, ${ITKT_${t}}, ${ITKT_${t}}>, rtk::Functor::ProjectedValueAccumulationAttenuated<${ITKT_${t}}, ${ITKT_${t}}>, rtk::Functor::ComputeAttenuationCorrection<${ITKT_${t}}, ${ITKT_${t}}>") - endforeach() -itk_end_wrap_class() - itk_wrap_class("rtk::JosephForwardAttenuatedProjectionImageFilter" POINTER) foreach(t ${WRAP_ITK_REAL}) itk_wrap_template("I${ITKM_${t}}3I${ITKM_${t}}3" "itk::Image<${ITKT_${t}}, 3>, itk::Image<${ITKT_${t}}, 3>") diff --git a/wrapping/rtkMaximumIntensityProjectionImageFilter.wrap b/wrapping/rtkMaximumIntensityProjectionImageFilter.wrap index b4b1acce6..93b4d4317 100644 --- a/wrapping/rtkMaximumIntensityProjectionImageFilter.wrap +++ b/wrapping/rtkMaximumIntensityProjectionImageFilter.wrap @@ -1,29 +1,4 @@ -set(WRAPPER_AUTO_INCLUDE_HEADERS OFF) -itk_wrap_named_class("rtk::Functor::InterpolationWeightMultiplication" "rtkInterpolationWeightMultiplication") - foreach(t ${WRAP_ITK_REAL}) - itk_wrap_template("${ITKM_${t}}${ITKM_${t}}" "${ITKT_${t}}, ${ITKT_${t}}") - endforeach() -itk_end_wrap_class() -itk_wrap_named_class("rtk::Functor::MaximumIntensityAlongRay" "rtkMaximumIntensityAlongRay") - foreach(t ${WRAP_ITK_REAL}) - itk_wrap_template("${ITKM_${t}}${ITKM_${t}}" "${ITKT_${t}}, ${ITKT_${t}}") - endforeach() -itk_end_wrap_class() -itk_wrap_named_class("rtk::Functor::MaximumIntensityProjectedValueAccumulation" "rtkMaximumIntensityProjectedValueAccumulation") - foreach(t ${WRAP_ITK_REAL}) - itk_wrap_template("${ITKM_${t}}${ITKM_${t}}" "${ITKT_${t}}, ${ITKT_${t}}") - endforeach() -itk_end_wrap_class() -set(WRAPPER_AUTO_INCLUDE_HEADERS ON) - -itk_wrap_class("rtk::JosephForwardProjectionImageFilter" POINTER) - foreach(t ${WRAP_ITK_REAL}) - itk_wrap_template("I${ITKM_${t}}3I${ITKM_${t}}3FWMI${ITKM_${t}}3I${ITKM_${t}}3FVAI${ITKM_${t}}3I${ITKM_${t}}3FIAI${ITKM_${t}}3I${ITKM_${t}}3" - "itk::Image<${ITKT_${t}}, 3>, itk::Image< ${ITKT_${t}}, 3>, rtk::Functor::InterpolationWeightMultiplication<${ITKT_${t}}, ${ITKT_${t}}>, rtk::Functor::MaximumIntensityProjectedValueAccumulation<${ITKT_${t}}, ${ITKT_${t}}>, rtk::Functor::MaximumIntensityAlongRay<${ITKT_${t}}, ${ITKT_${t}}>") - endforeach() -itk_end_wrap_class() - -itk_wrap_class("rtk::MaximumIntensityProjectionImageFilter" POINTER) +itk_wrap_class("rtk::MaximumIntensityProjectionImageFilter" POINTER) foreach(t ${WRAP_ITK_REAL}) itk_wrap_template("I${ITKM_${t}}3I${ITKM_${t}}3" "itk::Image<${ITKT_${t}}, 3>, itk::Image<${ITKT_${t}}, 3>") endforeach() From 3b94ce7a16307891189b4c1953a8dd4126620fda Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Wed, 22 Apr 2026 16:40:45 +0200 Subject: [PATCH 2/2] ENH: Use lambda functions instead of functor classes in rtkJosephBackProjectionFilter --- .../rtkprojectionmatrix.cxx | 28 ++- ...osephBackAttenuatedProjectionImageFilter.h | 182 +------------- ...ephBackAttenuatedProjectionImageFilter.hxx | 116 +++++---- include/rtkJosephBackProjectionImageFilter.h | 233 +++++------------- .../rtkJosephBackProjectionImageFilter.hxx | 198 +++++++-------- ...phBackAttenuatedProjectionImageFilter.wrap | 25 -- .../rtkJosephBackProjectionImageFilter.wrap | 8 - 7 files changed, 231 insertions(+), 559 deletions(-) diff --git a/applications/rtkprojectionmatrix/rtkprojectionmatrix.cxx b/applications/rtkprojectionmatrix/rtkprojectionmatrix.cxx index 8ee38c154..19cd47517 100644 --- a/applications/rtkprojectionmatrix/rtkprojectionmatrix.cxx +++ b/applications/rtkprojectionmatrix/rtkprojectionmatrix.cxx @@ -159,19 +159,28 @@ main(int argc, char * argv[]) if (args_info.verbose_flag) std::cout << "Backprojecting volume and recording matrix values..." << std::endl; - auto backProjection = rtk::JosephBackProjectionImageFilter< - OutputImageType, - OutputImageType, - rtk::Functor::InterpolationWeightMultiplicationBackProjection, - rtk::Functor::StoreSparseMatrixSplatWeightMultiplication>::New(); + auto backProjection = rtk::JosephBackProjectionImageFilter::New(); backProjection->SetInput(constantImageSource->GetOutput()); backProjection->SetInput(1, reader->GetOutput()); backProjection->SetGeometry(geometry); - backProjection->GetSplatWeightMultiplication().SetProjectionsBuffer(reader->GetOutput()->GetBufferPointer()); - backProjection->GetSplatWeightMultiplication().SetVolumeBuffer(constantImageSource->GetOutput()->GetBufferPointer()); - backProjection->GetSplatWeightMultiplication().GetVnlSparseMatrix().resize( + // Create and configure a StoreSparseMatrixSplatWeightMultiplication instance + using SplatFunctorType = + rtk::Functor::StoreSparseMatrixSplatWeightMultiplication; + auto splatFunctor = std::make_shared(); + splatFunctor->SetProjectionsBuffer(reader->GetOutput()->GetBufferPointer()); + splatFunctor->SetVolumeBuffer(constantImageSource->GetOutput()->GetBufferPointer()); + splatFunctor->GetVnlSparseMatrix().resize( reader->GetOutput()->GetLargestPossibleRegion().GetNumberOfPixels(), constantImageSource->GetOutput()->GetLargestPossibleRegion().GetNumberOfPixels()); + + using WeightCoordinateType = typename itk::PixelTraits::ValueType; + backProjection->SetSplatWeightMultiplication([splatFunctor](const OutputPixelType & rayValue, + OutputPixelType & output, + const double stepLengthInVoxel, + const double voxelSize, + const WeightCoordinateType weight) { + (*splatFunctor)(rayValue, output, stepLengthInVoxel, voxelSize, weight); + }); TRY_AND_EXIT_ON_ITK_EXCEPTION(backProjection->Update()) // Write matrix to disk @@ -183,8 +192,7 @@ main(int argc, char * argv[]) std::cerr << "Failed to open " << args_info.output_arg << std::endl; return EXIT_FAILURE; } - rtk::MatlabSparseMatrix matlabSparseMatrix(backProjection->GetSplatWeightMultiplication().GetVnlSparseMatrix(), - backProjection->GetOutput()); + rtk::MatlabSparseMatrix matlabSparseMatrix(splatFunctor->GetVnlSparseMatrix(), backProjection->GetOutput()); ofs << matlabSparseMatrix; ofs.close(); return EXIT_SUCCESS; diff --git a/include/rtkJosephBackAttenuatedProjectionImageFilter.h b/include/rtkJosephBackAttenuatedProjectionImageFilter.h index c98d51da1..43f517eaa 100644 --- a/include/rtkJosephBackAttenuatedProjectionImageFilter.h +++ b/include/rtkJosephBackAttenuatedProjectionImageFilter.h @@ -25,160 +25,7 @@ namespace rtk { -namespace Functor -{ -/** \class InterpolationWeightMultiplicationAttenuatedBackProjection - * \brief Function to multiply the interpolation weights with the projected - * volume values and attenuation map. - * - * \author Antoine Robert - * - * \ingroup RTK Functions - */ -template -class ITK_TEMPLATE_EXPORT InterpolationWeightMultiplicationAttenuatedBackProjection -{ -public: - InterpolationWeightMultiplicationAttenuatedBackProjection() { m_AttenuationPixel = 0; } - - ~InterpolationWeightMultiplicationAttenuatedBackProjection() = default; - bool - operator!=(const InterpolationWeightMultiplicationAttenuatedBackProjection &) const - { - return false; - } - - bool - operator==(const InterpolationWeightMultiplicationAttenuatedBackProjection & other) const - { - return !(*this != other); - } - - TOutput - operator()(const double stepLengthInVoxel, const TCoordinateType weight, const TInput * p, const int i) - { - const double w = weight * stepLengthInVoxel; - - m_AttenuationPixel += w * (p + m_AttenuationMinusEmissionMapsPtrDiff)[i]; - return w * (p + m_AttenuationMinusEmissionMapsPtrDiff)[i]; - } - void - SetAttenuationMinusEmissionMapsPtrDiff(std::ptrdiff_t pd) - { - m_AttenuationMinusEmissionMapsPtrDiff = pd; - } - TOutput * - GetAttenuationPixel() - { - return &m_AttenuationPixel; - } - -private: - std::ptrdiff_t m_AttenuationMinusEmissionMapsPtrDiff{}; - TInput m_AttenuationPixel; -}; - -/** \class ComputeAttenuationCorrectionBackProjection - * \brief Function to compute the attenuation correction on the projection. - * - * \author Antoine Robert - * - * \ingroup RTK Functions - */ -template -class ITK_TEMPLATE_EXPORT ComputeAttenuationCorrectionBackProjection -{ -public: - using VectorType = itk::Vector; - - ComputeAttenuationCorrectionBackProjection() { m_Ex1 = 1; } - - ~ComputeAttenuationCorrectionBackProjection() = default; - bool - operator!=(const ComputeAttenuationCorrectionBackProjection &) const - { - return false; - } - - bool - operator==(const ComputeAttenuationCorrectionBackProjection & other) const - { - return !(*this != other); - } - - TOutput - operator()(const TInput rayValue, const TInput attenuationRay, const VectorType & stepInMM, bool & isNewRay) - { - if (isNewRay) - { - m_Ex1 = 1; - isNewRay = false; - } - TInput ex2 = exp(-attenuationRay * stepInMM.GetNorm()); - TInput wf; - if (*m_AttenuationPixel > 0) - { - wf = (m_Ex1 - ex2) / *m_AttenuationPixel; - } - else - { - wf = m_Ex1 * stepInMM.GetNorm(); - } - - m_Ex1 = ex2; - *m_AttenuationPixel = 0; - return wf * rayValue; - } - - void - SetAttenuationPixel(TInput * attenuationPixel) - { - m_AttenuationPixel = attenuationPixel; - } - -private: - TInput m_Ex1; - TInput * m_AttenuationPixel; -}; - -/** \class SplatWeightMultiplicationAttenuated - * \brief Function to multiply the interpolation weights with the projection - * values. - * - * \author Cyril Mory - * - * \ingroup RTK Functions - */ -template -class ITK_TEMPLATE_EXPORT SplatWeightMultiplicationAttenuated -{ -public: - SplatWeightMultiplicationAttenuated() = default; - ~SplatWeightMultiplicationAttenuated() = default; - bool - operator!=(const SplatWeightMultiplicationAttenuated &) const - { - return false; - } - - bool - operator==(const SplatWeightMultiplicationAttenuated & other) const - { - return !(*this != other); - } - - void - operator()(const TInput & rayValue, - TOutput & output, - const double stepLengthInVoxel, - const double itkNotUsed(voxelSize), - const TCoordinateType weight) const - { - output += rayValue * weight * stepLengthInVoxel; - } -}; -} // end namespace Functor /** \class JosephBackAttenuatedProjectionImageFilter * \brief Attenuated Joseph back projection. @@ -194,37 +41,21 @@ class ITK_TEMPLATE_EXPORT SplatWeightMultiplicationAttenuated * \ingroup RTK Projector */ -template < - class TInputImage, - class TOutputImage, - class TInterpolationWeightMultiplication = Functor::InterpolationWeightMultiplicationAttenuatedBackProjection< - typename TInputImage::PixelType, - typename itk::PixelTraits::ValueType>, - class TSplatWeightMultiplication = Functor:: - SplatWeightMultiplicationAttenuated, - class TSumAlongRay = Functor::ComputeAttenuationCorrectionBackProjection> +template class ITK_TEMPLATE_EXPORT JosephBackAttenuatedProjectionImageFilter - : public JosephBackProjectionImageFilter + : public JosephBackProjectionImageFilter { public: ITK_DISALLOW_COPY_AND_MOVE(JosephBackAttenuatedProjectionImageFilter); /** Standard class type alias. */ using Self = JosephBackAttenuatedProjectionImageFilter; - using Superclass = JosephBackProjectionImageFilter; + using Superclass = JosephBackProjectionImageFilter; using Pointer = itk::SmartPointer; using ConstPointer = itk::SmartPointer; using InputPixelType = typename TInputImage::PixelType; using OutputPixelType = typename TOutputImage::PixelType; + using WeightCoordinateType = typename itk::PixelTraits::ValueType; using OutputImageRegionType = typename TOutputImage::RegionType; using CoordinateType = double; using VectorType = itk::Vector; @@ -244,6 +75,11 @@ class ITK_TEMPLATE_EXPORT JosephBackAttenuatedProjectionImageFilter JosephBackAttenuatedProjectionImageFilter(); ~JosephBackAttenuatedProjectionImageFilter() override = default; + std::ptrdiff_t m_AttenuationMinusEmissionMapsPtrDiff{ 0 }; + InputPixelType m_AttenuationPixel{ static_cast(0) }; + InputPixelType m_Ex1{ static_cast(1) }; + const InputPixelType * m_AttenuationMapBuffer{ nullptr }; + /** Apply changes to the input image requested region. */ void GenerateInputRequestedRegion() override; diff --git a/include/rtkJosephBackAttenuatedProjectionImageFilter.hxx b/include/rtkJosephBackAttenuatedProjectionImageFilter.hxx index 2d683e2e4..6dbebe53d 100644 --- a/include/rtkJosephBackAttenuatedProjectionImageFilter.hxx +++ b/include/rtkJosephBackAttenuatedProjectionImageFilter.hxx @@ -29,33 +29,67 @@ namespace rtk { -template -JosephBackAttenuatedProjectionImageFilter::JosephBackAttenuatedProjectionImageFilter() +template +JosephBackAttenuatedProjectionImageFilter::JosephBackAttenuatedProjectionImageFilter() + : JosephBackProjectionImageFilter() { this->m_InferiorClip = 0.; this->m_SuperiorClip = 1.; this->SetNumberOfRequiredInputs(3); + + /** \brief Interpolation: accumulate attenuation pixel along the ray */ + auto interpolationWeightMultiplicationAttenuatedFunc = [this](double stepLengthInVoxel, + const WeightCoordinateType weight, + const InputPixelType * p, + const int i) -> OutputPixelType { + const double w = weight * stepLengthInVoxel; + this->m_AttenuationPixel += w * (p + this->m_AttenuationMinusEmissionMapsPtrDiff)[i]; + return static_cast(w * (p + this->m_AttenuationMinusEmissionMapsPtrDiff)[i]); + }; + this->SetInterpolationWeightMultiplication(interpolationWeightMultiplicationAttenuatedFunc); + + /** \brief Sum along ray: compute attenuation correction and reset accumulator */ + auto computeAttenuationCorrectionFunc = [this](const InputPixelType rayValue, + const InputPixelType attenuationRay, + const VectorType & stepInMM, + bool & isNewRay) -> OutputPixelType { + if (isNewRay) + { + this->m_Ex1 = static_cast(1); + isNewRay = false; + } + InputPixelType ex2 = static_cast(exp(-attenuationRay * stepInMM.GetNorm())); + InputPixelType wf; + if (this->m_AttenuationPixel > 0) + { + wf = (this->m_Ex1 - ex2) / this->m_AttenuationPixel; + } + else + { + wf = this->m_Ex1 * stepInMM.GetNorm(); + } + this->m_Ex1 = ex2; + this->m_AttenuationPixel = static_cast(0); + return static_cast(wf * rayValue); + }; + this->SetSumAlongRay(computeAttenuationCorrectionFunc); + + /** \brief Splat multiplication: simple accumulation */ + auto splatWeightMultiplicationAttenuatedFunc = [](const InputPixelType & rayValue, + OutputPixelType & output, + const double stepLengthInVoxel, + const double itkNotUsed(voxelSize), + const WeightCoordinateType weight) { + output += static_cast(rayValue * weight * stepLengthInVoxel); + }; + this->SetSplatWeightMultiplication(splatWeightMultiplicationAttenuatedFunc); + + this->DynamicMultiThreadingOff(); } -template +template void -JosephBackAttenuatedProjectionImageFilter::GenerateInputRequestedRegion() +JosephBackAttenuatedProjectionImageFilter::GenerateInputRequestedRegion() { // Input 2 is the attenuation map relative to the volume typename Superclass::InputImagePointer inputPtr2 = const_cast(this->GetInput(2)); @@ -64,17 +98,9 @@ JosephBackAttenuatedProjectionImageFilter +template void -JosephBackAttenuatedProjectionImageFilter::VerifyInputInformation() const +JosephBackAttenuatedProjectionImageFilter::VerifyInputInformation() const { using ImageBaseType = const itk::ImageBase; @@ -154,38 +180,22 @@ JosephBackAttenuatedProjectionImageFilter +template void -JosephBackAttenuatedProjectionImageFilter::Init() +JosephBackAttenuatedProjectionImageFilter::Init() { if (!this->GetInput(2)) { itkExceptionMacro("Attenuation map (input 2) must be set before running the attenuated back projector."); } - this->m_InterpolationWeightMultiplication.SetAttenuationMinusEmissionMapsPtrDiff( - this->GetInput(2)->GetBufferPointer() - this->GetInput(0)->GetBufferPointer()); - this->m_SumAlongRay.SetAttenuationPixel(this->m_InterpolationWeightMultiplication.GetAttenuationPixel()); + this->m_AttenuationMinusEmissionMapsPtrDiff = + this->GetInput(2)->GetBufferPointer() - this->GetInput(0)->GetBufferPointer(); + this->m_AttenuationMapBuffer = this->GetInput(2)->GetBufferPointer(); } -template +template void -JosephBackAttenuatedProjectionImageFilter::GenerateData() +JosephBackAttenuatedProjectionImageFilter::GenerateData() { Init(); Superclass::GenerateData(); diff --git a/include/rtkJosephBackProjectionImageFilter.h b/include/rtkJosephBackProjectionImageFilter.h index bdcd42d14..aa9e78939 100644 --- a/include/rtkJosephBackProjectionImageFilter.h +++ b/include/rtkJosephBackProjectionImageFilter.h @@ -24,120 +24,10 @@ #include "rtkThreeDCircularProjectionGeometry.h" #include +#include namespace rtk { -namespace Functor -{ -/** \class InterpolationWeightMultiplicationBackProjection - * \brief Function to multiply the interpolation weights with the projected - * volume values. - * - * \author Simon Rit - * - * \ingroup RTK Functions - */ -template -class ITK_TEMPLATE_EXPORT InterpolationWeightMultiplicationBackProjection -{ -public: - InterpolationWeightMultiplicationBackProjection() = default; - ~InterpolationWeightMultiplicationBackProjection() = default; - bool - operator!=(const InterpolationWeightMultiplicationBackProjection &) const - { - return false; - } - bool - operator==(const InterpolationWeightMultiplicationBackProjection & other) const - { - return !(*this != other); - } - - TOutput - operator()(const double itkNotUsed(stepLengthInVoxel), - const TCoordinateType itkNotUsed(weight), - const TInput * itkNotUsed(p), - const int itkNotUsed(i)) const - { - return {}; - } -}; - -/** \class SumAlongRay - * \brief Function to compute the attenuation correction on the projection. - * - * \author Antoine Robert - * - * \ingroup RTK Functions - */ -template -class ITK_TEMPLATE_EXPORT ValueAlongRay -{ -public: - using VectorType = itk::Vector; - - ValueAlongRay() = default; - ~ValueAlongRay() = default; - bool - operator!=(const ValueAlongRay &) const - { - return false; - } - bool - operator==(const ValueAlongRay & other) const - { - return !(*this != other); - } - - const TOutput & - operator()(const TInput & rayValue, - const TInput itkNotUsed(attenuationRay), - const VectorType & itkNotUsed(stepInMM), - bool itkNotUsed(isEndRay)) const - { - return rayValue; - } -}; -/** \class SplatWeightMultiplication - * \brief Function to multiply the interpolation weights with the projection - * values. - * - * \author Cyril Mory - * - * \ingroup RTK Functions - */ -template -class ITK_TEMPLATE_EXPORT SplatWeightMultiplication -{ -public: - SplatWeightMultiplication() = default; - ~SplatWeightMultiplication() = default; - bool - operator!=(const SplatWeightMultiplication &) const - { - return false; - } - bool - operator==(const SplatWeightMultiplication & other) const - { - return !(*this != other); - } - - void - operator()(const TInput & rayValue, - TOutput & output, - const double stepLengthInVoxel, - const double voxelSize, - const TCoordinateType weight) const - { - output += rayValue * weight * voxelSize * stepLengthInVoxel; - } -}; - -} // end namespace Functor - - /** \class JosephBackProjectionImageFilter * \brief Joseph back projection. * @@ -151,16 +41,7 @@ class ITK_TEMPLATE_EXPORT SplatWeightMultiplication * * \ingroup RTK Projector */ - -template < - class TInputImage, - class TOutputImage, - class TInterpolationWeightMultiplication = Functor::InterpolationWeightMultiplicationBackProjection< - typename TInputImage::PixelType, - typename itk::PixelTraits::ValueType>, - class TSplatWeightMultiplication = - Functor::SplatWeightMultiplication, - class TSumAlongRay = Functor::ValueAlongRay> +template class ITK_TEMPLATE_EXPORT JosephBackProjectionImageFilter : public BackProjectionImageFilter { public: @@ -173,84 +54,86 @@ class ITK_TEMPLATE_EXPORT JosephBackProjectionImageFilter : public BackProjectio using ConstPointer = itk::SmartPointer; using InputPixelType = typename TInputImage::PixelType; using OutputPixelType = typename TOutputImage::PixelType; + using WeightCoordinateType = typename itk::PixelTraits::ValueType; + using OutputImageRegionType = typename TOutputImage::RegionType; using CoordinateType = double; using VectorType = itk::Vector; using GeometryType = rtk::ThreeDCircularProjectionGeometry; using GeometryPointer = typename GeometryType::Pointer; - /** Method for creation through the object factory. */ - itkNewMacro(Self); + /** \brief Function to multiply the interpolation weights with the sampled + * volume values. + */ + using InterpolationWeightMultiplicationFunc = + std::function; + + /** \brief Function to compute the value to backproject from a projection + * sample. + */ + using SumAlongRayFunc = + std::function; - /** Run-time type information (and related methods). */ + /** \brief Function to apply the splat weight multiplication and accumulate + * the contribution into the output voxel. + */ + using SplatWeightMultiplicationFunc = std::function< + void(const InputPixelType &, OutputPixelType &, const double, const double, const WeightCoordinateType)>; + + itkNewMacro(Self); itkOverrideGetNameOfClassMacro(JosephBackProjectionImageFilter); - /** Get/Set the functor that is used to multiply each interpolation value with a volume value */ - TInterpolationWeightMultiplication & + InterpolationWeightMultiplicationFunc & GetInterpolationWeightMultiplication() { return m_InterpolationWeightMultiplication; } - const TInterpolationWeightMultiplication & + const InterpolationWeightMultiplicationFunc & GetInterpolationWeightMultiplication() const { return m_InterpolationWeightMultiplication; } void - SetInterpolationWeightMultiplication(const TInterpolationWeightMultiplication & _arg) + SetInterpolationWeightMultiplication(const InterpolationWeightMultiplicationFunc & _arg) { - if (m_InterpolationWeightMultiplication != _arg) - { - m_InterpolationWeightMultiplication = _arg; - this->Modified(); - } + m_InterpolationWeightMultiplication = _arg; + this->Modified(); } - /** Get/Set the functor that is used to multiply each interpolation value with a volume value */ - TSplatWeightMultiplication & + SplatWeightMultiplicationFunc & GetSplatWeightMultiplication() { return m_SplatWeightMultiplication; } - const TSplatWeightMultiplication & + const SplatWeightMultiplicationFunc & GetSplatWeightMultiplication() const { return m_SplatWeightMultiplication; } void - SetSplatWeightMultiplication(const TSplatWeightMultiplication & _arg) + SetSplatWeightMultiplication(const SplatWeightMultiplicationFunc & _arg) { - if (m_SplatWeightMultiplication != _arg) - { - m_SplatWeightMultiplication = _arg; - this->Modified(); - } + m_SplatWeightMultiplication = _arg; + this->Modified(); } - /** Get/Set the functor that is used to compute the sum along the ray*/ - TSumAlongRay & + SumAlongRayFunc & GetSumAlongRay() { return m_SumAlongRay; } - const TSumAlongRay & + const SumAlongRayFunc & GetSumAlongRay() const { return m_SumAlongRay; } void - SetSumAlongRay(const TSumAlongRay & _arg) + SetSumAlongRay(const SumAlongRayFunc & _arg) { - if (m_SumAlongRay != _arg) - { - m_SumAlongRay = _arg; - this->Modified(); - } + m_SumAlongRay = _arg; + this->Modified(); } - /** Each ray is clipped from source+m_InferiorClip*(pixel-source) to - ** source+m_SuperiorClip*(pixel-source) with m_InferiorClip and - ** m_SuperiorClip equal 0 and 1 by default. */ itkGetMacro(InferiorClip, double); itkSetMacro(InferiorClip, double); itkGetMacro(SuperiorClip, double); @@ -263,8 +146,6 @@ class ITK_TEMPLATE_EXPORT JosephBackProjectionImageFilter : public BackProjectio void GenerateData() override; - /** The two inputs should not be in the same space so there is nothing - * to verify. */ void VerifyInputInformation() const override {} @@ -277,8 +158,8 @@ class ITK_TEMPLATE_EXPORT JosephBackProjectionImageFilter : public BackProjectio OutputPixelType * pxsyi, OutputPixelType * pxiys, OutputPixelType * pxsys, - double x, - double y, + CoordinateType x, + CoordinateType y, int ox, int oy); @@ -290,8 +171,8 @@ class ITK_TEMPLATE_EXPORT JosephBackProjectionImageFilter : public BackProjectio OutputPixelType * pxsyi, OutputPixelType * pxiys, OutputPixelType * pxsys, - double x, - double y, + CoordinateType x, + CoordinateType y, int ox, int oy, CoordinateType minx, @@ -305,8 +186,8 @@ class ITK_TEMPLATE_EXPORT JosephBackProjectionImageFilter : public BackProjectio const InputPixelType * pxsyi, const InputPixelType * pxiys, const InputPixelType * pxsys, - double x, - double y, + CoordinateType x, + CoordinateType y, int ox, int oy); @@ -316,21 +197,23 @@ class ITK_TEMPLATE_EXPORT JosephBackProjectionImageFilter : public BackProjectio const InputPixelType * pxsyi, const InputPixelType * pxiys, const InputPixelType * pxsys, - double x, - double y, + CoordinateType x, + CoordinateType y, int ox, int oy, - double minx, - double miny, - double maxx, - double maxy); - - /** Functor */ - TSplatWeightMultiplication m_SplatWeightMultiplication; - TInterpolationWeightMultiplication m_InterpolationWeightMultiplication; - TSumAlongRay m_SumAlongRay; - double m_InferiorClip{ 0. }; - double m_SuperiorClip{ 1. }; + CoordinateType minx, + CoordinateType miny, + CoordinateType maxx, + CoordinateType maxy); + + // lambdas or std::functions + // MUST BE initiated in constructor, or custom functions + // must be set by means of Set methods in application + SplatWeightMultiplicationFunc m_SplatWeightMultiplication; + InterpolationWeightMultiplicationFunc m_InterpolationWeightMultiplication; + SumAlongRayFunc m_SumAlongRay; + double m_InferiorClip{ 0. }; + double m_SuperiorClip{ 1. }; }; } // end namespace rtk diff --git a/include/rtkJosephBackProjectionImageFilter.hxx b/include/rtkJosephBackProjectionImageFilter.hxx index 7673d6fab..cfab9661f 100644 --- a/include/rtkJosephBackProjectionImageFilter.hxx +++ b/include/rtkJosephBackProjectionImageFilter.hxx @@ -33,28 +33,35 @@ namespace rtk { -template -JosephBackProjectionImageFilter::JosephBackProjectionImageFilter() = default; - -template +template +JosephBackProjectionImageFilter::JosephBackProjectionImageFilter() +{ + this->m_InterpolationWeightMultiplication = [](double itkNotUsed(stepLengthInVoxel), + const WeightCoordinateType weight, + const InputPixelType * p, + int i) -> OutputPixelType { + return static_cast(weight * p[i]); + }; + + this->m_SplatWeightMultiplication = [](const InputPixelType & rayValue, + OutputPixelType & output, + const double stepLengthInVoxel, + const double voxelSize, + const WeightCoordinateType weight) { + output += static_cast(rayValue * weight * voxelSize * stepLengthInVoxel); + }; + + this->m_SumAlongRay = [](const InputPixelType & rayValue, + const InputPixelType itkNotUsed(attenuationRay), + const VectorType & itkNotUsed(stepInMM), + bool & itkNotUsed(isEndRay)) -> OutputPixelType { + return static_cast(rayValue); + }; +} + +template void -JosephBackProjectionImageFilter::GenerateData() +JosephBackProjectionImageFilter::GenerateData() { // Allocate the output image this->AllocateOutputs(); @@ -322,27 +329,19 @@ JosephBackProjectionImageFilter +template void -JosephBackProjectionImageFilter::BilinearSplat(const InputPixelType & rayValue, - const double stepLengthInVoxel, - const double voxelSize, - OutputPixelType * pxiyi, - OutputPixelType * pxsyi, - OutputPixelType * pxiys, - OutputPixelType * pxsys, - const double x, - const double y, - const int ox, - const int oy) +JosephBackProjectionImageFilter::BilinearSplat(const InputPixelType & rayValue, + const double stepLengthInVoxel, + const double voxelSize, + OutputPixelType * pxiyi, + OutputPixelType * pxsyi, + OutputPixelType * pxiys, + OutputPixelType * pxsys, + const CoordinateType x, + const CoordinateType y, + const int ox, + const int oy) { int ix = itk::Math::floor(x); int iy = itk::Math::floor(y); @@ -358,31 +357,23 @@ JosephBackProjectionImageFilter +template void -JosephBackProjectionImageFilter::BilinearSplatOnBorders(const InputPixelType & rayValue, - const double stepLengthInVoxel, - const double voxelSize, - OutputPixelType * pxiyi, - OutputPixelType * pxsyi, - OutputPixelType * pxiys, - OutputPixelType * pxsys, - const double x, - const double y, - const int ox, - const int oy, - const CoordinateType minx, - const CoordinateType miny, - const CoordinateType maxx, - const CoordinateType maxy) +JosephBackProjectionImageFilter::BilinearSplatOnBorders(const InputPixelType & rayValue, + const double stepLengthInVoxel, + const double voxelSize, + OutputPixelType * pxiyi, + OutputPixelType * pxsyi, + OutputPixelType * pxiys, + OutputPixelType * pxsys, + const CoordinateType x, + const CoordinateType y, + const int ox, + const int oy, + const CoordinateType minx, + const CoordinateType miny, + const CoordinateType maxx, + const CoordinateType maxy) { int ix = itk::Math::floor(x); int iy = itk::Math::floor(y); @@ -412,29 +403,17 @@ JosephBackProjectionImageFilter -typename JosephBackProjectionImageFilter::OutputPixelType -JosephBackProjectionImageFilter::BilinearInterpolation(const double stepLengthInVoxel, - const InputPixelType * pxiyi, - const InputPixelType * pxsyi, - const InputPixelType * pxiys, - const InputPixelType * pxsys, - const CoordinateType x, - const CoordinateType y, - const int ox, - const int oy) +template +typename JosephBackProjectionImageFilter::OutputPixelType +JosephBackProjectionImageFilter::BilinearInterpolation(const double stepLengthInVoxel, + const InputPixelType * pxiyi, + const InputPixelType * pxsyi, + const InputPixelType * pxiys, + const InputPixelType * pxsys, + const CoordinateType x, + const CoordinateType y, + const int ox, + const int oy) { int ix = itk::Math::floor(x); int iy = itk::Math::floor(y); @@ -449,33 +428,22 @@ JosephBackProjectionImageFilter -typename JosephBackProjectionImageFilter::OutputPixelType -JosephBackProjectionImageFilter::BilinearInterpolationOnBorders(const double stepLengthInVoxel, - const InputPixelType * pxiyi, - const InputPixelType * pxsyi, - const InputPixelType * pxiys, - const InputPixelType * pxsys, - const CoordinateType x, - const CoordinateType y, - const int ox, - const int oy, - const CoordinateType minx, - const CoordinateType miny, - const CoordinateType maxx, - const CoordinateType maxy) +template +typename JosephBackProjectionImageFilter::OutputPixelType +JosephBackProjectionImageFilter::BilinearInterpolationOnBorders( + const double stepLengthInVoxel, + const InputPixelType * pxiyi, + const InputPixelType * pxsyi, + const InputPixelType * pxiys, + const InputPixelType * pxsys, + const CoordinateType x, + const CoordinateType y, + const int ox, + const int oy, + const CoordinateType minx, + const CoordinateType miny, + const CoordinateType maxx, + const CoordinateType maxy) { int ix = itk::Math::floor(x); int iy = itk::Math::floor(y); diff --git a/wrapping/rtkJosephBackAttenuatedProjectionImageFilter.wrap b/wrapping/rtkJosephBackAttenuatedProjectionImageFilter.wrap index ce3f7da17..5556638a9 100644 --- a/wrapping/rtkJosephBackAttenuatedProjectionImageFilter.wrap +++ b/wrapping/rtkJosephBackAttenuatedProjectionImageFilter.wrap @@ -1,28 +1,3 @@ -set(WRAPPER_AUTO_INCLUDE_HEADERS OFF) -itk_wrap_named_class("rtk::Functor::InterpolationWeightMultiplicationAttenuatedBackProjection" "rtkFunctorInterpolationWeightMultiplicationAttenuatedBackProjection") - foreach(t ${WRAP_ITK_REAL}) - itk_wrap_template("${ITKM_${t}}${ITKM_${t}}${ITKM_${t}}" "${ITKT_${t}}, ${ITKT_${t}}, ${ITKT_${t}}") - endforeach() -itk_end_wrap_class() -itk_wrap_named_class("rtk::Functor::SplatWeightMultiplicationAttenuated" "rtkFunctorSplatWeightMultiplicationAttenuated") - foreach(t ${WRAP_ITK_REAL}) - itk_wrap_template("${ITKM_${t}}D${ITKM_${t}}" "${ITKT_${t}}, double, ${ITKT_${t}}") - endforeach() -itk_end_wrap_class() -itk_wrap_named_class("rtk::Functor::ComputeAttenuationCorrectionBackProjection" "rtkComputeAttenuationCorrectionBackProjection") - foreach(t ${WRAP_ITK_REAL}) - itk_wrap_template("${ITKM_${t}}${ITKM_${t}}" "${ITKT_${t}}, ${ITKT_${t}}") - endforeach() -itk_end_wrap_class() -set(WRAPPER_AUTO_INCLUDE_HEADERS ON) - -itk_wrap_class("rtk::JosephBackProjectionImageFilter" POINTER) - foreach(t ${WRAP_ITK_REAL}) - itk_wrap_template("I${ITKM_${t}}3I${ITKM_${t}}3SWM${ITKM_${t}}D${ITKM_${t}}IS" - "itk::Image<${ITKT_${t}}, 3>, itk::Image< ${ITKT_${t}}, 3>, rtk::Functor::InterpolationWeightMultiplicationAttenuatedBackProjection<${ITKT_${t}}, ${ITKT_${t}}, ${ITKT_${t}}>, rtk::Functor::SplatWeightMultiplicationAttenuated<${ITKT_${t}}, double, ${ITKT_${t}}>, rtk::Functor::ComputeAttenuationCorrectionBackProjection<${ITKT_${t}}, ${ITKT_${t}}>") - endforeach() -itk_end_wrap_class() - itk_wrap_class("rtk::JosephBackAttenuatedProjectionImageFilter" POINTER) foreach(t ${WRAP_ITK_REAL}) itk_wrap_template("I${ITKM_${t}}3I${ITKM_${t}}3SWM${ITKM_${t}}D${ITKM_${t}}" diff --git a/wrapping/rtkJosephBackProjectionImageFilter.wrap b/wrapping/rtkJosephBackProjectionImageFilter.wrap index d1d437fe2..ceb71265c 100644 --- a/wrapping/rtkJosephBackProjectionImageFilter.wrap +++ b/wrapping/rtkJosephBackProjectionImageFilter.wrap @@ -1,11 +1,3 @@ -set(WRAPPER_AUTO_INCLUDE_HEADERS OFF) -itk_wrap_named_class("rtk::Functor::SplatWeightMultiplication" "rtkFunctorSplatWeightMultiplication") - foreach(t ${WRAP_ITK_REAL}) - itk_wrap_template("${ITKM_${t}}D${ITKM_${t}}" "${ITKT_${t}}, double, ${ITKT_${t}}") - endforeach() -itk_end_wrap_class() -set(WRAPPER_AUTO_INCLUDE_HEADERS ON) - itk_wrap_class("rtk::JosephBackProjectionImageFilter" POINTER) foreach(t ${WRAP_ITK_REAL}) itk_wrap_template("I${ITKM_${t}}3I${ITKM_${t}}3SWM${ITKM_${t}}D${ITKM_${t}}"