rtkFourDROOSTERConeBeamReconstructionFilter.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 rtkFourDROOSTERConeBeamReconstructionFilter_h
#define rtkFourDROOSTERConeBeamReconstructionFilter_h

#include "rtkFourDConjugateGradientConeBeamReconstructionFilter.h"
#include "rtkTotalVariationDenoiseSequenceImageFilter.h"
#ifdef RTK_USE_CUDA
#  include "rtkCudaLastDimensionTVDenoisingImageFilter.h"
#  include "rtkCudaAverageOutOfROIImageFilter.h"
#else
#  include "rtkTotalVariationDenoisingBPDQImageFilter.h"
#  include "rtkAverageOutOfROIImageFilter.h"
#endif
#include "rtkDaubechiesWaveletsDenoiseSequenceImageFilter.h"
#include "rtkWarpSequenceImageFilter.h"
#include "rtkUnwarpSequenceImageFilter.h"
#include "rtkLastDimensionL0GradientDenoisingImageFilter.h"
#include "rtkTotalNuclearVariationDenoisingBPDQImageFilter.h"

#include <itkThresholdImageFilter.h>
#include <itkSubtractImageFilter.h>
#include <itkAddImageFilter.h>

#include <itkResampleImageFilter.h>
#include <itkNearestNeighborInterpolateImageFunction.h>
#include <itkIdentityTransform.h>


namespace rtk
{
template <typename VolumeSeriesType, typename ProjectionStackType>
class ITK_TEMPLATE_EXPORT FourDROOSTERConeBeamReconstructionFilter
  : public rtk::IterativeConeBeamReconstructionFilter<VolumeSeriesType, ProjectionStackType>
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(FourDROOSTERConeBeamReconstructionFilter);

  using Self = FourDROOSTERConeBeamReconstructionFilter;
  using Superclass = rtk::IterativeConeBeamReconstructionFilter<VolumeSeriesType, ProjectionStackType>;
  using Pointer = itk::SmartPointer<Self>;
  using VolumeType = ProjectionStackType;
  using CovariantVectorForSpatialGradient =
    itk::CovariantVector<typename VolumeSeriesType::ValueType, VolumeSeriesType::ImageDimension - 1>;
  using CovariantVectorForTemporalGradient = itk::CovariantVector<typename VolumeSeriesType::ValueType, 1>;
  using DVFVectorType = CovariantVectorForSpatialGradient;

  using CPUVolumeSeriesType =
    typename itk::Image<typename VolumeSeriesType::PixelType, VolumeSeriesType::ImageDimension>;
#ifdef RTK_USE_CUDA
  typedef
    typename std::conditional<std::is_same<VolumeSeriesType, CPUVolumeSeriesType>::value,
                              itk::Image<CovariantVectorForSpatialGradient, VolumeSeriesType::ImageDimension>,
                              itk::CudaImage<CovariantVectorForSpatialGradient, VolumeSeriesType::ImageDimension>>::type
                                                                                            SpatialGradientImageType;
  typedef typename std::conditional<std::is_same<VolumeSeriesType, CPUVolumeSeriesType>::value,
                                    itk::Image<CovariantVectorForTemporalGradient, VolumeSeriesType::ImageDimension>,
                                    itk::CudaImage<CovariantVectorForTemporalGradient,
                                                   VolumeSeriesType::ImageDimension>>::type TemporalGradientImageType;
  typedef typename std::conditional<std::is_same<VolumeSeriesType, CPUVolumeSeriesType>::value,
                                    itk::Image<DVFVectorType, VolumeSeriesType::ImageDimension>,
                                    itk::CudaImage<DVFVectorType, VolumeSeriesType::ImageDimension>>::type
    DVFSequenceImageType;
  typedef
    typename std::conditional<std::is_same<VolumeSeriesType, CPUVolumeSeriesType>::value,
                              itk::Image<DVFVectorType, VolumeSeriesType::ImageDimension - 1>,
                              itk::CudaImage<DVFVectorType, VolumeSeriesType::ImageDimension - 1>>::type DVFImageType;
  typedef typename std::conditional<std::is_same<VolumeSeriesType, CPUVolumeSeriesType>::value,
                                    AverageOutOfROIImageFilter<VolumeSeriesType, VolumeType>,
                                    CudaAverageOutOfROIImageFilter>::type          AverageOutOfROIFilterType;
  typedef typename std::conditional<std::is_same<VolumeSeriesType, CPUVolumeSeriesType>::value,
                                    TotalVariationDenoisingBPDQImageFilter<VolumeSeriesType, TemporalGradientImageType>,
                                    CudaLastDimensionTVDenoisingImageFilter>::type TemporalTVDenoisingFilterType;
#else
  using SpatialGradientImageType = itk::Image<CovariantVectorForSpatialGradient, VolumeSeriesType::ImageDimension>;
  using TemporalGradientImageType = itk::Image<CovariantVectorForTemporalGradient, VolumeSeriesType::ImageDimension>;
  using DVFSequenceImageType = itk::Image<DVFVectorType, VolumeSeriesType::ImageDimension>;
  using DVFImageType = itk::Image<DVFVectorType, VolumeSeriesType::ImageDimension - 1>;
  using AverageOutOfROIFilterType = AverageOutOfROIImageFilter<VolumeSeriesType, VolumeType>;
  using TemporalTVDenoisingFilterType =
    TotalVariationDenoisingBPDQImageFilter<VolumeSeriesType, TemporalGradientImageType>;
#endif

