ENH: MIP image filter for DRR image calculation

MIPForwardProjectionImageFilter is derived from
JosephForwardProjectionImageFilter and computes maximum
intensity step along the ray casting on the projection

include/rtkJosephForwardProjectionImageFilter.h

@@ -93,13 +93,56 @@ class ITK_TEMPLATE_EXPORT SumAlongRay
     return !(*this != other);
   }
 
-  inline TOutput
+  virtual inline TOutput
   operator()(const ThreadIdType itkNotUsed(threadId), const TInput volumeValue, const VectorType & itkNotUsed(stepInMM))
   {
     return volumeValue;
   }
 };
 
+/** \class MIPAlongRay
+ * \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 MIPAlongRay : public SumAlongRay< TInput, TOutput >
+{
+public:
+  using VectorType = itk::Vector<double, 3>;
+
+  MIPAlongRay() = default;
+  ~MIPAlongRay() = default;
+  bool
+  operator!=(const MIPAlongRay &) const
+  {
+    return false;
+  }
+  bool
+  operator==(const MIPAlongRay & other) const
+  {
+    return !(*this != other);
+  }
+
+  inline TOutput
+  operator()(const ThreadIdType itkNotUsed(threadId), const TInput volumeValue, const VectorType & itkNotUsed(stepInMM)) override
+  {
+    return volumeValue;
+  }
+
+  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 ProjectedValueAccumulation
  * \brief Function to accumulate the ray casting on the projection.
  *
@@ -126,7 +169,7 @@ class ITK_TEMPLATE_EXPORT ProjectedValueAccumulation
     return !(*this != other);
   }
 
-  inline void
+  virtual inline void
   operator()(const ThreadIdType itkNotUsed(threadId),
              const TInput &     input,
              TOutput &          output,
@@ -141,6 +184,46 @@ class ITK_TEMPLATE_EXPORT ProjectedValueAccumulation
   }
 };
 
+/** \class MIPProjectedValueAccumulation
+ * \brief Function to calculate maximum intensity projection (MIP) or
+ * minimum intensity projection (MinIP) the ray casting on the projection.
+ *
+ * \author Simon Rit
+ *
+ * \ingroup RTK Functions
+ */
+template <class TInput, class TOutput>
+class ITK_TEMPLATE_EXPORT MIPProjectedValueAccumulation : public ProjectedValueAccumulation< TInput, TOutput >
+{
+public:
+  using VectorType = itk::Vector<double, 3>;
+
+  bool
+  operator!=(const MIPProjectedValueAccumulation &) const
+  {
+    return false;
+  }
+  bool
+  operator==(const MIPProjectedValueAccumulation & other) const
+  {
+    return !(*this != other);
+  }
+
+  inline void
+  operator()(const ThreadIdType itkNotUsed(threadId),
+             const TInput &     itkNotUsed(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 override
+  {
+    output = rayCastValue * stepInMM.GetNorm();
+  }
+};
+
 } // end namespace Functor
 
 
@@ -342,7 +425,7 @@ class ITK_TEMPLATE_EXPORT JosephForwardProjectionImageFilter
                                  const double           maxx,
                                  const double           maxy);
 
-private:
+protected:
   // Functors
   TInterpolationWeightMultiplication m_InterpolationWeightMultiplication;
   TProjectedValueAccumulation        m_ProjectedValueAccumulation;
@@ -351,6 +434,52 @@ class ITK_TEMPLATE_EXPORT JosephForwardProjectionImageFilter
   double                             m_SuperiorClip{ 1. };
 };
 
+template <class TInputImage,
+          class TOutputImage,
+          class TInterpolationWeightMultiplication = Functor::InterpolationWeightMultiplication<
+            typename TInputImage::PixelType,
+            typename itk::PixelTraits<typename TInputImage::PixelType>::ValueType>,
+          class TProjectedValueAccumulation =
+            Functor::MIPProjectedValueAccumulation<typename TInputImage::PixelType, typename TOutputImage::PixelType>,
+          class TSumAlongRay = Functor::MIPAlongRay<typename TInputImage::PixelType, typename TOutputImage::PixelType>>
+class ITK_TEMPLATE_EXPORT MIPForwardProjectionImageFilter
+  : public JosephForwardProjectionImageFilter<TInputImage, TOutputImage, TInterpolationWeightMultiplication, TProjectedValueAccumulation, TSumAlongRay>
+{
+public:
+  ITK_DISALLOW_COPY_AND_MOVE(MIPForwardProjectionImageFilter);
+
+  /** Standard class type alias. */
+  using Self = MIPForwardProjectionImageFilter;
+  using Superclass = JosephForwardProjectionImageFilter<TInputImage, TOutputImage, TInterpolationWeightMultiplication, TProjectedValueAccumulation, TSumAlongRay>;
+  using Pointer = itk::SmartPointer<Self>;
+  using ConstPointer = itk::SmartPointer<const Self>;
+  using InputPixelType = typename TInputImage::PixelType;
+  using OutputPixelType = typename TOutputImage::PixelType;
+  using OutputImageRegionType = typename TOutputImage::RegionType;
+  using CoordRepType = double;
+  using VectorType = itk::Vector<CoordRepType, TInputImage::ImageDimension>;
+  using TClipImageType = itk::Image<double, TOutputImage::ImageDimension>;
+  using TClipImagePointer = typename TClipImageType::Pointer;
+
+  /** Method for creation through the object factory. */
+  itkNewMacro(Self);
+
+  /** Run-time type information (and related methods). */
+#ifdef itkOverrideGetNameOfClassMacro
+  itkOverrideGetNameOfClassMacro(MIPForwardProjectionImageFilter);
+#else
+  itkTypeMacro(MIPForwardProjectionImageFilter, JosephForwardProjectionImageFilter);
+#endif
+
+protected:
+  MIPForwardProjectionImageFilter() = default;
+  ~MIPForwardProjectionImageFilter() override = default;
+
+  void
+  ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, ThreadIdType threadId) override;
+
+};
+
 } // end namespace rtk
 
 #ifndef ITK_MANUAL_INSTANTIATION