rtkJosephForwardProjectionImageFilter.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  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 rtkJosephForwardProjectionImageFilter_h
#define rtkJosephForwardProjectionImageFilter_h

#include "rtkConfiguration.h"
#include "rtkForwardProjectionImageFilter.h"
#include "rtkMacro.h"
#include <itkPixelTraits.h>

#include "rtkProjectionsRegionConstIteratorRayBased.h"

#include <itkVectorImage.h>
namespace rtk
{
namespace Functor
{
template <class TInput, class TCoordRepType, class TOutput = TInput>
class ITK_TEMPLATE_EXPORT InterpolationWeightMultiplication
{
public:
  InterpolationWeightMultiplication() = default;
  ~InterpolationWeightMultiplication() = default;
  bool
  operator!=(const InterpolationWeightMultiplication &) const
  {
    return false;
  }
  bool
  operator==(const InterpolationWeightMultiplication & other) const
  {
    return !(*this != other);
  }

  inline TOutput
  operator()(const ThreadIdType  itkNotUsed(threadId),
             const double        itkNotUsed(stepLengthInVoxel),
             const TCoordRepType weight,
             const TInput *      p,
             const int           i) const
  {
    return weight * p[i];
  }
};

template <class TInput, class TOutput>
class ITK_TEMPLATE_EXPORT SumAlongRay
{
public:
  using VectorType = itk::Vector<double, 3>;

  SumAlongRay() = default;
  ~SumAlongRay() = default;
  bool
  operator!=(const SumAlongRay &) const
  {
    return false;
  }
  bool
  operator==(const SumAlongRay & other) const
  {
    return !(*this != other);
  }

  inline void
  operator()(const ThreadIdType itkNotUsed(threadId),
             TOutput &          sumValue,
             const TInput       volumeValue,
             const VectorType & itkNotUsed(stepInMM))
  {
    sumValue += static_cast<TOutput>(volumeValue);
  }
};

template <class TInput, class TOutput>
class ITK_TEMPLATE_EXPORT ProjectedValueAccumulation
{
public:
  using VectorType = itk::Vector<double, 3>;

  ProjectedValueAccumulation() = default;
  ~ProjectedValueAccumulation() = default;
  bool
  operator!=(const ProjectedValueAccumulation &) const
  {
    return false;
  }
  bool
  operator==(const ProjectedValueAccumulation & 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
  {
    output = input + rayCastValue * stepInMM.GetNorm();
  }
};

} // end namespace Functor


template <class TInputImage,
          class TOutputImage,
          class TInterpolationWeightMultiplication = Functor::InterpolationWeightMultiplication<
            typename TInputImage::PixelType,
            typename itk::PixelTraits<typename TInputImage::PixelType>::ValueType>,
          class TProjectedValueAccumulation =
            Functor::ProjectedValueAccumulation<typename TInputImage::PixelType, typename TOutputImage::PixelType>,
          class TSumAlongRay = Functor::SumAlongRay<typename TInputImage::PixelType, typename TOutputImage::PixelType>>
class ITK_TEMPLATE_EXPORT JosephForwardProjectionImageFilter
  : public ForwardProjectionImageFilter<TInputImage, TOutputImage>
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(JosephForwardProjectionImageFilter);

  using Self = JosephForwardProjectionImageFilter;
  using Superclass = ForwardProjectionImageFilter<TInputImage, TOutputImage>;
  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;

  itkNewMacro(Self);

#ifdef itkOverrideGetNameOfClassMacro
  itkOverrideGetNameOfClassMacro(JosephForwardProjectionImageFilter);
#else
  itkTypeMacro(JosephForwardProjectionImageFilter, ForwardProjectionImageFilter);
#endif


  TInterpolationWeightMultiplication &
  GetInterpolationWeightMultiplication()
  {
    return m_InterpolationWeightMultiplication;
  }
  const TInterpolationWeightMultiplication &
  GetInterpolationWeightMultiplication() const
  {
    return m_InterpolationWeightMultiplication;
  }
  void
  SetInterpolationWeightMultiplication(const TInterpolationWeightMultiplication & _arg)
  {
    if (m_InterpolationWeightMultiplication != _arg)
    {
      m_InterpolationWeightMultiplication = _arg;
      this->Modified();
    }
  }

