rtkMechlemOneStepSpectralReconstructionFilter.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 rtkMechlemOneStepSpectralReconstructionFilter_h
#define rtkMechlemOneStepSpectralReconstructionFilter_h

#include "rtkIterativeConeBeamReconstructionFilter.h"
#include "rtkThreeDCircularProjectionGeometry.h"
#include "rtkWeidingerForwardModelImageFilter.h"
#include "rtkGetNewtonUpdateImageFilter.h"
#include "rtkConstantImageSource.h"
#include "rtkNesterovUpdateImageFilter.h"
#include "rtkSeparableQuadraticSurrogateRegularizationImageFilter.h"
#include "rtkAddMatrixAndDiagonalImageFilter.h"
#include "rtkReorderProjectionsImageFilter.h"

#include <itkExtractImageFilter.h>
#include <itkAddImageFilter.h>
#include <itkMultiplyImageFilter.h>

#ifdef RTK_USE_CUDA
#  include "rtkCudaWeidingerForwardModelImageFilter.h"
#endif

namespace rtk
{
template <typename TOutputImage, typename TPhotonCounts, typename TSpectrum>
class ITK_TEMPLATE_EXPORT MechlemOneStepSpectralReconstructionFilter
  : public rtk::IterativeConeBeamReconstructionFilter<TOutputImage, TOutputImage>
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(MechlemOneStepSpectralReconstructionFilter);

  using Self = MechlemOneStepSpectralReconstructionFilter;
  using Superclass = IterativeConeBeamReconstructionFilter<TOutputImage, TOutputImage>;
  using Pointer = itk::SmartPointer<Self>;

  itkNewMacro(Self);

  itkOverrideGetNameOfClassMacro(MechlemOneStepSpectralReconstructionFilter);

  static constexpr unsigned int nBins = TPhotonCounts::PixelType::Dimension;
  static constexpr unsigned int nMaterials = TOutputImage::PixelType::Dimension;
  using dataType = typename TOutputImage::PixelType::ValueType;

  using CPUOutputImageType = typename itk::Image<typename TOutputImage::PixelType, TOutputImage::ImageDimension>;
#ifdef RTK_USE_CUDA
  typedef
    typename std::conditional<std::is_same<TOutputImage, CPUOutputImageType>::value,
                              itk::Image<itk::Vector<dataType, nMaterials * nMaterials>, TOutputImage::ImageDimension>,
                              itk::CudaImage<itk::Vector<dataType, nMaterials * nMaterials>,
                                             TOutputImage::ImageDimension>>::type HessiansImageType;
  typedef
    typename std::conditional<std::is_same<TOutputImage, CPUOutputImageType>::value,
                              itk::Image<dataType, TOutputImage::ImageDimension>,
                              itk::CudaImage<dataType, TOutputImage::ImageDimension>>::type SingleComponentImageType;
#else
  using HessiansImageType =
    typename itk::Image<itk::Vector<dataType, nMaterials * nMaterials>, TOutputImage::ImageDimension>;
  using SingleComponentImageType = typename itk::Image<dataType, TOutputImage::ImageDimension>;
#endif

#if !defined(ITK_WRAPPING_PARSER)
#  ifdef RTK_USE_CUDA
  typedef typename std::conditional<std::is_same<TOutputImage, CPUOutputImageType>::value,
                                    WeidingerForwardModelImageFilter<TOutputImage, TPhotonCounts, TSpectrum>,
                                    CudaWeidingerForwardModelImageFilter<TOutputImage, TPhotonCounts, TSpectrum>>::type
    WeidingerForwardModelType;
  typedef
    typename std::conditional<std::is_same<TOutputImage, CPUOutputImageType>::value,
                              JosephForwardProjectionImageFilter<SingleComponentImageType, SingleComponentImageType>,
                              CudaForwardProjectionImageFilter<SingleComponentImageType, SingleComponentImageType>>::
      type CudaSingleComponentForwardProjectionImageFilterType;
  typedef typename std::conditional<std::is_same<TOutputImage, CPUOutputImageType>::value,
                                    BackProjectionImageFilter<HessiansImageType, HessiansImageType>,
                                    CudaBackProjectionImageFilter<HessiansImageType>>::type
    CudaHessiansBackProjectionImageFilterType;
#  else
  using WeidingerForwardModelType = WeidingerForwardModelImageFilter<TOutputImage, TPhotonCounts, TSpectrum>;
  using CudaSingleComponentForwardProjectionImageFilterType =
    JosephForwardProjectionImageFilter<SingleComponentImageType, SingleComponentImageType>;
  using CudaHessiansBackProjectionImageFilterType = BackProjectionImageFilter<HessiansImageType, HessiansImageType>;
#  endif
  using GradientsImageType = TOutputImage;
#endif

