itkMRIBiasFieldCorrectionFilter.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
 *
 *         http://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 itkMRIBiasFieldCorrectionFilter_h
#define itkMRIBiasFieldCorrectionFilter_h
 
#include <ctime>
 
#include "itkImageToImageFilter.h"
#include "itkArray2D.h"
#include "itkMRASlabIdentifier.h"
#include "itkCompositeValleyFunction.h"
#include "itkMultivariateLegendrePolynomial.h"
#include "itkNormalVariateGenerator.h"
#include "itkOnePlusOneEvolutionaryOptimizer.h"
#include "itkImageRegionIterator.h"
 
namespace itk
{
template <typename TImage, typename TImageMask, typename TBiasField>
class ITK_TEMPLATE_EXPORT MRIBiasEnergyFunction : public SingleValuedCostFunction
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(MRIBiasEnergyFunction);
 
  using Self = MRIBiasEnergyFunction;
  using Superclass = SingleValuedCostFunction;
  using Pointer = SmartPointer<Self>;
  using ConstPointer = SmartPointer<const Self>;
 
  itkTypeMacro(MRIBiasEnergyFunction, SingleValuedCostFunction);
 
  itkNewMacro(Self);
 
  using ImageType = TImage;
  using ImagePointer = typename ImageType::Pointer;
  using ImageElementType = typename ImageType::PixelType;
  using ImageIndexType = typename ImageType::IndexType;
  using ImageRegionType = typename ImageType::RegionType;
  using MaskType = TImageMask;
  using MaskPointer = typename MaskType::Pointer;
  using MaskElementType = typename MaskType::PixelType;
 
  using BiasFieldType = TBiasField;
 
  using ParametersType = typename Superclass::ParametersType;
 
  using DerivativeType = Superclass::DerivativeType;
 
  using MeasureType = Superclass::MeasureType;
 
  static constexpr unsigned int SpaceDimension = 3;
 
  using InternalEnergyFunction = CompositeValleyFunction;
 
  using SamplingFactorType = unsigned int[SpaceDimension];
 
  itkSetObjectMacro(Image, ImageType);
 
  itkSetObjectMacro(Mask, MaskType);
 
  itkSetMacro(Region, ImageRegionType);
 
  void
  SetBiasField(BiasFieldType * bias)
  {
    m_BiasField = bias;
  }
 
  void
  SetSamplingFactors(const SamplingFactorType factor)
  {
    for (unsigned int i = 0; i < SpaceDimension; i++)
    {
      m_SamplingFactor[i] = factor[i];
    }
  }
 
  double
  GetEnergy0(double diff)
  {
    return (*m_InternalEnergyFunction)(diff);
  }
 
  MeasureType
  GetValue(const ParametersType & parameters) const override;
 
  void
  GetDerivative(const ParametersType & itkNotUsed(parameters), DerivativeType & itkNotUsed(derivative)) const override
  {}
 
  void
  InitializeDistributions(Array<double> classMeans, Array<double> classSigmas);
 
  unsigned int
  GetNumberOfParameters() const override;
 
protected:
  MRIBiasEnergyFunction();
 
  ~MRIBiasEnergyFunction() override;
 
private:
  BiasFieldType * m_BiasField;
 
  ImagePointer m_Image;
 
  MaskPointer m_Mask;
 
  ImageRegionType m_Region;
 
  InternalEnergyFunction * m_InternalEnergyFunction;
 
  SamplingFactorType m_SamplingFactor;
}; // end of class
 
