diff --git a/applications/rtkforwardprojections/rtkforwardprojections.cxx b/applications/rtkforwardprojections/rtkforwardprojections.cxx index d34405306..590d46a84 100644 --- a/applications/rtkforwardprojections/rtkforwardprojections.cxx +++ b/applications/rtkforwardprojections/rtkforwardprojections.cxx @@ -27,6 +27,7 @@ #ifdef RTK_USE_CUDA # include "rtkCudaForwardProjectionImageFilter.h" #endif +#include "rtkAttenuatedToExponentialCorrectionImageFilter.h" #include #include @@ -81,6 +82,14 @@ main(int argc, char * argv[]) attenuationMap = itk::ReadImage(args_info.attenuationmap_arg); } + OutputImageType::Pointer KRegion; + if (args_info.kregion_given) + { + if (args_info.verbose_flag) + std::cout << "Reading K region " << args_info.kregion_arg << "..." << std::endl; + KRegion = itk::ReadImage(args_info.kregion_arg); + } + using ClipImageType = itk::Image; ClipImageType::Pointer inferiorClipImage, superiorClipImage; if (args_info.inferiorclipimage_given) @@ -109,6 +118,10 @@ main(int argc, char * argv[]) case (fp_arg_Joseph): forwardProjection = rtk::JosephForwardProjectionImageFilter::New(); break; + case (fp_arg_AttToExp): + forwardProjection = rtk::AttenuatedToExponentialCorrectionImageFilter::New(); + constantImageSource->SetConstant(1.); + break; case (fp_arg_JosephAttenuated): forwardProjection = rtk::JosephForwardAttenuatedProjectionImageFilter::New(); break; @@ -137,6 +150,8 @@ main(int argc, char * argv[]) forwardProjection->SetInput(1, inputVolume); if (args_info.attenuationmap_given) forwardProjection->SetInput(2, attenuationMap); + if (args_info.kregion_arg) + forwardProjection->SetInput(2, KRegion); if (args_info.inferiorclipimage_given) { if (args_info.fp_arg == fp_arg_Joseph) diff --git a/applications/rtkforwardprojections/rtkforwardprojections.ggo b/applications/rtkforwardprojections/rtkforwardprojections.ggo index d867b7690..84d3ad390 100644 --- a/applications/rtkforwardprojections/rtkforwardprojections.ggo +++ b/applications/rtkforwardprojections/rtkforwardprojections.ggo @@ -10,9 +10,10 @@ option "step" s "Step size along ray (for CudaRayCast only)" option "lowmem" l "Compute only one projection at a time" flag off section "Projectors" -option "fp" f "Forward projection method" values="Joseph","JosephAttenuated","CudaRayCast","Zeng","MIP" enum no default="Joseph" -option "attenuationmap" - "Attenuation map relative to the volume to perfom the attenuation correction" string no -option "sigmazero" - "PSF value at a distance of 0 meter of the detector" double no -option "alphapsf" - "Slope of the PSF against the detector distance" double no +option "fp" f "Forward projection method" values="Joseph","JosephAttenuated","CudaRayCast","Zeng","MIP","AttToExp" enum no default="Joseph" +option "attenuationmap" - "Attenuation map relative to the volume to perfom the attenuation correction" string no +option "sigmazero" - "PSF value at a distance of 0 meter of the detector" double no +option "alphapsf" - "Slope of the PSF against the detector distance" double no option "inferiorclipimage" - "Value of the inferior clip of the ray for each pixel of the projections (only with Joseph-based projector)" string no option "superiorclipimage" - "Value of the superior clip of the ray for each pixel of the projections (only with Joseph-based projector)" string no +option "kregion" - "K region file name for AttToExp, the correction from attenuated to exponential Radon transform" string no diff --git a/include/rtkAttenuatedToExponentialCorrectionImageFilter.h b/include/rtkAttenuatedToExponentialCorrectionImageFilter.h new file mode 100644 index 000000000..f671378fc --- /dev/null +++ b/include/rtkAttenuatedToExponentialCorrectionImageFilter.h @@ -0,0 +1,356 @@ +/*========================================================================= + * + * 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 rtkAttenuatedToExponentialCorrectionImageFilter_h +#define rtkAttenuatedToExponentialCorrectionImageFilter_h + +#include "rtkConfiguration.h" +#include "rtkForwardProjectionImageFilter.h" +#include "rtkJosephForwardProjectionImageFilter.h" +#include "rtkMacro.h" +#include +#include +#include + +namespace rtk +{ +namespace Functor +{ +/** \class InterpolationWeightMultiplicationAttToExp + * \brief Function to check if pixel along is in the K region and to record the + * step length. + * + * \author Simon Rit + * + * \ingroup RTK Functions + */ +template +class ITK_TEMPLATE_EXPORT InterpolationWeightMultiplicationAttToExp +{ +public: + InterpolationWeightMultiplicationAttToExp() = default; + + ~InterpolationWeightMultiplicationAttToExp() = default; + bool + operator!=(const InterpolationWeightMultiplicationAttToExp &) const + { + return false; + } + + bool + operator==(const InterpolationWeightMultiplicationAttToExp & other) const + { + return !(*this != other); + } + + inline TOutput + operator()(const ThreadIdType threadId, + const double stepLengthInVoxel, + const TCoordRepType weight, + const TInput * p, + const int i) + { + const TInput kRegion = (p + m_KRegionMinusAttenuationMapPtrDiff)[i]; + const int bNotInKRegion = (kRegion == itk::NumericTraits::ZeroValue()); + m_BeforeKRegion[threadId] = m_BeforeKRegion[threadId] * bNotInKRegion; + m_LatestStepLength[threadId] = stepLengthInVoxel; + return weight * p[i]; + } + + void + SetKRegionMinusAttenuationMapPtrDiff(std::ptrdiff_t pd) + { + m_KRegionMinusAttenuationMapPtrDiff = pd; + } + + int * + GetBeforeKRegion() + { + return m_BeforeKRegion; + } + + double * + GetLatestStepLength() + { + return m_LatestStepLength; + } + +private: + std::ptrdiff_t m_KRegionMinusAttenuationMapPtrDiff{}; + int m_BeforeKRegion[itk::ITK_MAX_THREADS]{ 1 }; + double m_LatestStepLength[itk::ITK_MAX_THREADS]{}; +}; + +/** \class SumAttenuationForCorrection + * \brief Function to calculate the correction factor from attenuation to + * exponential Radon transform. + * + * \author Simon Rit + * + * \ingroup RTK Functions + */ +template +class ITK_TEMPLATE_EXPORT SumAttenuationForCorrection +{ +public: + using VectorType = itk::Vector; + + SumAttenuationForCorrection() = default; + + ~SumAttenuationForCorrection() = default; + bool + operator!=(const SumAttenuationForCorrection &) const + { + return false; + } + + bool + operator==(const SumAttenuationForCorrection & other) const + { + return !(*this != other); + } + + inline void + operator()(const ThreadIdType threadId, + TOutput & sumValue, + const TInput volumeValue, + const VectorType & itkNotUsed(stepInMM)) + { + sumValue += static_cast(volumeValue) * m_BeforeKRegion[threadId]; + m_TraveledBeforeKRegion[threadId] += m_LatestStepLength[threadId] * m_BeforeKRegion[threadId]; + } + + double * + GetTraveledBeforeKRegion() + { + return m_TraveledBeforeKRegion; + } + + void + SetBeforeKRegion(int * pixelInKRegion) + { + m_BeforeKRegion = pixelInKRegion; + } + + void + SetLatestStepLength(double * latestStepLength) + { + m_LatestStepLength = latestStepLength; + } + +private: + int * m_BeforeKRegion{}; + double m_TraveledBeforeKRegion[itk::ITK_MAX_THREADS]{ 0. }; + double * m_LatestStepLength{}; +}; + +/** \class ProjectedValueAccumulationAttToExp + * \brief Function to calculate the correction factor from attenuation to + * exponential Radon transform using the integral of the attenuation before the + * K region and the length traveled in the attenuation image before the K + * region. + * + * \author Simon Rit + * + * \ingroup RTK Functions + */ +template +class ITK_TEMPLATE_EXPORT ProjectedValueAccumulationAttToExp +{ +public: + using VectorType = itk::Vector; + + ProjectedValueAccumulationAttToExp() = default; + ~ProjectedValueAccumulationAttToExp() = default; + bool + operator!=(const ProjectedValueAccumulationAttToExp &) const + { + return false; + } + + bool + operator==(const ProjectedValueAccumulationAttToExp & 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(pixel), + const VectorType & pixelToSource, + const VectorType & nearestPoint, + const VectorType & itkNotUsed(farthestPoint)) + { + // If we have not hit the K region, we set the correction factor to 0. as + // there shouldn't be any emission outside the K region. + if (m_BeforeKRegion[threadId]) + { + output = 0.; + m_TraveledBeforeKRegion[threadId] = 0.; + return; + } + + // Calculate tau, the distance from the entrance point of the K region to + // the origin. We assume that the detector is orthogonal to the + // pixelToSource line. + VectorType originToNearest = nearestPoint - m_Origin; + VectorType originToNearestInMM = itk::MakeVector( + originToNearest[0] * m_Spacing[0], originToNearest[1] * m_Spacing[1], originToNearest[2] * m_Spacing[2]); + VectorType nearestToKRegionInMM = m_TraveledBeforeKRegion[threadId] * stepInMM; + VectorType originToKRegionInMM = originToNearestInMM + nearestToKRegionInMM; + VectorType pixelToSourceInMM = itk::MakeVector( + pixelToSource[0] * m_Spacing[0], pixelToSource[1] * m_Spacing[1], pixelToSource[2] * m_Spacing[2]); + double tau = originToKRegionInMM * pixelToSourceInMM / -pixelToSourceInMM.GetNorm(); + output = input * std::exp(rayCastValue * stepInMM.GetNorm() + tau * m_Mu0); + + // Reinitialize for next ray + m_BeforeKRegion[threadId] = 1; + m_TraveledBeforeKRegion[threadId] = 0.; + } + + void + SetBeforeKRegion(int * pixelInKRegion) + { + m_BeforeKRegion = pixelInKRegion; + } + + void + SetTraveledBeforeKRegion(double * traveledBeforeKRegion) + { + m_TraveledBeforeKRegion = traveledBeforeKRegion; + } + + void + SetMu0(TOutput mu0) + { + m_Mu0 = mu0; + } + + void + SetOrigin(VectorType origin) + { + m_Origin = origin; + } + + void + SetSpacing(VectorType spacing) + { + m_Spacing = spacing; + } + +private: + int * m_BeforeKRegion{}; + double * m_TraveledBeforeKRegion{}; + TOutput m_Mu0{}; + VectorType m_Origin{}; + VectorType m_Spacing{}; +}; +} // end namespace Functor + +/** \class AttenuatedToExponentialCorrectionImageFilter + * \brief Converts input projections from the attenuation to exponential Radon transform + * + * The conversion is explained p47 of Natterer's book, "The mathematics of + * computerized tomography", 1986. The conversion assumes that the emission is + * contained in a K region (referred to as $\Omega$ in the book) of constant + * attenuation $\mu_0$. The projector therefore uses a mask as third input with + * the same voxel lattice as the attenuation map in the second input. $\mu_0$ is + * calculated as the average in the K region. + * + * \test TODO + * + * \author Simon Rit + * + * \ingroup RTK Projector + */ + +template < + class TInputImage, + class TOutputImage, + class TInterpolationWeightMultiplication = Functor::InterpolationWeightMultiplicationAttToExp< + typename TInputImage::PixelType, + typename itk::PixelTraits::ValueType>, + class TProjectedValueAccumulation = + Functor::ProjectedValueAccumulationAttToExp, + class TSumAlongRay = + Functor::SumAttenuationForCorrection> +class ITK_TEMPLATE_EXPORT AttenuatedToExponentialCorrectionImageFilter + : public JosephForwardProjectionImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(AttenuatedToExponentialCorrectionImageFilter); + + /** Standard class type alias. */ + using Self = AttenuatedToExponentialCorrectionImageFilter; + 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 CoordRepType = double; + using VectorType = itk::Vector; + + /** ImageDimension constants */ + static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ +#ifdef itkOverrideGetNameOfClassMacro + itkOverrideGetNameOfClassMacro(AttenuatedToExponentialCorrectionImageFilter); +#else + itkTypeMacro(AttenuatedToExponentialCorrectionImageFilter, JosephForwardProjectionImageFilter); +#endif + +protected: + AttenuatedToExponentialCorrectionImageFilter(); + ~AttenuatedToExponentialCorrectionImageFilter() override = default; + + /** Apply changes to the input image requested region. */ + void + GenerateInputRequestedRegion() override; + + void + BeforeThreadedGenerateData() override; + + /** Only the last two inputs should be in the same space so we need + * to overwrite the method. */ + void + VerifyInputInformation() const override; +}; +} // end namespace rtk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "rtkAttenuatedToExponentialCorrectionImageFilter.hxx" +#endif + +#endif diff --git a/include/rtkAttenuatedToExponentialCorrectionImageFilter.hxx b/include/rtkAttenuatedToExponentialCorrectionImageFilter.hxx new file mode 100644 index 000000000..d1732d375 --- /dev/null +++ b/include/rtkAttenuatedToExponentialCorrectionImageFilter.hxx @@ -0,0 +1,209 @@ +/*========================================================================= + * + * 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 rtkAttenuatedToExponentialCorrectionImageFilter_hxx +#define rtkAttenuatedToExponentialCorrectionImageFilter_hxx + + +#include "rtkHomogeneousMatrix.h" +#include "rtkBoxShape.h" +#include "rtkProjectionsRegionConstIteratorRayBased.h" + +#include +#include +#include +#include + +namespace rtk +{ +template +AttenuatedToExponentialCorrectionImageFilter::AttenuatedToExponentialCorrectionImageFilter() +{ + this->SetNumberOfRequiredInputs(3); + this->DynamicMultiThreadingOff(); +} + +template +void +AttenuatedToExponentialCorrectionImageFilter::GenerateInputRequestedRegion() +{ + Superclass::GenerateInputRequestedRegion(); + // Input 2 is the attenuation map relative to the volume + typename Superclass::InputImagePointer inputPtr2 = const_cast(this->GetInput(2)); + if (!inputPtr2) + return; + + typename TInputImage::RegionType reqRegion2 = inputPtr2->GetLargestPossibleRegion(); + inputPtr2->SetRequestedRegion(reqRegion2); +} + +template +void +AttenuatedToExponentialCorrectionImageFilter::VerifyInputInformation() const +{ + Superclass::VerifyInputInformation(); + using ImageBaseType = const itk::ImageBase; + + ImageBaseType * inputPtr1 = nullptr; + + itk::InputDataObjectConstIterator it(this); + for (; !it.IsAtEnd(); ++it) + { + // Check whether the output is an image of the appropriate + // dimension (use ProcessObject's version of the GetInput() + // method since it returns the input as a pointer to a + // DataObject as opposed to the subclass version which + // static_casts the input to an TInputImage). + if (it.GetName() == "_1") + { + inputPtr1 = dynamic_cast(it.GetInput()); + } + if (inputPtr1) + { + break; + } + } + + for (; !it.IsAtEnd(); ++it) + { + if (it.GetName() == "_2") + { + auto * inputPtrN = dynamic_cast(it.GetInput()); + // Physical space computation only matters if we're using two + // images, and not an image and a constant. + if (inputPtrN) + { + // check that the image occupy the same physical space, and that + // each index is at the same physical location + + // tolerance for origin and spacing depends on the size of pixel + // tolerance for directions a fraction of the unit cube. + const double coordinateTol = itk::Math::abs(Self::GetGlobalDefaultCoordinateTolerance() * + inputPtr1->GetSpacing()[0]); // use first dimension spacing + + if (!inputPtr1->GetOrigin().GetVnlVector().is_equal(inputPtrN->GetOrigin().GetVnlVector(), coordinateTol) || + !inputPtr1->GetSpacing().GetVnlVector().is_equal(inputPtrN->GetSpacing().GetVnlVector(), coordinateTol) || + !inputPtr1->GetDirection().GetVnlMatrix().as_ref().is_equal( + inputPtrN->GetDirection().GetVnlMatrix().as_ref(), Self::GetGlobalDefaultDirectionTolerance())) + { + std::ostringstream originString, spacingString, directionString; + if (!inputPtr1->GetOrigin().GetVnlVector().is_equal(inputPtrN->GetOrigin().GetVnlVector(), coordinateTol)) + { + originString.setf(std::ios::scientific); + originString.precision(7); + originString << "InputImage Origin: " << inputPtr1->GetOrigin() << ", InputImage" << it.GetName() + << " Origin: " << inputPtrN->GetOrigin() << std::endl; + originString << "\tTolerance: " << coordinateTol << std::endl; + } + if (!inputPtr1->GetSpacing().GetVnlVector().is_equal(inputPtrN->GetSpacing().GetVnlVector(), coordinateTol)) + { + spacingString.setf(std::ios::scientific); + spacingString.precision(7); + spacingString << "InputImage Spacing: " << inputPtr1->GetSpacing() << ", InputImage" << it.GetName() + << " Spacing: " << inputPtrN->GetSpacing() << std::endl; + spacingString << "\tTolerance: " << coordinateTol << std::endl; + } + if (!inputPtr1->GetDirection().GetVnlMatrix().as_ref().is_equal( + inputPtrN->GetDirection().GetVnlMatrix().as_ref(), Self::GetGlobalDefaultDirectionTolerance())) + { + directionString.setf(std::ios::scientific); + directionString.precision(7); + directionString << "InputImage Direction: " << inputPtr1->GetDirection() << ", InputImage" << it.GetName() + << " Direction: " << inputPtrN->GetDirection() << std::endl; + directionString << "\tTolerance: " << Self::GetGlobalDefaultDirectionTolerance() << std::endl; + } + itkExceptionMacro(<< "Inputs do not occupy the same physical space! " << std::endl + << originString.str() << spacingString.str() << directionString.str()); + } + } + } + } +} + +template +void +AttenuatedToExponentialCorrectionImageFilter::BeforeThreadedGenerateData() +{ + this->GetInterpolationWeightMultiplication().SetKRegionMinusAttenuationMapPtrDiff( + this->GetInput(2)->GetBufferPointer() - this->GetInput(1)->GetBufferPointer()); + this->GetSumAlongRay().SetBeforeKRegion(this->GetInterpolationWeightMultiplication().GetBeforeKRegion()); + this->GetSumAlongRay().SetLatestStepLength(this->GetInterpolationWeightMultiplication().GetLatestStepLength()); + this->GetProjectedValueAccumulation().SetBeforeKRegion( + this->GetInterpolationWeightMultiplication().GetBeforeKRegion()); + this->GetProjectedValueAccumulation().SetTraveledBeforeKRegion(this->GetSumAlongRay().GetTraveledBeforeKRegion()); + this->GetProjectedValueAccumulation().SetSpacing(this->GetInput(1)->GetSpacing()); + + // Get origin, i.e. point (0.,0.,0.), in voxel coordinates + typename Superclass::Superclass::GeometryType::ThreeDHomogeneousMatrixType volPPToIndex; + volPPToIndex = GetPhysicalPointToIndexMatrix(this->GetInput(1)).GetVnlMatrix(); + VectorType origin = itk::MakeVector(volPPToIndex[0][3], volPPToIndex[1][3], volPPToIndex[2][3]); + this->GetProjectedValueAccumulation().SetOrigin(origin); + + // Calculate mu0 as average in K region + itk::ImageRegionConstIterator itAtt(this->GetInput(1), this->GetInput(1)->GetLargestPossibleRegion()); + itk::ImageRegionConstIterator itKRegion(this->GetInput(2), + this->GetInput(2)->GetLargestPossibleRegion()); + unsigned long nVal = 0; + typename TOutputImage::PixelType val = 0.; + while (!itAtt.IsAtEnd()) + { + if (itKRegion.Get() != 0.) + { + val += itAtt.Get(); + nVal++; + } + ++itAtt; + ++itKRegion; + } + this->GetProjectedValueAccumulation().SetMu0(val / nVal); +} +} // end namespace rtk + +#endif diff --git a/include/rtkJosephForwardProjectionImageFilter.hxx b/include/rtkJosephForwardProjectionImageFilter.hxx index 9d0f7485b..06ffb15ae 100644 --- a/include/rtkJosephForwardProjectionImageFilter.hxx +++ b/include/rtkJosephForwardProjectionImageFilter.hxx @@ -240,7 +240,6 @@ JosephForwardProjectionImageFilterNext(), ++itOut) { typename InputRegionIterator::PointType pixelPosition = itIn->GetPixelPosition(); @@ -281,8 +280,8 @@ JosephForwardProjectionImageFilter nearDist) { // Compute and sort intersections: (n)earest and (f)arthest (p)points - np = pixelPosition + nearDist * dirVox; - fp = pixelPosition + farDist * dirVox; + typename BoxShape::VectorType np = pixelPosition + nearDist * dirVox; + typename BoxShape::VectorType fp = pixelPosition + farDist * dirVox; // Compute main nearest and farthest slice indices const int ns = itk::Math::rnd(np[mainDir]); @@ -315,6 +314,10 @@ JosephForwardProjectionImageFilter::One / dirVox[mainDir]; CoordRepType stepx = dirVox[notMainDirInf] * norm; CoordRepType stepy = dirVox[notMainDirSup] * norm; + + // Voxel to millimeters conversion vector + typename BoxShape::VectorType stepMM = this->GetInput(1)->GetSpacing(); + if (np[mainDir] > fp[mainDir]) { residualB *= -1; @@ -322,14 +325,13 @@ JosephForwardProjectionImageFilterGetInput(1)->GetSpacing()[notMainDirInf] * stepx; - stepMM[notMainDirSup] = this->GetInput(1)->GetSpacing()[notMainDirSup] * stepy; - stepMM[mainDir] = this->GetInput(1)->GetSpacing()[mainDir]; + stepMM[notMainDirInf] *= stepx; + stepMM[notMainDirSup] *= stepy; // Initialize the accumulation typename TOutputImage::PixelType sum = itk::NumericTraits::ZeroValue(); diff --git a/wrapping/rtkAttenuatedToExponentialCorrectionImageFilter.wrap b/wrapping/rtkAttenuatedToExponentialCorrectionImageFilter.wrap new file mode 100644 index 000000000..34a82d822 --- /dev/null +++ b/wrapping/rtkAttenuatedToExponentialCorrectionImageFilter.wrap @@ -0,0 +1,33 @@ +set(WRAPPER_AUTO_INCLUDE_HEADERS OFF) + +itk_wrap_named_class("rtk::Functor::InterpolationWeightMultiplicationAttToExp" "rtkFunctorInterpolationWeightMultiplicationAttToExp") + 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::ProjectedValueAccumulationAttToExp" "rtkProjectedValueAccumulationAttToExp") + 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::SumAttenuationForCorrection" "rtkSumAttenuationForCorrection") + 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}}3IWMATE${ITKM_${t}}${ITKM_${t}}PVAATE${ITKM_${t}}${ITKM_${t}}SAFC${ITKM_${t}}${ITKM_${t}}" + "itk::Image<${ITKT_${t}}, 3>, itk::Image< ${ITKT_${t}}, 3>, rtk::Functor::InterpolationWeightMultiplicationAttToExp<${ITKT_${t}}, ${ITKT_${t}}, ${ITKT_${t}}>, rtk::Functor::ProjectedValueAccumulationAttToExp<${ITKT_${t}}, ${ITKT_${t}}>, rtk::Functor::SumAttenuationForCorrection<${ITKT_${t}}, ${ITKT_${t}}>") + endforeach() +itk_end_wrap_class() + +itk_wrap_class("rtk::AttenuatedToExponentialCorrectionImageFilter" 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() +itk_end_wrap_class()