rtkRegularizedConjugateGradientConeBeamReconstructionFilter.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 rtkRegularizedConjugateGradientConeBeamReconstructionFilter_h
#define rtkRegularizedConjugateGradientConeBeamReconstructionFilter_h

#include "rtkConjugateGradientConeBeamReconstructionFilter.h"
#ifdef RTK_USE_CUDA
#  include "rtkCudaTotalVariationDenoisingBPDQImageFilter.h"
#else
#  include "rtkTotalVariationDenoisingBPDQImageFilter.h"
#endif
#include "rtkDeconstructSoftThresholdReconstructImageFilter.h"

#include <itkThresholdImageFilter.h>

namespace rtk
{
template <typename TImage>
class ITK_TEMPLATE_EXPORT RegularizedConjugateGradientConeBeamReconstructionFilter
  : public rtk::IterativeConeBeamReconstructionFilter<TImage, TImage>
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(RegularizedConjugateGradientConeBeamReconstructionFilter);

  using Self = RegularizedConjugateGradientConeBeamReconstructionFilter;
  using Superclass = rtk::IterativeConeBeamReconstructionFilter<TImage, TImage>;
  using Pointer = itk::SmartPointer<Self>;
  using CovariantVectorForSpatialGradient = itk::CovariantVector<typename TImage::ValueType, TImage::ImageDimension>;

#ifdef RTK_USE_CUDA
  using GradientImageType = itk::CudaImage<CovariantVectorForSpatialGradient, TImage::ImageDimension>;
#else
  using GradientImageType = itk::Image<CovariantVectorForSpatialGradient, TImage::ImageDimension>;
#endif

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

  itkNewMacro(Self);

#ifdef itkOverrideGetNameOfClassMacro
  itkOverrideGetNameOfClassMacro(RegularizedConjugateGradientConeBeamReconstructionFilter);
#else
  itkTypeMacro(RegularizedConjugateGradientConeBeamReconstructionFilter, itk::ImageToImageFilter);
#endif


  void
  SetInputVolume(const TImage * Volume);
  typename TImage::ConstPointer
  GetInputVolume();

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

  void
  SetInputWeights(const TImage * Weights);
  typename TImage::Pointer
  GetInputWeights();

  void
  SetSupportMask(const TImage * SupportMask);
  typename TImage::ConstPointer
  GetSupportMask();

  using CGFilterType = rtk::ConjugateGradientConeBeamReconstructionFilter<TImage>;
  using ThresholdFilterType = itk::ThresholdImageFilter<TImage>;
  using TVDenoisingFilterType = rtk::TotalVariationDenoisingBPDQImageFilter<TImage, GradientImageType>;
  using WaveletsDenoisingFilterType = rtk::DeconstructSoftThresholdReconstructImageFilter<TImage>;
  using SoftThresholdFilterType = rtk::SoftThresholdImageFilter<TImage, TImage>;

  // Regularization steps to perform
  itkSetMacro(PerformPositivity, bool);
  itkGetMacro(PerformPositivity, bool);
  itkSetMacro(PerformTVSpatialDenoising, bool);
  itkGetMacro(PerformTVSpatialDenoising, bool);
  itkSetMacro(PerformWaveletsSpatialDenoising, bool);
  itkGetMacro(PerformWaveletsSpatialDenoising, bool);
  itkSetMacro(PerformSoftThresholdOnImage, bool);
  itkGetMacro(PerformSoftThresholdOnImage, bool);

  // Regularization parameters
  itkSetMacro(GammaTV, float);
  itkGetMacro(GammaTV, float);
  itkSetMacro(SoftThresholdWavelets, float);
  itkGetMacro(SoftThresholdWavelets, float);
  itkSetMacro(SoftThresholdOnImage, float);
  itkGetMacro(SoftThresholdOnImage, 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);

  // Geometry
  itkSetObjectMacro(Geometry, ThreeDCircularProjectionGeometry);
  itkGetModifiableObjectMacro(Geometry, ThreeDCircularProjectionGeometry);

  itkSetMacro(Preconditioned, bool);
  itkGetMacro(Preconditioned, bool);

  itkSetMacro(Tikhonov, float);
  itkGetMacro(Tikhonov, float);
  itkSetMacro(Gamma, float);
  itkGetMacro(Gamma, float);

  itkSetMacro(CudaConjugateGradient, bool);
  itkGetMacro(CudaConjugateGradient, bool);

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

protected:
  RegularizedConjugateGradientConeBeamReconstructionFilter();
  ~RegularizedConjugateGradientConeBeamReconstructionFilter() 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 CGFilterType::Pointer                m_CGFilter;
  typename ThresholdFilterType::Pointer         m_PositivityFilter;
  typename TVDenoisingFilterType::Pointer       m_TVDenoising;
  typename WaveletsDenoisingFilterType::Pointer m_WaveletsDenoising;
  typename SoftThresholdFilterType::Pointer     m_SoftThresholdFilter;

  // Booleans for each regularization (should it be performed or not)
  // as well as to choose whether CG should be on GPU or not
  bool m_PerformPositivity;
  bool m_PerformTVSpatialDenoising;
  bool m_PerformWaveletsSpatialDenoising;
  bool m_CudaConjugateGradient;
  bool m_PerformSoftThresholdOnImage;

  // Regularization parameters
  float m_GammaTV;
  float m_Gamma;
  float m_Tikhonov;
  float m_SoftThresholdWavelets;
  float m_SoftThresholdOnImage;
  bool  m_DimensionsProcessedForTV[TImage::ImageDimension];
  bool  m_Preconditioned;
  bool  m_RegularizedCG;

  unsigned int m_Order;
  unsigned int m_NumberOfLevels;

  bool m_DisableDisplacedDetectorFilter;

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

  // Geometry
  typename rtk::ThreeDCircularProjectionGeometry::Pointer m_Geometry;
};
} // namespace rtk


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

#endif