template <typename TInputImage, typename TOutputImage, typename TMaskImage>
class ITK_TEMPLATE_EXPORT MRIBiasFieldCorrectionFilter : public ImageToImageFilter<TInputImage, TOutputImage>
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(MRIBiasFieldCorrectionFilter);
 
  using Self = MRIBiasFieldCorrectionFilter;
  using Superclass = ImageToImageFilter<TInputImage, TOutputImage>;
  using Pointer = SmartPointer<Self>;
  using ConstPointer = SmartPointer<const Self>;
 
  itkNewMacro(Self);
 
  itkTypeMacro(MRIBiasFieldCorrectionFilter, ImageToImageFilter);
 
  static constexpr unsigned int ImageDimension = TOutputImage::ImageDimension;
 
  using OutputImageType = TOutputImage;
  using OutputImagePointer = typename TOutputImage::Pointer;
  using OutputImageIndexType = typename TOutputImage::IndexType;
  using OutputImagePixelType = typename TOutputImage::PixelType;
  using OutputImageSizeType = typename TOutputImage::SizeType;
  using OutputImageRegionType = typename TOutputImage::RegionType;
 
  using InputImageType = TInputImage;
  using InputImagePointer = typename TInputImage::Pointer;
  using InputImageIndexType = typename TInputImage::IndexType;
  using InputImagePixelType = typename TInputImage::PixelType;
  using InputImageSizeType = typename TInputImage::SizeType;
  using InputImageRegionType = typename TInputImage::RegionType;
 
  using ImageMaskType = TMaskImage;
  using ImageMaskPointer = typename ImageMaskType::Pointer;
  using ImageMaskRegionType = typename ImageMaskType::RegionType;
 
  using InternalImageType = Image<float, Self::ImageDimension>;
  using InternalImagePixelType = typename InternalImageType::PixelType;
  using InternalImagePointer = typename InternalImageType::Pointer;
  using InternalImageRegionType = typename InternalImageType::RegionType;
 
  using MRASlabIdentifierType = MRASlabIdentifier<InputImageType>;
  using SlabRegionVectorType = typename MRASlabIdentifierType::SlabRegionVectorType;
  using SlabRegionVectorIteratorType = typename SlabRegionVectorType::iterator;
 
  using BiasFieldType = MultivariateLegendrePolynomial;
 
  using EnergyFunctionType = MRIBiasEnergyFunction<InternalImageType, ImageMaskType, BiasFieldType>;
  using EnergyFunctionPointer = typename EnergyFunctionType::Pointer;
 
  using NormalVariateGeneratorType = Statistics::NormalVariateGenerator;
 
  using OptimizerType = OnePlusOneEvolutionaryOptimizer;
 
  using ScheduleType = Array2D<unsigned int>;
 
  void
  SetInputMask(ImageMaskType * inputMask);
  itkGetModifiableObjectMacro(InputMask, ImageMaskType);
 
  void
  SetOutputMask(ImageMaskType * outputMask);
  itkGetModifiableObjectMacro(OutputMask, ImageMaskType);
 
#if defined(ITK_LEGACY_REMOVE)
 
  void
  IsBiasFieldMultiplicative(bool flag)
  {
    m_BiasFieldMultiplicative = flag;
  }
#endif // defined ( ITK_LEGACY_REMOVE )
 
  itkSetMacro(BiasFieldMultiplicative, bool);
  itkGetConstMacro(BiasFieldMultiplicative, bool);
  itkBooleanMacro(BiasFieldMultiplicative);
 
#if defined(ITK_LEGACY_REMOVE)
 
  bool
  IsBiasFieldMultiplicative()
  {
    return m_BiasFieldMultiplicative;
  }
#endif // defined ( ITK_LEGACY_REMOVE )
 
  itkSetMacro(UsingInterSliceIntensityCorrection, bool);
  itkGetConstMacro(UsingInterSliceIntensityCorrection, bool);
  itkBooleanMacro(UsingInterSliceIntensityCorrection);
 
  itkSetMacro(UsingSlabIdentification, bool);
  itkGetConstMacro(UsingSlabIdentification, bool);
  itkBooleanMacro(UsingSlabIdentification);
 
  itkSetMacro(SlabBackgroundMinimumThreshold, InputImagePixelType);
  itkGetConstReferenceMacro(SlabBackgroundMinimumThreshold, InputImagePixelType);
 
  itkSetMacro(SlabNumberOfSamples, unsigned int);
  itkGetConstReferenceMacro(SlabNumberOfSamples, unsigned int);
 
  itkSetMacro(SlabTolerance, double);
  itkGetConstReferenceMacro(SlabTolerance, double);
 
  itkSetMacro(UsingBiasFieldCorrection, bool);
  itkGetConstMacro(UsingBiasFieldCorrection, bool);
  itkBooleanMacro(UsingBiasFieldCorrection);
 
  itkSetMacro(GeneratingOutput, bool);
  itkGetConstMacro(GeneratingOutput, bool);
  itkBooleanMacro(GeneratingOutput);
 
  itkSetMacro(SlicingDirection, int);
  itkGetConstMacro(SlicingDirection, int);
 
  itkSetMacro(BiasFieldDegree, int);
  itkGetConstMacro(BiasFieldDegree, int);
 
  void
  SetInitialBiasFieldCoefficients(const BiasFieldType::CoefficientArrayType & coefficients)
  {
    this->Modified();
    m_BiasFieldCoefficients = coefficients;
  }
 
  BiasFieldType::CoefficientArrayType
  GetEstimatedBiasFieldCoefficients()
  {
    return m_EstimatedBiasFieldCoefficients;
  }
 
  void
  SetTissueClassStatistics(const Array<double> & means, const Array<double> & sigmas);
 
  itkSetMacro(VolumeCorrectionMaximumIteration, int);
  itkGetConstMacro(VolumeCorrectionMaximumIteration, int);
 
  itkSetMacro(InterSliceCorrectionMaximumIteration, int);
  itkGetConstMacro(InterSliceCorrectionMaximumIteration, int);
 
  itkSetMacro(OptimizerInitialRadius, double);
  itkGetConstMacro(OptimizerInitialRadius, double);
 
  itkSetMacro(OptimizerGrowthFactor, double);
  itkGetConstMacro(OptimizerGrowthFactor, double);
 
  itkSetMacro(OptimizerShrinkFactor, double);
  itkGetConstMacro(OptimizerShrinkFactor, double);
 
  void
  SetNumberOfLevels(unsigned int num);
 
  itkGetConstMacro(NumberOfLevels, unsigned int);
 
  void
  SetSchedule(const ScheduleType & schedule);
 
  itkGetConstReferenceMacro(Schedule, ScheduleType);
 
  void
  SetStartingShrinkFactors(unsigned int factor);
 
  void
  SetStartingShrinkFactors(const unsigned int * factors);
 
  const unsigned int *
  GetStartingShrinkFactors() const;
 
  static bool
  IsScheduleDownwardDivisible(const ScheduleType & schedule);
 
  void
  Initialize();
 
  BiasFieldType
  EstimateBiasField(InputImageRegionType region, unsigned int degree, int maximumIteration);
 
  void
  CorrectImage(BiasFieldType & bias, InputImageRegionType region);
 
  void
  CorrectInterSliceIntensityInhomogeneity(InputImageRegionType region);
 