  using ForwardProjectionType = typename Superclass::ForwardProjectionType;
  using BackProjectionType = typename Superclass::BackProjectionType;

#if !defined(ITK_WRAPPING_PARSER)

  using ExtractPhotonCountsFilterType = itk::ExtractImageFilter<TPhotonCounts, TPhotonCounts>;
  using AddFilterType = itk::AddImageFilter<GradientsImageType>;
  using SingleComponentForwardProjectionFilterType =
    rtk::ForwardProjectionImageFilter<SingleComponentImageType, SingleComponentImageType>;
  using ForwardProjectionFilterType = rtk::ForwardProjectionImageFilter<TOutputImage, TOutputImage>;
  using GradientsBackProjectionFilterType = rtk::BackProjectionImageFilter<GradientsImageType, GradientsImageType>;
  using HessiansBackProjectionFilterType = rtk::BackProjectionImageFilter<HessiansImageType, HessiansImageType>;
  using NesterovFilterType = rtk::NesterovUpdateImageFilter<TOutputImage>;
  using SingleComponentImageSourceType = rtk::ConstantImageSource<SingleComponentImageType>;
  using MaterialProjectionsSourceType = rtk::ConstantImageSource<TOutputImage>;
  using GradientsSourceType = rtk::ConstantImageSource<GradientsImageType>;
  using HessiansSourceType = rtk::ConstantImageSource<HessiansImageType>;
  using SQSRegularizationType = rtk::SeparableQuadraticSurrogateRegularizationImageFilter<GradientsImageType>;
  using AddMatrixAndDiagonalFilterType = rtk::AddMatrixAndDiagonalImageFilter<GradientsImageType, HessiansImageType>;
  using NewtonFilterType = rtk::GetNewtonUpdateImageFilter<GradientsImageType, HessiansImageType>;
  using MultiplyFilterType = itk::MultiplyImageFilter<TOutputImage, SingleComponentImageType>;
  using MultiplyGradientFilterType = itk::MultiplyImageFilter<GradientsImageType, SingleComponentImageType>;
  using ReorderProjectionsFilterPhotonCountsType = rtk::ReorderProjectionsImageFilter<TPhotonCounts>;
  using ReorderProjectionsFilterProjectionsWeightsType = rtk::ReorderProjectionsImageFilter<SingleComponentImageType>;
#endif

  itkSetConstObjectMacro(Geometry, ThreeDCircularProjectionGeometry);

  itkSetMacro(NumberOfIterations, int);
  itkGetMacro(NumberOfIterations, int);

  itkSetMacro(NumberOfSubsets, int);
  itkGetMacro(NumberOfSubsets, int);

  itkSetMacro(ResetNesterovEvery, int);
  itkGetMacro(ResetNesterovEvery, int);

  void
  SetInputMaterialVolumes(const TOutputImage * materialVolumes);
  void
  SetInputPhotonCounts(const TPhotonCounts * photonCounts);
  void
  SetInputSpectrum(const TSpectrum * spectrum);
  void
  SetSupportMask(const SingleComponentImageType * support);
  void
  SetSpatialRegularizationWeights(const SingleComponentImageType * regweights);
  void
  SetProjectionWeights(const SingleComponentImageType * weiprojections);

  itkSetMacro(RegularizationWeights, typename TOutputImage::PixelType);
  itkGetMacro(RegularizationWeights, typename TOutputImage::PixelType);

  itkSetMacro(RegularizationRadius, typename TOutputImage::RegionType::SizeType);
  itkGetMacro(RegularizationRadius, typename TOutputImage::RegionType::SizeType);

