diff --git a/applications/rtkforwardprojections/rtkforwardprojections.cxx b/applications/rtkforwardprojections/rtkforwardprojections.cxx index 4cebb80f4..d34405306 100644 --- a/applications/rtkforwardprojections/rtkforwardprojections.cxx +++ b/applications/rtkforwardprojections/rtkforwardprojections.cxx @@ -21,6 +21,7 @@ #include "rtkThreeDCircularProjectionGeometryXMLFile.h" #include "rtkJosephForwardProjectionImageFilter.h" +#include "rtkMaximumIntensityProjectionImageFilter.h" #include "rtkJosephForwardAttenuatedProjectionImageFilter.h" #include "rtkZengForwardProjectionImageFilter.h" #ifdef RTK_USE_CUDA @@ -114,6 +115,9 @@ main(int argc, char * argv[]) case (fp_arg_Zeng): forwardProjection = rtk::ZengForwardProjectionImageFilter::New(); break; + case (fp_arg_MIP): + forwardProjection = rtk::MaximumIntensityProjectionImageFilter::New(); + break; case (fp_arg_CudaRayCast): #ifdef RTK_USE_CUDA forwardProjection = rtk::CudaForwardProjectionImageFilter::New(); @@ -147,6 +151,12 @@ main(int argc, char * argv[]) forwardProjection.GetPointer()) ->SetInferiorClipImage(inferiorClipImage); } + else if (args_info.fp_arg == fp_arg_MIP) + { + dynamic_cast *>( + forwardProjection.GetPointer()) + ->SetInferiorClipImage(inferiorClipImage); + } } if (args_info.superiorclipimage_given) { @@ -162,6 +172,12 @@ main(int argc, char * argv[]) forwardProjection.GetPointer()) ->SetSuperiorClipImage(superiorClipImage); } + else if (args_info.fp_arg == fp_arg_MIP) + { + dynamic_cast *>( + forwardProjection.GetPointer()) + ->SetSuperiorClipImage(superiorClipImage); + } } if (args_info.sigmazero_given && args_info.fp_arg == fp_arg_Zeng) dynamic_cast *>( diff --git a/applications/rtkforwardprojections/rtkforwardprojections.ggo b/applications/rtkforwardprojections/rtkforwardprojections.ggo index 49085c552..d867b7690 100644 --- a/applications/rtkforwardprojections/rtkforwardprojections.ggo +++ b/applications/rtkforwardprojections/rtkforwardprojections.ggo @@ -10,9 +10,9 @@ 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" enum no default="Joseph" +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 "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 \ No newline at end of file +option "superiorclipimage" - "Value of the superior clip of the ray for each pixel of the projections (only with Joseph-based projector)" string no diff --git a/include/rtkJosephForwardAttenuatedProjectionImageFilter.h b/include/rtkJosephForwardAttenuatedProjectionImageFilter.h index 56729146a..f1e17283d 100644 --- a/include/rtkJosephForwardAttenuatedProjectionImageFilter.h +++ b/include/rtkJosephForwardAttenuatedProjectionImageFilter.h @@ -135,8 +135,8 @@ class ITK_TEMPLATE_EXPORT ComputeAttenuationCorrection return !(*this != other); } - inline TOutput - operator()(const ThreadIdType threadId, const TInput volumeValue, const VectorType & stepInMM) + inline void + operator()(const ThreadIdType threadId, TOutput & sumValue, const TInput volumeValue, const VectorType & stepInMM) { TInput ex2 = exp(-m_AttenuationRay[threadId] * stepInMM.GetNorm()); TInput wf; @@ -152,7 +152,7 @@ class ITK_TEMPLATE_EXPORT ComputeAttenuationCorrection m_Ex1[threadId] = ex2; m_AttenuationPixel[threadId] = 0; - return wf * volumeValue; + sumValue += wf * volumeValue; } void diff --git a/include/rtkJosephForwardProjectionImageFilter.h b/include/rtkJosephForwardProjectionImageFilter.h index 2b5115ed8..87a29dd55 100644 --- a/include/rtkJosephForwardProjectionImageFilter.h +++ b/include/rtkJosephForwardProjectionImageFilter.h @@ -93,10 +93,13 @@ class ITK_TEMPLATE_EXPORT SumAlongRay return !(*this != other); } - inline TOutput - operator()(const ThreadIdType itkNotUsed(threadId), const TInput volumeValue, const VectorType & itkNotUsed(stepInMM)) + inline void + operator()(const ThreadIdType itkNotUsed(threadId), + TOutput & sumValue, + const TInput volumeValue, + const VectorType & itkNotUsed(stepInMM)) { - return volumeValue; + sumValue += static_cast(volumeValue); } }; diff --git a/include/rtkJosephForwardProjectionImageFilter.hxx b/include/rtkJosephForwardProjectionImageFilter.hxx index 3f706542f..9d0f7485b 100644 --- a/include/rtkJosephForwardProjectionImageFilter.hxx +++ b/include/rtkJosephForwardProjectionImageFilter.hxx @@ -351,7 +351,7 @@ JosephForwardProjectionImageFilterGet(), itOut.Value(), sum, stepMM, pixelPosition, dirVox, np, fp); @@ -538,7 +538,6 @@ JosephForwardProjectionImageFilter +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); + } + + inline void + operator()(const ThreadIdType itkNotUsed(threadId), + TOutput & mipValue, + const TInput volumeValue, + const VectorType & itkNotUsed(stepInMM)) + { + TOutput 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; + + bool + operator!=(const MaximumIntensityProjectedValueAccumulation &) const + { + return false; + } + bool + operator==(const MaximumIntensityProjectedValueAccumulation & other) const + { + return !(*this != other); + } + + inline void + operator()(const ThreadIdType itkNotUsed(threadId), + const TInput & input, + TOutput & output, + const TOutput & rayCastValue, + const VectorType & stepInMM, + const VectorType & itkNotUsed(source), + const VectorType & itkNotUsed(sourceToPixel), + const VectorType & itkNotUsed(nearestPoint), + const VectorType & itkNotUsed(farthestPoint)) const + { + TOutput tmp = static_cast(input); + if (tmp < rayCastValue) + { + tmp = rayCastValue; + } + output = tmp * stepInMM.GetNorm(); + } +}; + +} // end namespace Functor + + +/** \class MaximumIntensityProjectionImageFilter + * \brief MIP filter. + * + * Performs a MIP forward projection, i.e. calculation of a maximum intensity + * step along the x-ray line. + * + * \author Mikhail Polkovnikov + * + * \ingroup RTK Projector + */ + +template ::ValueType>, + class TProjectedValueAccumulation = + Functor::MaximumIntensityProjectedValueAccumulation, + class TSumAlongRay = + Functor::MaximumIntensityAlongRay> +class ITK_TEMPLATE_EXPORT MaximumIntensityProjectionImageFilter + : public JosephForwardProjectionImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(MaximumIntensityProjectionImageFilter); + + /** Standard class type alias. */ + using Self = MaximumIntensityProjectionImageFilter; + using Pointer = itk::SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ +#ifdef itkOverrideGetNameOfClassMacro + itkOverrideGetNameOfClassMacro(MaximumIntensityProjectionImageFilter); +#else + itkTypeMacro(MaximumIntensityProjectionImageFilter, JosephForwardProjectionImageFilter); +#endif + +protected: + MaximumIntensityProjectionImageFilter() = default; + ~MaximumIntensityProjectionImageFilter() override = default; +}; +} // end namespace rtk + +#endif diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index d0a1a8b1a..197562c55 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -89,6 +89,7 @@ rtk_add_test(rtkForwardProjectionTest rtkforwardprojectiontest.cxx) rtk_add_cuda_test(rtkForwardProjectionCudaTest rtkforwardprojectiontest.cxx) rtk_add_test(rtkForwardAttenuatedProjectionTest rtkforwardattenuatedprojectiontest.cxx) rtk_add_test(rtkZengProjectionTest rtkzengforwardprojectiontest.cxx) +rtk_add_test(rtkMaximumIntensityProjectionTest rtkmaximumintensityprojectiontest.cxx) rtk_add_test(rtkGeometryFileTest rtkgeometryfiletest.cxx) @@ -301,4 +302,5 @@ rtk_add_test(rtkBioscanTest rtkbioscantest.cxx if(ITK_WRAP_PYTHON) itk_python_add_test(NAME rtkFirstReconstructionPythonTest COMMAND rtkFirstReconstruction.py ${CMAKE_CURRENT_BINARY_DIR}/rtkFirstReconstruction.mha) + itk_python_add_test(NAME rtkMaximumIntensityPythonTest COMMAND rtkMaximumIntensity.py ${CMAKE_CURRENT_BINARY_DIR}/rtkFirstReconstruction.mha) endif() diff --git a/test/rtkMaximumIntensity.py b/test/rtkMaximumIntensity.py new file mode 100644 index 000000000..58492e169 --- /dev/null +++ b/test/rtkMaximumIntensity.py @@ -0,0 +1,40 @@ +#!/usr/bin/env python +from __future__ import print_function +import itk +from itk import RTK as rtk +import sys +import os + +TImageType = itk.Image[itk.F,3] +# Defines the RTK geometry object +geometry = rtk.ThreeDCircularProjectionGeometry.New() +numberOfProjections = 1 +geometry.AddProjection(700, 800, 0) +# Constant image sources +# Create MIP Forward Projector volume +volInput = rtk.ConstantImageSource[TImageType].New() +origin = [ 0., 0., 0. ] +size = [ 64, 64, 64 ] +sizeOutput = [ 200, 200, numberOfProjections ] +spacing = [ 4.0, 4.0, 4.0 ] +spacingOutput = [ 1.0, 1.0, 1.0 ] +volInput.SetOrigin( origin ) +volInput.SetSpacing( spacing ) +volInput.SetSize( size ) +volInput.SetConstant(1.0) +volInput.UpdateOutputInformation() +volInputSource = volInput.GetOutput() +# Initialization Imager volume +projInput = rtk.ConstantImageSource[TImageType].New() +projInput.SetOrigin( origin ) +projInput.SetSpacing( spacingOutput ) +projInput.SetSize( sizeOutput ) +projInput.SetConstant(0.) +projInput.Update() +projInputSource = projInput.GetOutput() +# MIP Forward Projection filter +mip = rtk.MaximumIntensityProjectionImageFilter[TImageType, TImageType].New() +mip.SetGeometry( geometry ) +mip.SetInput(volInputSource) +mip.SetInput(1, projInputSource) +mipImage = mip.GetOutput() diff --git a/test/rtkmaximumintensityprojectiontest.cxx b/test/rtkmaximumintensityprojectiontest.cxx new file mode 100644 index 000000000..1bdaae9ad --- /dev/null +++ b/test/rtkmaximumintensityprojectiontest.cxx @@ -0,0 +1,109 @@ +#include "rtkTest.h" +#include "rtkThreeDCircularProjectionGeometryXMLFile.h" +#include "rtkRayBoxIntersectionImageFilter.h" +#include "rtkConstantImageSource.h" +#include "rtkMaximumIntensityProjectionImageFilter.h" +#include +#include +#include +#include +#include + + +/** + * \file rtkmaximumintensityprojectiontest.cxx + * + * \brief Functional test for MIP forward filter + * + * The test projects a volume filled with ones. + * + * \author Mikhail Polkovnikov + */ + +int +main(int, char **) +{ + constexpr unsigned int Dimension = 3; + using OutputPixelType = float; + using OutputImageType = itk::Image; + + using VectorType = itk::Vector; + constexpr unsigned int NumberOfProjectionImages = 1; + + // Constant image sources + using ConstantImageSourceType = rtk::ConstantImageSource; + ConstantImageSourceType::PointType origin; + ConstantImageSourceType::SizeType size; + ConstantImageSourceType::SizeType sizeOutput; + ConstantImageSourceType::SpacingType spacing; + ConstantImageSourceType::SpacingType spacingOutput; + + // Create MIP Forward Projector volume input. + const ConstantImageSourceType::Pointer volInput = ConstantImageSourceType::New(); + origin[0] = 0; + origin[1] = 0; + origin[2] = 0; + size[0] = 64; + size[1] = 64; + size[2] = 64; + spacing[0] = 4.; + spacing[1] = 4.; + spacing[2] = 4.; + + sizeOutput[0] = 200; + sizeOutput[1] = 200; + sizeOutput[2] = 200; + spacingOutput[0] = 1.; + spacingOutput[1] = 1.; + spacingOutput[2] = 1.; + + volInput->SetOrigin(origin); + volInput->SetSpacing(spacing); + volInput->SetSize(size); + volInput->SetConstant(1.); + volInput->UpdateOutputInformation(); + + // Initialization Imager Volume + const ConstantImageSourceType::Pointer projInput = ConstantImageSourceType::New(); + sizeOutput[2] = NumberOfProjectionImages; + projInput->SetOrigin(origin); + projInput->SetSpacing(spacingOutput); + projInput->SetSize(sizeOutput); + projInput->SetConstant(0.); + projInput->Update(); + + // MIP Forward Projection filter + using MIPType = rtk::MaximumIntensityProjectionImageFilter; + MIPType::Pointer mipfp = MIPType::New(); + mipfp->SetInput(projInput->GetOutput()); + mipfp->SetInput(1, volInput->GetOutput()); + + using GeometryType = rtk::ThreeDCircularProjectionGeometry; + GeometryType::Pointer geometry = GeometryType::New(); + geometry->AddProjection(700, 800, 0); + + mipfp->SetGeometry(geometry); + mipfp->Update(); + + using ConstIteratorType = itk::ImageRegionConstIterator; + ConstIteratorType inputIt(mipfp->GetOutput(), mipfp->GetOutput()->GetRequestedRegion()); + + inputIt.GoToBegin(); + + bool res = false; + while (!inputIt.IsAtEnd()) + { + OutputPixelType pixel = inputIt.Get(); + if (pixel < 4. || pixel > 4.25) + { + res = true; + } + ++inputIt; + } + if (res) + { + return EXIT_FAILURE; + } + std::cout << "\n\nTest PASSED! " << std::endl; + return EXIT_SUCCESS; +} diff --git a/wrapping/rtkMaximumIntensityProjectionImageFilter.wrap b/wrapping/rtkMaximumIntensityProjectionImageFilter.wrap new file mode 100644 index 000000000..b4b1acce6 --- /dev/null +++ b/wrapping/rtkMaximumIntensityProjectionImageFilter.wrap @@ -0,0 +1,30 @@ +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) + 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()