itkDiffusionTensor3DReconstructionImageFilter.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  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 itkDiffusionTensor3DReconstructionImageFilter_h
#define itkDiffusionTensor3DReconstructionImageFilter_h
 
#include "itkImageToImageFilter.h"
#include "itkSpatialObject.h"
#include "itkDiffusionTensor3D.h"
#include "vnl/vnl_matrix.h"
#include "vnl/vnl_vector_fixed.h"
#include "vnl/vnl_matrix_fixed.h"
#include "itkMathSVD.h"
#if !defined(ITK_LEGACY_REMOVE) && !defined(ITK_FUTURE_LEGACY_REMOVE)
#  include "vnl/algo/vnl_svd.h" // transitional transitive include; dropped on ITK legacy removal
#endif
#include "itkVectorContainer.h"
#include "itkVectorImage.h"
#include "ITKDiffusionTensorImageExport.h"
 
namespace itk
{
class DiffusionTensor3DReconstructionImageFilterEnums
{
public:
  enum class GradientImageFormat : uint8_t
  {
    GradientIsInASingleImage = 1,
    GradientIsInManyImages,
    Else
  };
};
// Define how to print enumeration
extern ITKDiffusionTensorImage_EXPORT std::ostream &
operator<<(std::ostream & out, const DiffusionTensor3DReconstructionImageFilterEnums::GradientImageFormat value);

 
template <typename TReferenceImagePixelType,
          typename TGradientImagePixelType = TReferenceImagePixelType,
          typename TTensorPixelType = double,
          typename TMaskImageType = Image<unsigned char, 3>>
class ITK_TEMPLATE_EXPORT DiffusionTensor3DReconstructionImageFilter
  : public ImageToImageFilter<Image<TReferenceImagePixelType, 3>, Image<DiffusionTensor3D<TTensorPixelType>, 3>>
{
public:
  using Self = DiffusionTensor3DReconstructionImageFilter;
  using Pointer = SmartPointer<Self>;
  using ConstPointer = SmartPointer<const Self>;
  using Superclass =
    ImageToImageFilter<Image<TReferenceImagePixelType, 3>, Image<DiffusionTensor3D<TTensorPixelType>, 3>>;

  itkNewMacro(Self);

  itkOverrideGetNameOfClassMacro(DiffusionTensor3DReconstructionImageFilter);
 
  using ReferencePixelType = TReferenceImagePixelType;
 
  using GradientPixelType = TGradientImagePixelType;
 
  using TensorPixelType = DiffusionTensor3D<TTensorPixelType>;

  using ReferenceImageType = typename Superclass::InputImageType;
 
  using TensorImageType = Image<TensorPixelType, 3>;
 
  using OutputImageType = TensorImageType;
 
  using typename Superclass::OutputImageRegionType;

  using GradientImageType = Image<GradientPixelType, 3>;

  using GradientImagesType = VectorImage<GradientPixelType, 3>;

  using MaskSpatialObjectType = SpatialObject<3>;

  using MaskImageType = TMaskImageType;

  using TensorBasisMatrixType = vnl_matrix_fixed<double, 6, 6>;
 
  using CoefficientMatrixType = vnl_matrix<double>;

  using GradientDirectionType = vnl_vector_fixed<double, 3>;

  using GradientDirectionContainerType = VectorContainer<unsigned int, GradientDirectionType>;

  void
  AddGradientImage(const GradientDirectionType &, const GradientImageType * gradientImage);
  const GradientImageType *
  GetGradientImage(unsigned int index) const;

  void
  SetGradientImage(GradientDirectionContainerType *, const GradientImagesType * gradientImage);

  void
  SetReferenceImage(ReferenceImageType * referenceImage)
  {
    if (m_GradientImageTypeEnumeration ==
        DiffusionTensor3DReconstructionImageFilterEnums::GradientImageFormat::GradientIsInASingleImage)
    {
      itkExceptionStringMacro(
        "Cannot call both methods:AddGradientImage and SetGradientImage. Please call only one of them.");
    }
 
    this->ProcessObject::SetNthInput(0, referenceImage);
 
    m_GradientImageTypeEnumeration =
      DiffusionTensor3DReconstructionImageFilterEnums::GradientImageFormat::GradientIsInManyImages;
  }

  virtual ReferenceImageType *
  GetReferenceImage()
  {
    return (static_cast<ReferenceImageType *>(this->ProcessObject::GetInput(0)));
  }

  [[nodiscard]] virtual GradientDirectionType
  GetGradientDirection(unsigned int idx) const
  {
    if (idx >= m_NumberOfGradientDirections)
    {
      itkExceptionMacro("Gradient direction " << idx << " does not exist");
    }
    return m_GradientDirectionContainer->ElementAt(idx);
  }

  void
  SetMaskImage(MaskImageType * maskImage);

  void
  SetMaskSpatialObject(MaskSpatialObjectType * maskSpatialObject);
 

  itkSetMacro(Threshold, ReferencePixelType);
  itkGetConstMacro(Threshold, ReferencePixelType);

  itkSetMacro(BValue, TTensorPixelType);
#ifdef GetBValue
#  undef GetBValue
#endif
  itkGetConstReferenceMacro(BValue, TTensorPixelType);
 
  itkConceptMacro(ReferenceEqualityComparableCheck, (Concept::EqualityComparable<ReferencePixelType>));
  itkConceptMacro(TensorEqualityComparableCheck, (Concept::EqualityComparable<TensorPixelType>));
  itkConceptMacro(GradientConvertibleToDoubleCheck, (Concept::Convertible<GradientPixelType, double>));
  itkConceptMacro(DoubleConvertibleToTensorCheck, (Concept::Convertible<double, TensorPixelType>));
  itkConceptMacro(GradientReferenceAdditiveOperatorsCheck,
                  (Concept::AdditiveOperators<GradientPixelType, GradientPixelType, ReferencePixelType>));
  itkConceptMacro(GradientReferenceAdditiveAndAssignOperatorsCheck,
                  (Concept::AdditiveAndAssignOperators<GradientPixelType, ReferencePixelType>));
 
  itkConceptMacro(ReferenceOStreamWritableCheck, (Concept::OStreamWritable<ReferencePixelType>));
  itkConceptMacro(TensorOStreamWritableCheck, (Concept::OStreamWritable<TensorPixelType>));
 
protected:
  DiffusionTensor3DReconstructionImageFilter();
  ~DiffusionTensor3DReconstructionImageFilter() override = default;
  void
  PrintSelf(std::ostream & os, Indent indent) const override;
 
  void
  ComputeTensorBasis();
 
  void
  BeforeThreadedGenerateData() override;
 
  void
  DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread) override;
 
