rtkJosephBackAttenuatedProjectionImageFilter.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 rtkJosephBackAttenuatedProjectionImageFilter_h
#define rtkJosephBackAttenuatedProjectionImageFilter_h

#include "rtkConfiguration.h"
#include "rtkJosephBackProjectionImageFilter.h"
#include "rtkThreeDCircularProjectionGeometry.h"

namespace rtk
{
namespace Functor
{
template <class TInput, class TCoordRepType, class TOutput = TInput>
class ITK_TEMPLATE_EXPORT InterpolationWeightMultiplicationAttenuatedBackProjection
{
public:
  InterpolationWeightMultiplicationAttenuatedBackProjection() { m_AttenuationPixel = 0; }

  ~InterpolationWeightMultiplicationAttenuatedBackProjection() = default;
  bool
  operator!=(const InterpolationWeightMultiplicationAttenuatedBackProjection &) const
  {
    return false;
  }

  bool
  operator==(const InterpolationWeightMultiplicationAttenuatedBackProjection & other) const
  {
    return !(*this != other);
  }

  inline TOutput
  operator()(const double stepLengthInVoxel, const TCoordRepType weight, const TInput * p, const int i)
  {
    const double w = weight * stepLengthInVoxel;

    m_AttenuationPixel += w * (p + m_AttenuationMinusEmissionMapsPtrDiff)[i];
    return w * (p + m_AttenuationMinusEmissionMapsPtrDiff)[i];
  }

  void
  SetAttenuationMinusEmissionMapsPtrDiff(std::ptrdiff_t pd)
  {
    m_AttenuationMinusEmissionMapsPtrDiff = pd;
  }
  TOutput *
  GetAttenuationPixel()
  {
    return &m_AttenuationPixel;
  }

private:
  std::ptrdiff_t m_AttenuationMinusEmissionMapsPtrDiff;
  TInput         m_AttenuationPixel;
};

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

  ComputeAttenuationCorrectionBackProjection() { m_Ex1 = 1; }

  ~ComputeAttenuationCorrectionBackProjection() = default;
  bool
  operator!=(const ComputeAttenuationCorrectionBackProjection &) const
  {
    return false;
  }

  bool
  operator==(const ComputeAttenuationCorrectionBackProjection & other) const
  {
    return !(*this != other);
  }

  inline TOutput
  operator()(const TInput rayValue, const TInput attenuationRay, const VectorType & stepInMM, bool & isNewRay)
  {
    if (isNewRay)
    {
      m_Ex1 = 1;
      isNewRay = false;
    }
    TInput ex2 = exp(-attenuationRay * stepInMM.GetNorm());
    TInput wf;
    if (*m_AttenuationPixel > 0)
    {
      wf = (m_Ex1 - ex2) / *m_AttenuationPixel;
    }
    else
    {
      wf = m_Ex1 * stepInMM.GetNorm();
    }

    m_Ex1 = ex2;
    *m_AttenuationPixel = 0;
    return wf * rayValue;
  }

  void
  SetAttenuationPixel(TInput * attenuationPixel)
  {
    m_AttenuationPixel = attenuationPixel;
  }

private:
  TInput   m_Ex1;
  TInput * m_AttenuationPixel;
};

template <class TInput, class TCoordRepType, class TOutput = TCoordRepType>
class ITK_TEMPLATE_EXPORT SplatWeightMultiplicationAttenuated
{
public:
  SplatWeightMultiplicationAttenuated() = default;
  ~SplatWeightMultiplicationAttenuated() = default;
  bool
  operator!=(const SplatWeightMultiplicationAttenuated &) const
  {
    return false;
  }

  bool
  operator==(const SplatWeightMultiplicationAttenuated & other) const
  {
    return !(*this != other);
  }

  inline void
  operator()(const TInput &      rayValue,
             TOutput &           output,
             const double        stepLengthInVoxel,
             const double        itkNotUsed(voxelSize),
             const TCoordRepType weight) const
  {
    output += rayValue * weight * stepLengthInVoxel;
  }
};
} // end namespace Functor

template <
  class TInputImage,
  class TOutputImage,
  class TInterpolationWeightMultiplication = Functor::InterpolationWeightMultiplicationAttenuatedBackProjection<
    typename TInputImage::PixelType,
    typename itk::PixelTraits<typename TInputImage::PixelType>::ValueType>,
  class TSplatWeightMultiplication = Functor::
    SplatWeightMultiplicationAttenuated<typename TInputImage::PixelType, double, typename TOutputImage::PixelType>,
  class TSumAlongRay = Functor::ComputeAttenuationCorrectionBackProjection<typename TInputImage::PixelType,
                                                                           typename TOutputImage::PixelType>>
class ITK_TEMPLATE_EXPORT JosephBackAttenuatedProjectionImageFilter
  : public JosephBackProjectionImageFilter<TInputImage,
                                           TOutputImage,
                                           TInterpolationWeightMultiplication,
                                           TSplatWeightMultiplication,
                                           TSumAlongRay>
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(JosephBackAttenuatedProjectionImageFilter);

  using Self = JosephBackAttenuatedProjectionImageFilter;
  using Superclass = JosephBackProjectionImageFilter<TInputImage,
                                                     TOutputImage,
                                                     TInterpolationWeightMultiplication,
                                                     TSplatWeightMultiplication,
                                                     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 GeometryType = rtk::ThreeDCircularProjectionGeometry;
  using GeometryPointer = typename GeometryType::Pointer;

  static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension;

  itkNewMacro(Self);

#ifdef itkOverrideGetNameOfClassMacro
  itkOverrideGetNameOfClassMacro(JosephBackAttenuatedProjectionImageFilter);
#else
  itkTypeMacro(JosephBackAttenuatedProjectionImageFilter, JosephBackProjectionImageFilter);
#endif


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

  void
  GenerateInputRequestedRegion() override;

  void
  VerifyInputInformation() const override;

  void
  GenerateData() override;

  void
  Init();
};
} // end namespace rtk

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

#endif