  itkNewMacro(Self);

  itkOverrideGetNameOfClassMacro(FourDROOSTERConeBeamReconstructionFilter);

  void
  SetInputVolumeSeries(const VolumeSeriesType * VolumeSeries);
  typename VolumeSeriesType::ConstPointer
  GetInputVolumeSeries();

  void
  SetInputProjectionStack(const ProjectionStackType * Projection);
  typename ProjectionStackType::Pointer
  GetInputProjectionStack();

  void
  SetMotionMask(const VolumeType * mask);
  typename VolumeType::Pointer
  GetMotionMask();

  void
  SetDisplacementField(const DVFSequenceImageType * DVFs);
  void
  SetInverseDisplacementField(const DVFSequenceImageType * DVFs);
  typename DVFSequenceImageType::Pointer
  GetDisplacementField();
  typename DVFSequenceImageType::Pointer
  GetInverseDisplacementField();

  using FourDCGFilterType =
    rtk::FourDConjugateGradientConeBeamReconstructionFilter<VolumeSeriesType, ProjectionStackType>;
  using ThresholdFilterType = itk::ThresholdImageFilter<VolumeSeriesType>;
  using ResampleFilterType = itk::ResampleImageFilter<VolumeType, VolumeType>;
  using SpatialTVDenoisingFilterType = rtk::TotalVariationDenoiseSequenceImageFilter<VolumeSeriesType>;
  using SpatialWaveletsDenoisingFilterType = rtk::DaubechiesWaveletsDenoiseSequenceImageFilter<VolumeSeriesType>;
  using WarpSequenceFilterType =
    rtk::WarpSequenceImageFilter<VolumeSeriesType, DVFSequenceImageType, VolumeType, DVFImageType>;
  using UnwarpSequenceFilterType =
    rtk::UnwarpSequenceImageFilter<VolumeSeriesType, DVFSequenceImageType, VolumeType, DVFImageType>;
  using SubtractFilterType = itk::SubtractImageFilter<VolumeSeriesType, VolumeSeriesType>;
  using AddFilterType = itk::AddImageFilter<VolumeSeriesType, VolumeSeriesType>;
  using TemporalL0DenoisingFilterType = rtk::LastDimensionL0GradientDenoisingImageFilter<VolumeSeriesType>;
  using TNVDenoisingFilterType =
    rtk::TotalNuclearVariationDenoisingBPDQImageFilter<VolumeSeriesType, SpatialGradientImageType>;

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

  virtual void
  SetWeights(const itk::Array2D<float> _arg);

  itkSetMacro(DisableDisplacedDetectorFilter, bool);
  itkGetMacro(DisableDisplacedDetectorFilter, bool);

  // Regularization steps to perform
  itkSetMacro(PerformPositivity, bool);
  itkGetMacro(PerformPositivity, bool);
  itkSetMacro(PerformMotionMask, bool);
  itkGetMacro(PerformMotionMask, bool);
  itkSetMacro(PerformTVSpatialDenoising, bool);
  itkGetMacro(PerformTVSpatialDenoising, bool);
  itkSetMacro(PerformWaveletsSpatialDenoising, bool);
  itkGetMacro(PerformWaveletsSpatialDenoising, bool);
  itkSetMacro(PerformWarping, bool);
  itkGetMacro(PerformWarping, bool);
  itkSetMacro(PerformTVTemporalDenoising, bool);
  itkGetMacro(PerformTVTemporalDenoising, bool);
  itkSetMacro(PerformL0TemporalDenoising, bool);
  itkGetMacro(PerformL0TemporalDenoising, bool);
  itkSetMacro(PerformTNVDenoising, bool);
  itkGetMacro(PerformTNVDenoising, bool);
  itkSetMacro(ComputeInverseWarpingByConjugateGradient, bool);
  itkGetMacro(ComputeInverseWarpingByConjugateGradient, bool);
  itkSetMacro(UseNearestNeighborInterpolationInWarping, bool);
  itkGetMacro(UseNearestNeighborInterpolationInWarping, bool);
  itkGetMacro(CudaConjugateGradient, bool);
  itkSetMacro(CudaConjugateGradient, bool);

  itkSetMacro(UseCudaCyclicDeformation, bool);
  itkGetMacro(UseCudaCyclicDeformation, bool);