  void
  VerifyPreconditions() const override;

  using GradientImageTypeEnumeration = DiffusionTensor3DReconstructionImageFilterEnums::GradientImageFormat;
#if !defined(ITK_LEGACY_REMOVE)
  // We need to expose the enum values at the class level
  // for backwards compatibility
  static constexpr GradientImageTypeEnumeration GradientIsInASingleImage =
    GradientImageTypeEnumeration::GradientIsInASingleImage;
  static constexpr GradientImageTypeEnumeration GradientIsInManyImages =
    GradientImageTypeEnumeration::GradientIsInManyImages;
  static constexpr GradientImageTypeEnumeration Else = GradientImageTypeEnumeration::Else;
#endif
 
private:
  /* Tensor basis coeffs */
  TensorBasisMatrixType m_TensorBasis{};
 
  /* Pseudo-inverse of m_TensorBasis (dual tensor basis), precomputed in BeforeThreadedGenerateData. */
  vnl_matrix<double> m_TensorBasisInverse{};
 
  CoefficientMatrixType m_BMatrix{};

  GradientDirectionContainerType::Pointer m_GradientDirectionContainer{};

  unsigned int m_NumberOfGradientDirections{};

  unsigned int m_NumberOfBaselineImages{};

  ReferencePixelType m_Threshold{};

  TTensorPixelType m_BValue{};

  GradientImageTypeEnumeration m_GradientImageTypeEnumeration{};

  bool m_MaskImagePresent{};
};
} // namespace itk
 
#ifndef ITK_MANUAL_INSTANTIATION
#  include "itkDiffusionTensor3DReconstructionImageFilter.hxx"
#endif
 
#endif