From fa52d03736909642bac254341bc730856f5b2f06 Mon Sep 17 00:00:00 2001 From: Simon Rit Date: Tue, 5 Nov 2024 22:21:05 +0100 Subject: [PATCH 1/5] BUG: Step vector in mm along ray was wrongly calculated Fixes a problem introduced in 581f3b3fd14230b50bbff697f0ee9b108f43bf7c. --- .../rtkJosephForwardProjectionImageFilter.hxx | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) 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(); From 3a272446df419cf285aacfe3f612a64e039578d5 Mon Sep 17 00:00:00 2001 From: Simon Rit Date: Wed, 6 Nov 2024 07:03:42 +0100 Subject: [PATCH 2/5] ENH: Add conversion from the attenuated 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 first input. $\mu_0$ is calculated as the average in the K region. --- .../rtkforwardprojections.cxx | 15 + .../rtkforwardprojections.ggo | 9 +- ...nuatedToExponentialCorrectionImageFilter.h | 377 ++++++++++++++++++ ...atedToExponentialCorrectionImageFilter.hxx | 210 ++++++++++ 4 files changed, 607 insertions(+), 4 deletions(-) create mode 100644 include/rtkAttenuatedToExponentialCorrectionImageFilter.h create mode 100644 include/rtkAttenuatedToExponentialCorrectionImageFilter.hxx 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..a0d51dee3 --- /dev/null +++ b/include/rtkAttenuatedToExponentialCorrectionImageFilter.h @@ -0,0 +1,377 @@ +/*========================================================================= + * + * 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 bool bInKRegion = ((p + m_KRegionMinusAttenuationMapPtrDiff)[i] != 0.); + m_PixelInKRegion[threadId] = m_PixelInKRegion[threadId] || bInKRegion; + m_LatestStepLength[threadId] = stepLengthInVoxel; + return weight * p[i]; + } + + void + SetKRegionMinusAttenuationMapPtrDiff(std::ptrdiff_t pd) + { + m_KRegionMinusAttenuationMapPtrDiff = pd; + } + + bool * + GetPixelInKRegion() + { + return m_PixelInKRegion; + } + + double * + GetLatestStepLength() + { + return m_LatestStepLength; + } + +private: + std::ptrdiff_t m_KRegionMinusAttenuationMapPtrDiff{}; + bool m_PixelInKRegion[itk::ITK_MAX_THREADS]{ false }; + 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 & stepInMM) + { + // Do nothing if we have passed the K region border + if (!m_BeforeKRegion[threadId]) + return; + + // Calculate the correction factor if we have reached the K region + if (m_PixelInKRegion[threadId]) + { + m_BeforeKRegion[threadId] = false; + } + else + { + sumValue += static_cast(volumeValue); + m_TraveledBeforeKRegion[threadId] += m_LatestStepLength[threadId]; + } + } + + double * + GetTraveledBeforeKRegion() + { + return m_TraveledBeforeKRegion; + } + + bool * + GetBeforeKRegion() + { + return m_BeforeKRegion; + } + + void + SetPixelInKRegion(bool * pixelInKRegion) + { + m_PixelInKRegion = pixelInKRegion; + } + + void + SetLatestStepLength(double * latestStepLength) + { + m_LatestStepLength = latestStepLength; + } + +private: + bool * m_PixelInKRegion{}; + bool m_BeforeKRegion[itk::ITK_MAX_THREADS]{ true }; + 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 & pixel, + const VectorType & pixelToSource, + const VectorType & nearestPoint, + const VectorType & 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 + tau * m_Mu0); + m_PixelInKRegion[threadId] = false; + m_BeforeKRegion[threadId] = true; + m_TraveledBeforeKRegion[threadId] = 0.; + } + + void + SetPixelInKRegion(bool * pixelInKRegion) + { + m_PixelInKRegion = pixelInKRegion; + } + + void + SetBeforeKRegion(bool * beforeKRegion) + { + m_BeforeKRegion = beforeKRegion; + } + + 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: + bool * m_PixelInKRegion{}; + bool * 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 first 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..44b13ff9a --- /dev/null +++ b/include/rtkAttenuatedToExponentialCorrectionImageFilter.hxx @@ -0,0 +1,210 @@ +/*========================================================================= + * + * 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().SetPixelInKRegion(this->GetInterpolationWeightMultiplication().GetPixelInKRegion()); + this->GetSumAlongRay().SetLatestStepLength(this->GetInterpolationWeightMultiplication().GetLatestStepLength()); + this->GetProjectedValueAccumulation().SetPixelInKRegion( + this->GetInterpolationWeightMultiplication().GetPixelInKRegion()); + this->GetProjectedValueAccumulation().SetBeforeKRegion(this->GetSumAlongRay().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 From ab7f8d344e13b119bd777a56698ed4d8b89193ed Mon Sep 17 00:00:00 2001 From: Simon Rit Date: Wed, 6 Nov 2024 23:27:25 +0100 Subject: [PATCH 3/5] ENH: Optimize the conversion from attenuated to exponential --- ...nuatedToExponentialCorrectionImageFilter.h | 71 +++++++------------ ...atedToExponentialCorrectionImageFilter.hxx | 7 +- 2 files changed, 28 insertions(+), 50 deletions(-) diff --git a/include/rtkAttenuatedToExponentialCorrectionImageFilter.h b/include/rtkAttenuatedToExponentialCorrectionImageFilter.h index a0d51dee3..beb490413 100644 --- a/include/rtkAttenuatedToExponentialCorrectionImageFilter.h +++ b/include/rtkAttenuatedToExponentialCorrectionImageFilter.h @@ -65,8 +65,9 @@ class ITK_TEMPLATE_EXPORT InterpolationWeightMultiplicationAttToExp const TInput * p, const int i) { - const bool bInKRegion = ((p + m_KRegionMinusAttenuationMapPtrDiff)[i] != 0.); - m_PixelInKRegion[threadId] = m_PixelInKRegion[threadId] || bInKRegion; + 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]; } @@ -77,10 +78,10 @@ class ITK_TEMPLATE_EXPORT InterpolationWeightMultiplicationAttToExp m_KRegionMinusAttenuationMapPtrDiff = pd; } - bool * - GetPixelInKRegion() + int * + GetBeforeKRegion() { - return m_PixelInKRegion; + return m_BeforeKRegion; } double * @@ -91,7 +92,7 @@ class ITK_TEMPLATE_EXPORT InterpolationWeightMultiplicationAttToExp private: std::ptrdiff_t m_KRegionMinusAttenuationMapPtrDiff{}; - bool m_PixelInKRegion[itk::ITK_MAX_THREADS]{ false }; + int m_BeforeKRegion[itk::ITK_MAX_THREADS]{ 1 }; double m_LatestStepLength[itk::ITK_MAX_THREADS]{}; }; @@ -125,22 +126,13 @@ class ITK_TEMPLATE_EXPORT SumAttenuationForCorrection } inline void - operator()(const ThreadIdType threadId, TOutput & sumValue, const TInput volumeValue, const VectorType & stepInMM) + operator()(const ThreadIdType threadId, + TOutput & sumValue, + const TInput volumeValue, + const VectorType & itkNotUsed(stepInMM)) { - // Do nothing if we have passed the K region border - if (!m_BeforeKRegion[threadId]) - return; - - // Calculate the correction factor if we have reached the K region - if (m_PixelInKRegion[threadId]) - { - m_BeforeKRegion[threadId] = false; - } - else - { - sumValue += static_cast(volumeValue); - m_TraveledBeforeKRegion[threadId] += m_LatestStepLength[threadId]; - } + sumValue += static_cast(volumeValue) * m_BeforeKRegion[threadId]; + m_TraveledBeforeKRegion[threadId] += m_LatestStepLength[threadId] * m_BeforeKRegion[threadId]; } double * @@ -149,16 +141,10 @@ class ITK_TEMPLATE_EXPORT SumAttenuationForCorrection return m_TraveledBeforeKRegion; } - bool * - GetBeforeKRegion() - { - return m_BeforeKRegion; - } - void - SetPixelInKRegion(bool * pixelInKRegion) + SetBeforeKRegion(int * pixelInKRegion) { - m_PixelInKRegion = pixelInKRegion; + m_BeforeKRegion = pixelInKRegion; } void @@ -168,8 +154,7 @@ class ITK_TEMPLATE_EXPORT SumAttenuationForCorrection } private: - bool * m_PixelInKRegion{}; - bool m_BeforeKRegion[itk::ITK_MAX_THREADS]{ true }; + int * m_BeforeKRegion{}; double m_TraveledBeforeKRegion[itk::ITK_MAX_THREADS]{ 0. }; double * m_LatestStepLength{}; }; @@ -210,10 +195,10 @@ class ITK_TEMPLATE_EXPORT ProjectedValueAccumulationAttToExp TOutput & output, const TOutput & rayCastValue, const VectorType & stepInMM, - const VectorType & pixel, + const VectorType & itkNotUsed(pixel), const VectorType & pixelToSource, const VectorType & nearestPoint, - const VectorType & farthestPoint) + 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. @@ -236,21 +221,16 @@ class ITK_TEMPLATE_EXPORT ProjectedValueAccumulationAttToExp 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 + tau * m_Mu0); - m_PixelInKRegion[threadId] = false; - m_BeforeKRegion[threadId] = true; - m_TraveledBeforeKRegion[threadId] = 0.; - } - void - SetPixelInKRegion(bool * pixelInKRegion) - { - m_PixelInKRegion = pixelInKRegion; + // Reinitialize for next ray + m_BeforeKRegion[threadId] = 1; + m_TraveledBeforeKRegion[threadId] = 0.; } void - SetBeforeKRegion(bool * beforeKRegion) + SetBeforeKRegion(int * pixelInKRegion) { - m_BeforeKRegion = beforeKRegion; + m_BeforeKRegion = pixelInKRegion; } void @@ -278,8 +258,7 @@ class ITK_TEMPLATE_EXPORT ProjectedValueAccumulationAttToExp } private: - bool * m_PixelInKRegion{}; - bool * m_BeforeKRegion{}; + int * m_BeforeKRegion{}; double * m_TraveledBeforeKRegion{}; TOutput m_Mu0{}; VectorType m_Origin{}; @@ -294,7 +273,7 @@ class ITK_TEMPLATE_EXPORT ProjectedValueAccumulationAttToExp * 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 first input. $\mu_0$ is + * 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 diff --git a/include/rtkAttenuatedToExponentialCorrectionImageFilter.hxx b/include/rtkAttenuatedToExponentialCorrectionImageFilter.hxx index 44b13ff9a..d1732d375 100644 --- a/include/rtkAttenuatedToExponentialCorrectionImageFilter.hxx +++ b/include/rtkAttenuatedToExponentialCorrectionImageFilter.hxx @@ -173,11 +173,10 @@ AttenuatedToExponentialCorrectionImageFilterGetInterpolationWeightMultiplication().SetKRegionMinusAttenuationMapPtrDiff( this->GetInput(2)->GetBufferPointer() - this->GetInput(1)->GetBufferPointer()); - this->GetSumAlongRay().SetPixelInKRegion(this->GetInterpolationWeightMultiplication().GetPixelInKRegion()); + this->GetSumAlongRay().SetBeforeKRegion(this->GetInterpolationWeightMultiplication().GetBeforeKRegion()); this->GetSumAlongRay().SetLatestStepLength(this->GetInterpolationWeightMultiplication().GetLatestStepLength()); - this->GetProjectedValueAccumulation().SetPixelInKRegion( - this->GetInterpolationWeightMultiplication().GetPixelInKRegion()); - this->GetProjectedValueAccumulation().SetBeforeKRegion(this->GetSumAlongRay().GetBeforeKRegion()); + this->GetProjectedValueAccumulation().SetBeforeKRegion( + this->GetInterpolationWeightMultiplication().GetBeforeKRegion()); this->GetProjectedValueAccumulation().SetTraveledBeforeKRegion(this->GetSumAlongRay().GetTraveledBeforeKRegion()); this->GetProjectedValueAccumulation().SetSpacing(this->GetInput(1)->GetSpacing()); From 3f14eb95176c3d6b0b10c030c4d81af2d5a7824e Mon Sep 17 00:00:00 2001 From: Simon Rit Date: Thu, 5 Dec 2024 23:04:18 +0100 Subject: [PATCH 4/5] BUG: Add missing step length conversion to mm --- include/rtkAttenuatedToExponentialCorrectionImageFilter.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/rtkAttenuatedToExponentialCorrectionImageFilter.h b/include/rtkAttenuatedToExponentialCorrectionImageFilter.h index beb490413..f671378fc 100644 --- a/include/rtkAttenuatedToExponentialCorrectionImageFilter.h +++ b/include/rtkAttenuatedToExponentialCorrectionImageFilter.h @@ -220,7 +220,7 @@ class ITK_TEMPLATE_EXPORT ProjectedValueAccumulationAttToExp 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 + tau * m_Mu0); + output = input * std::exp(rayCastValue * stepInMM.GetNorm() + tau * m_Mu0); // Reinitialize for next ray m_BeforeKRegion[threadId] = 1; From c454a1e7809e5e90563c9a01abe03f8c127b3ffe Mon Sep 17 00:00:00 2001 From: Simon Rit Date: Mon, 9 Dec 2024 22:37:01 +0100 Subject: [PATCH 5/5] ENH: Python wrapping of AttenuatedToExponentialCorrectionImageFilter --- ...tedToExponentialCorrectionImageFilter.wrap | 33 +++++++++++++++++++ 1 file changed, 33 insertions(+) create mode 100644 wrapping/rtkAttenuatedToExponentialCorrectionImageFilter.wrap 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()