  // Regularization parameters
  itkSetMacro(GammaTVSpace, float);
  itkGetMacro(GammaTVSpace, float);
  itkSetMacro(GammaTVTime, float);
  itkGetMacro(GammaTVTime, float);
  itkSetMacro(GammaTNV, float);
  itkGetMacro(GammaTNV, float);
  itkSetMacro(LambdaL0Time, float);
  itkGetMacro(LambdaL0Time, float);
  itkSetMacro(SoftThresholdWavelets, float);
  itkGetMacro(SoftThresholdWavelets, float);
  itkSetMacro(PhaseShift, float);
  itkGetMacro(PhaseShift, float);

  itkGetMacro(NumberOfLevels, unsigned int);
  itkSetMacro(NumberOfLevels, unsigned int);

  itkGetMacro(Order, unsigned int);
  itkSetMacro(Order, unsigned int);

  // Iterations
  itkSetMacro(MainLoop_iterations, int);
  itkGetMacro(MainLoop_iterations, int);
  itkSetMacro(CG_iterations, int);
  itkGetMacro(CG_iterations, int);
  itkSetMacro(TV_iterations, int);
  itkGetMacro(TV_iterations, int);
  itkSetMacro(L0_iterations, int);
  itkGetMacro(L0_iterations, int);

  // Geometry
  itkSetConstObjectMacro(Geometry, ThreeDCircularProjectionGeometry);
  itkGetConstObjectMacro(Geometry, ThreeDCircularProjectionGeometry);

  virtual void
  SetSignal(const std::vector<double> signal);

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

  void
  VerifyPreconditions() const override;

  void
  GenerateData() override;

  void
  GenerateOutputInformation() override;

  void
  GenerateInputRequestedRegion() override;

  // Inputs are not supposed to occupy the same physical space,
  // so there is nothing to verify
  void
  VerifyInputInformation() const override
  {}

  typename FourDCGFilterType::Pointer                  m_FourDCGFilter;
  typename ThresholdFilterType::Pointer                m_PositivityFilter;
  typename ResampleFilterType::Pointer                 m_ResampleFilter;
  typename AverageOutOfROIFilterType::Pointer          m_AverageOutOfROIFilter;
  typename SpatialTVDenoisingFilterType::Pointer       m_TVDenoisingSpace;
  typename SpatialWaveletsDenoisingFilterType::Pointer m_WaveletsDenoisingSpace;
  typename WarpSequenceFilterType::Pointer             m_Warp;
  typename TemporalTVDenoisingFilterType::Pointer      m_TVDenoisingTime;
  typename UnwarpSequenceFilterType::Pointer           m_Unwarp;
  typename WarpSequenceFilterType::Pointer             m_InverseWarp;
  typename SubtractFilterType::Pointer                 m_SubtractFilter;
  typename AddFilterType::Pointer                      m_AddFilter;
  typename TemporalL0DenoisingFilterType::Pointer      m_L0DenoisingTime;
  typename TNVDenoisingFilterType::Pointer             m_TNVDenoising;

  // Booleans :
  // should warping be performed ?
  // should conjugate gradient be performed on GPU ?
  // should wavelets replace TV in spatial denoising ?
  bool m_PerformPositivity;
  bool m_PerformMotionMask;
  bool m_PerformTVSpatialDenoising;
  bool m_PerformWaveletsSpatialDenoising;
  bool m_PerformWarping;
  bool m_PerformTVTemporalDenoising;
  bool m_PerformL0TemporalDenoising;
  bool m_PerformTNVDenoising;
  bool m_ComputeInverseWarpingByConjugateGradient;
  bool m_UseNearestNeighborInterpolationInWarping; // Default is false, linear interpolation is used instead
  bool m_CudaConjugateGradient;
  bool m_UseCudaCyclicDeformation;
  bool m_DisableDisplacedDetectorFilter;

  // Regularization parameters
  float m_GammaTVSpace;
  float m_GammaTVTime;
  float m_GammaTNV;
  float m_LambdaL0Time;
  float m_SoftThresholdWavelets;
  float m_PhaseShift;
  bool  m_DimensionsProcessedForTVSpace[VolumeSeriesType::ImageDimension];
  bool  m_DimensionsProcessedForTVTime[VolumeSeriesType::ImageDimension];

  typename itk::ImageToImageFilter<VolumeSeriesType, VolumeSeriesType>::Pointer m_DownstreamFilter;

  unsigned int m_Order;
  unsigned int m_NumberOfLevels;

  // Iterations
  int m_MainLoop_iterations;
  int m_CG_iterations;
  int m_TV_iterations;
  int m_L0_iterations;

  // Geometry
  typename rtk::ThreeDCircularProjectionGeometry::ConstPointer m_Geometry;

  // Signal
  std::vector<double> m_Signal;
};
} // namespace rtk


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

#endif