protected:
  MRIBiasFieldCorrectionFilter();
  ~MRIBiasFieldCorrectionFilter() override = default;
  void
  PrintSelf(std::ostream & os, Indent indent) const override;
 
  bool
  CheckMaskImage(ImageMaskType * mask);
 
protected:
  void
  Log1PImage(InternalImageType * source, InternalImageType * target);
 
  void
  ExpImage(InternalImageType * source, InternalImageType * target);
 
  template <typename TSource, typename TTarget>
  void
  CopyAndConvertImage(const TSource * source, TTarget * target, typename TTarget::RegionType requestedRegion)
  {
    using SourceIterator = ImageRegionConstIterator<TSource>;
    using TargetIterator = ImageRegionIterator<TTarget>;
    using TargetPixelType = typename TTarget::PixelType;
 
    SourceIterator s_iter(source, requestedRegion);
    TargetIterator t_iter(target, requestedRegion);
 
    s_iter.GoToBegin();
    t_iter.GoToBegin();
    while (!s_iter.IsAtEnd())
    {
      t_iter.Set(static_cast<TargetPixelType>(s_iter.Get()));
      ++s_iter;
      ++t_iter;
    }
  }
 
  void
  GetBiasFieldSize(InputImageRegionType region, BiasFieldType::DomainSizeType & biasSize);
 
  void
  AdjustSlabRegions(SlabRegionVectorType & slabs, OutputImageRegionType requestedRegion);
 
  void
  GenerateData() override;
 
private:
  EnergyFunctionPointer m_EnergyFunction;
 
  NormalVariateGeneratorType::Pointer m_NormalVariateGenerator;
 
  ImageMaskPointer m_InputMask;
 
  ImageMaskPointer m_OutputMask;
 
  InternalImagePointer m_InternalInput;
 
  SlabRegionVectorType m_Slabs;
 
  int m_SlicingDirection;
 
  bool m_BiasFieldMultiplicative{ true };
 
  bool m_UsingInterSliceIntensityCorrection;
  bool m_UsingSlabIdentification{ false };
  bool m_UsingBiasFieldCorrection{ true };
  bool m_GeneratingOutput{ true };
 
  unsigned int        m_SlabNumberOfSamples{ 200 };
  InputImagePixelType m_SlabBackgroundMinimumThreshold;
  double              m_SlabTolerance{ 0.0 };
 
  int m_BiasFieldDegree{ 3 };
 
  unsigned int m_NumberOfLevels{ 0 };
 
  ScheduleType m_Schedule;
 
  BiasFieldType::CoefficientArrayType m_BiasFieldCoefficients;
 
  BiasFieldType::CoefficientArrayType m_EstimatedBiasFieldCoefficients;
 
  int m_VolumeCorrectionMaximumIteration{ 2000 };
 
  int m_InterSliceCorrectionMaximumIteration{ 4000 };
 
  double m_OptimizerInitialRadius{ 1.01 };
 
  double m_OptimizerGrowthFactor{ 1.05 };
 
  double m_OptimizerShrinkFactor;
 
  Array<double> m_TissueClassMeans;
 
  Array<double> m_TissueClassSigmas;
};
 
} // end namespace itk
 
#ifndef ITK_MANUAL_INSTANTIATION
#  include "itkMRIBiasFieldCorrectionFilter.hxx"
#endif
 
#endif