Skip to content

Commit

Permalink
ENH: MIP image filter for DRR image calculation
Browse files Browse the repository at this point in the history
The MaximumIntensityForwardProjectionImageFilter is derived from
JosephForwardProjectionImageFilter and computes maximum
intensity step along the ray casting on the projection.
  • Loading branch information
Mikhail Polkovnikov authored and SimonRit committed Mar 13, 2024
1 parent 52c5d9f commit 081e017
Show file tree
Hide file tree
Showing 10 changed files with 379 additions and 13 deletions.
16 changes: 16 additions & 0 deletions applications/rtkforwardprojections/rtkforwardprojections.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@

#include "rtkThreeDCircularProjectionGeometryXMLFile.h"
#include "rtkJosephForwardProjectionImageFilter.h"
#include "rtkMaximumIntensityProjectionImageFilter.h"
#include "rtkJosephForwardAttenuatedProjectionImageFilter.h"
#include "rtkZengForwardProjectionImageFilter.h"
#ifdef RTK_USE_CUDA
Expand Down Expand Up @@ -114,6 +115,9 @@ main(int argc, char * argv[])
case (fp_arg_Zeng):
forwardProjection = rtk::ZengForwardProjectionImageFilter<OutputImageType, OutputImageType>::New();
break;
case (fp_arg_MIP):
forwardProjection = rtk::MaximumIntensityProjectionImageFilter<OutputImageType, OutputImageType>::New();
break;
case (fp_arg_CudaRayCast):
#ifdef RTK_USE_CUDA
forwardProjection = rtk::CudaForwardProjectionImageFilter<OutputImageType, OutputImageType>::New();
Expand Down Expand Up @@ -147,6 +151,12 @@ main(int argc, char * argv[])
forwardProjection.GetPointer())
->SetInferiorClipImage(inferiorClipImage);
}
else if (args_info.fp_arg == fp_arg_MIP)
{
dynamic_cast<rtk::MaximumIntensityProjectionImageFilter<OutputImageType, OutputImageType> *>(
forwardProjection.GetPointer())
->SetInferiorClipImage(inferiorClipImage);
}
}
if (args_info.superiorclipimage_given)
{
Expand All @@ -162,6 +172,12 @@ main(int argc, char * argv[])
forwardProjection.GetPointer())
->SetSuperiorClipImage(superiorClipImage);
}
else if (args_info.fp_arg == fp_arg_MIP)
{
dynamic_cast<rtk::MaximumIntensityProjectionImageFilter<OutputImageType, OutputImageType> *>(
forwardProjection.GetPointer())
->SetSuperiorClipImage(superiorClipImage);
}
}
if (args_info.sigmazero_given && args_info.fp_arg == fp_arg_Zeng)
dynamic_cast<rtk::ZengForwardProjectionImageFilter<OutputImageType, OutputImageType> *>(
Expand Down
4 changes: 2 additions & 2 deletions applications/rtkforwardprojections/rtkforwardprojections.ggo
Original file line number Diff line number Diff line change
Expand Up @@ -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
option "superiorclipimage" - "Value of the superior clip of the ray for each pixel of the projections (only with Joseph-based projector)" string no
6 changes: 3 additions & 3 deletions include/rtkJosephForwardAttenuatedProjectionImageFilter.h
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -152,7 +152,7 @@ class ITK_TEMPLATE_EXPORT ComputeAttenuationCorrection

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

void
Expand Down
9 changes: 6 additions & 3 deletions include/rtkJosephForwardProjectionImageFilter.h
Original file line number Diff line number Diff line change
Expand Up @@ -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<TOutput>(volumeValue);
}
};

Expand Down
9 changes: 4 additions & 5 deletions include/rtkJosephForwardProjectionImageFilter.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -351,7 +351,7 @@ JosephForwardProjectionImageFilter<TInputImage,
miny,
maxx,
maxy);
sum += m_SumAlongRay(threadId, volumeValue, stepMM);
m_SumAlongRay(threadId, sum, volumeValue, stepMM);
}
else
{
Expand All @@ -370,7 +370,7 @@ JosephForwardProjectionImageFilter<TInputImage,
miny,
maxx,
maxy);
sum += m_SumAlongRay(threadId, volumeValue, stepMM);
m_SumAlongRay(threadId, sum, volumeValue, stepMM);

// Move to next main direction slice
pxiyi += offsetz;
Expand All @@ -385,7 +385,7 @@ JosephForwardProjectionImageFilter<TInputImage,
{
volumeValue =
BilinearInterpolation(threadId, 1., pxiyi, pxsyi, pxiys, pxsys, currentx, currenty, offsetx, offsety);
sum += m_SumAlongRay(threadId, volumeValue, stepMM);
m_SumAlongRay(threadId, sum, volumeValue, stepMM);

// Move to next main direction slice
pxiyi += offsetz;
Expand All @@ -411,7 +411,7 @@ JosephForwardProjectionImageFilter<TInputImage,
miny,
maxx,
maxy);
sum += m_SumAlongRay(threadId, volumeValue, stepMM);
m_SumAlongRay(threadId, sum, volumeValue, stepMM);
}
// Accumulate
m_ProjectedValueAccumulation(threadId, itIn->Get(), itOut.Value(), sum, stepMM, pixelPosition, dirVox, np, fp);
Expand Down Expand Up @@ -538,7 +538,6 @@ JosephForwardProjectionImageFilter<TInputImage,

return (result);
}

} // end namespace rtk

#endif
167 changes: 167 additions & 0 deletions include/rtkMaximumIntensityProjectionImageFilter.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,167 @@
/*=========================================================================
*
* 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_h
#define rtkMaximumIntensityProjectionImageFilter_h

#include "rtkJosephForwardProjectionImageFilter.h"

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 TInput, class TOutput>
class ITK_TEMPLATE_EXPORT MaximumIntensityAlongRay
{
public:
using VectorType = itk::Vector<double, 3>;

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<TOutput>(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 TInput, class TOutput>
class ITK_TEMPLATE_EXPORT MaximumIntensityProjectedValueAccumulation
{
public:
using VectorType = itk::Vector<double, 3>;

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<TOutput>(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 <class TInputImage,
class TOutputImage,
class TInterpolationWeightMultiplication = Functor::InterpolationWeightMultiplication<
typename TInputImage::PixelType,
typename itk::PixelTraits<typename TInputImage::PixelType>::ValueType>,
class TProjectedValueAccumulation =
Functor::MaximumIntensityProjectedValueAccumulation<typename TInputImage::PixelType,
typename TOutputImage::PixelType>,
class TSumAlongRay =
Functor::MaximumIntensityAlongRay<typename TInputImage::PixelType, typename TOutputImage::PixelType>>
class ITK_TEMPLATE_EXPORT MaximumIntensityProjectionImageFilter
: public JosephForwardProjectionImageFilter<TInputImage,
TOutputImage,
TInterpolationWeightMultiplication,
TProjectedValueAccumulation,
TSumAlongRay>
{
public:
ITK_DISALLOW_COPY_AND_MOVE(MaximumIntensityProjectionImageFilter);

/** Standard class type alias. */
using Self = MaximumIntensityProjectionImageFilter;
using Pointer = itk::SmartPointer<Self>;

/** 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
2 changes: 2 additions & 0 deletions test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down Expand Up @@ -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()
40 changes: 40 additions & 0 deletions test/rtkMaximumIntensity.py
Original file line number Diff line number Diff line change
@@ -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()
Loading

0 comments on commit 081e017

Please sign in to comment.