  using BinnedDetectorResponseType = vnl_matrix<dataType>;
  using MaterialAttenuationsType = vnl_matrix<dataType>;
  virtual void
  SetBinnedDetectorResponse(const BinnedDetectorResponseType & detResp);
  virtual void
  SetMaterialAttenuations(const MaterialAttenuationsType & matAtt);

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

  void
  VerifyPreconditions() const override;

  void
  GenerateData() override;

#if !defined(ITK_WRAPPING_PARSER)

  typename ExtractPhotonCountsFilterType::Pointer                  m_ExtractPhotonCountsFilter;
  typename AddFilterType::Pointer                                  m_AddGradients;
  typename SingleComponentForwardProjectionFilterType::Pointer     m_SingleComponentForwardProjectionFilter;
  typename MaterialProjectionsSourceType::Pointer                  m_ProjectionsSource;
  typename SingleComponentImageSourceType::Pointer                 m_SingleComponentProjectionsSource;
  typename SingleComponentImageSourceType::Pointer                 m_SingleComponentVolumeSource;
  typename GradientsSourceType::Pointer                            m_GradientsSource;
  typename HessiansSourceType::Pointer                             m_HessiansSource;
  typename WeidingerForwardModelType::Pointer                      m_WeidingerForward;
  typename SQSRegularizationType::Pointer                          m_SQSRegul;
  typename AddMatrixAndDiagonalFilterType::Pointer                 m_AddHessians;
  typename NewtonFilterType::Pointer                               m_NewtonFilter;
  typename NesterovFilterType::Pointer                             m_NesterovFilter;
  typename ForwardProjectionFilterType::Pointer                    m_ForwardProjectionFilter;
  typename GradientsBackProjectionFilterType::Pointer              m_GradientsBackProjectionFilter;
  typename HessiansBackProjectionFilterType::Pointer               m_HessiansBackProjectionFilter;
  typename MultiplyFilterType::Pointer                             m_MultiplySupportFilter;
  typename MultiplyGradientFilterType::Pointer                     m_MultiplyRegulGradientsFilter;
  typename MultiplyGradientFilterType::Pointer                     m_MultiplyRegulHessiansFilter;
  typename MultiplyGradientFilterType::Pointer                     m_MultiplyGradientToBeBackprojectedFilter;
  typename ReorderProjectionsFilterPhotonCountsType::Pointer       m_ReorderPhotonCountsFilter;
  typename ReorderProjectionsFilterProjectionsWeightsType::Pointer m_ReorderProjectionsWeightsFilter;
#endif

  void
  VerifyInputInformation() const override
  {}

  void
  GenerateInputRequestedRegion() override;
  void
  GenerateOutputInformation() override;

  typename TOutputImage::ConstPointer
  GetInputMaterialVolumes();
  typename TPhotonCounts::ConstPointer
  GetInputPhotonCounts();
  typename TSpectrum::ConstPointer
  GetInputSpectrum();
  typename SingleComponentImageType::ConstPointer
  GetSupportMask();
  typename SingleComponentImageType::ConstPointer
  GetSpatialRegularizationWeights();
  typename SingleComponentImageType::ConstPointer
  GetProjectionWeights();

#if !defined(ITK_WRAPPING_PARSER)

  typename SingleComponentForwardProjectionFilterType::Pointer
  InstantiateSingleComponentForwardProjectionFilter(int fwtype);
  typename HessiansBackProjectionFilterType::Pointer
  InstantiateHessiansBackProjectionFilter(int bptype);
#endif


  ThreeDCircularProjectionGeometry::ConstPointer m_Geometry;

  int              m_NumberOfIterations;
  int              m_NumberOfProjectionsPerSubset;
  int              m_NumberOfSubsets;
  std::vector<int> m_NumberOfProjectionsInSubset;
  int              m_NumberOfProjections;
  int              m_ResetNesterovEvery;

  typename TOutputImage::PixelType            m_RegularizationWeights;
  typename TOutputImage::RegionType::SizeType m_RegularizationRadius;
};
} // namespace rtk


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

#endif