  TProjectedValueAccumulation &
  GetProjectedValueAccumulation()
  {
    return m_ProjectedValueAccumulation;
  }
  const TProjectedValueAccumulation &
  GetProjectedValueAccumulation() const
  {
    return m_ProjectedValueAccumulation;
  }
  void
  SetProjectedValueAccumulation(const TProjectedValueAccumulation & _arg)
  {
    if (m_ProjectedValueAccumulation != _arg)
    {
      m_ProjectedValueAccumulation = _arg;
      this->Modified();
    }
  }

  TSumAlongRay &
  GetSumAlongRay()
  {
    return m_SumAlongRay;
  }
  const TSumAlongRay &
  GetSumAlongRay() const
  {
    return m_SumAlongRay;
  }
  void
  SetSumAlongRay(const TSumAlongRay & _arg)
  {
    if (m_SumAlongRay != _arg)
    {
      m_SumAlongRay = _arg;
      this->Modified();
    }
  }

  void
  SetInferiorClipImage(const TClipImageType * inferiorClipImage)
  {
    // Process object is not const-correct so the const casting is required.
    this->SetInput("InferiorClipImage", const_cast<TClipImageType *>(inferiorClipImage));
  }
  typename TClipImageType::ConstPointer
  GetInferiorClipImage()
  {
    return static_cast<const TClipImageType *>(this->itk::ProcessObject::GetInput("InferiorClipImage"));
  }

  void
  SetSuperiorClipImage(const TClipImageType * superiorClipImage)
  {
    // Process object is not const-correct so the const casting is required.
    this->SetInput("SuperiorClipImage", const_cast<TClipImageType *>(superiorClipImage));
  }
  typename TClipImageType::ConstPointer
  GetSuperiorClipImage()
  {
    return static_cast<const TClipImageType *>(this->itk::ProcessObject::GetInput("SuperiorClipImage"));
  }

  itkGetMacro(InferiorClip, double);
  itkSetMacro(InferiorClip, double);
  itkGetMacro(SuperiorClip, double);
  itkSetMacro(SuperiorClip, double);

protected:
  JosephForwardProjectionImageFilter();
  ~JosephForwardProjectionImageFilter() override = default;

  void
  GenerateInputRequestedRegion() override;

  void
  ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, ThreadIdType threadId) override;

  void
  VerifyInputInformation() const override;

  inline OutputPixelType
  BilinearInterpolation(const ThreadIdType     threadId,
                        const double           stepLengthInVoxel,
                        const InputPixelType * pxiyi,
                        const InputPixelType * pxsyi,
                        const InputPixelType * pxiys,
                        const InputPixelType * pxsys,
                        const double           x,
                        const double           y,
                        const int              ox,
                        const int              oy);

  inline OutputPixelType
  BilinearInterpolationOnBorders(const ThreadIdType     threadId,
                                 const double           stepLengthInVoxel,
                                 const InputPixelType * pxiyi,
                                 const InputPixelType * pxsyi,
                                 const InputPixelType * pxiys,
                                 const InputPixelType * pxsys,
                                 const double           x,
                                 const double           y,
                                 const int              ox,
                                 const int              oy,
                                 const double           minx,
                                 const double           miny,
                                 const double           maxx,
                                 const double           maxy);

private:
  // Functors
  TInterpolationWeightMultiplication m_InterpolationWeightMultiplication;
  TProjectedValueAccumulation        m_ProjectedValueAccumulation;
  TSumAlongRay                       m_SumAlongRay;
  double                             m_InferiorClip{ 0. };
  double                             m_SuperiorClip{ 1. };
};

} // end namespace rtk

#ifndef ITK_MANUAL_INSTANTIATION
#  include "rtkJosephForwardProjectionImageFilter.hxx"
#endif

#endif