itkN4BiasFieldCorrectionImageFilter.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 itkN4BiasFieldCorrectionImageFilter_h
#define itkN4BiasFieldCorrectionImageFilter_h
 
#include "itkImageToImageFilter.h"
 
#include "itkArray.h"
#include "itkBSplineControlPointImageFilter.h"
#include "itkBSplineScatteredDataPointSetToImageFilter.h"
#include "itkPointSet.h"
#include "itkVector.h"
 
#include <vector>
 
#include "vnl/vnl_vector.h"
 
namespace itk
{

 
template <typename TInputImage,
          typename TMaskImage = Image<unsigned char, TInputImage::ImageDimension>,
          class TOutputImage = TInputImage>
class ITK_TEMPLATE_EXPORT N4BiasFieldCorrectionImageFilter : public ImageToImageFilter<TInputImage, TOutputImage>
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(N4BiasFieldCorrectionImageFilter);

  using Self = N4BiasFieldCorrectionImageFilter;
  using Superclass = ImageToImageFilter<TInputImage, TOutputImage>;
  using Pointer = SmartPointer<Self>;
  using ConstPointer = SmartPointer<const Self>;

  itkOverrideGetNameOfClassMacro(N4BiasFieldCorrectionImageFilter);

  itkNewMacro(Self);

  static constexpr unsigned int ImageDimension = TInputImage::ImageDimension;

  using InputImageType = TInputImage;
  using OutputImageType = TOutputImage;
  using MaskImageType = TMaskImage;
  using MaskPixelType = typename MaskImageType::PixelType;
 
  using RealType = float;
  using RealImageType = Image<RealType, ImageDimension>;
  using RealImagePointer = typename RealImageType::Pointer;
  using VariableSizeArrayType = Array<unsigned int>;

  using ScalarType = Vector<RealType, 1>;
  using PointSetType = PointSet<ScalarType, Self::ImageDimension>;
  using ScalarImageType = Image<ScalarType, Self::ImageDimension>;
  using PointSetPointer = typename PointSetType::Pointer;
  using PointType = typename PointSetType::PointType;

  using BSplineFilterType = BSplineScatteredDataPointSetToImageFilter<PointSetType, ScalarImageType>;
  using BiasFieldControlPointLatticeType = typename BSplineFilterType::PointDataImageType;
  using ArrayType = typename BSplineFilterType::ArrayType;

  void
  EnlargeOutputRequestedRegion(DataObject *) override;

  void
  SetInput1(const InputImageType * image)
  {
    this->SetInput(image);
  }

  itkSetInputMacro(MaskImage, MaskImageType);
  void
  SetInput2(const MaskImageType * mask)
  {
    this->SetMaskImage(mask);
  }


  itkGetInputMacro(MaskImage, MaskImageType);

  itkSetMacro(MaskLabel, MaskPixelType);
  itkGetConstMacro(MaskLabel, MaskPixelType);

  itkSetMacro(UseMaskLabel, bool);
  itkGetConstMacro(UseMaskLabel, bool);
  itkBooleanMacro(UseMaskLabel);

  itkSetInputMacro(ConfidenceImage, RealImageType);
  void
  SetInput3(const RealImageType * image)
  {
    this->SetConfidenceImage(image);
  }


  itkGetInputMacro(ConfidenceImage, RealImageType);
 
  // Sharpen histogram parameters: in estimating the bias field, the
  // first step is to sharpen the intensity histogram by Wiener deconvolution
  // with a 1-D Gaussian.  The following parameters define this operation.
  // These default values in N4 match the default values in N3.

  itkSetMacro(NumberOfHistogramBins, unsigned int);

  itkGetConstMacro(NumberOfHistogramBins, unsigned int);

  itkSetMacro(WienerFilterNoise, RealType);

  itkGetConstMacro(WienerFilterNoise, RealType);

  itkSetMacro(BiasFieldFullWidthAtHalfMaximum, RealType);

  itkGetConstMacro(BiasFieldFullWidthAtHalfMaximum, RealType);
 
  // B-spline parameters governing the fitting routine

  itkSetMacro(SplineOrder, unsigned int);

  itkGetConstMacro(SplineOrder, unsigned int);

  itkSetMacro(NumberOfControlPoints, ArrayType);

  itkGetConstMacro(NumberOfControlPoints, ArrayType);

  itkSetMacro(NumberOfFittingLevels, ArrayType);

  void
  SetNumberOfFittingLevels(unsigned int n)
  {
    auto nlevels = MakeFilled<ArrayType>(n);
    this->SetNumberOfFittingLevels(nlevels);
  }

  itkGetConstMacro(NumberOfFittingLevels, ArrayType);

  itkSetMacro(MaximumNumberOfIterations, VariableSizeArrayType);

  itkGetConstMacro(MaximumNumberOfIterations, VariableSizeArrayType);

  itkSetMacro(ConvergenceThreshold, RealType);

  itkGetConstMacro(ConvergenceThreshold, RealType);

  itkGetConstObjectMacro(LogBiasFieldControlPointLattice, BiasFieldControlPointLatticeType);

  itkGetConstMacro(ElapsedIterations, unsigned int);

  itkGetConstMacro(CurrentConvergenceMeasurement, RealType);

  itkGetConstMacro(CurrentLevel, unsigned int);

  RealImagePointer
  ReconstructBiasField(const BiasFieldControlPointLatticeType *);
 
protected:
  N4BiasFieldCorrectionImageFilter();
  ~N4BiasFieldCorrectionImageFilter() override = default;
  void
  PrintSelf(std::ostream & os, Indent indent) const override;
 
  void
  GenerateData() override;
 
private:
  // N4 algorithm functions:  The basic algorithm iterates between sharpening
  // the intensity histogram of the corrected input image and spatially
  // smoothing those results with a B-spline scalar field estimate of the
  // bias field.  The former is handled by the function SharpenImage()
  // whereas the latter is handled by the function UpdateBiasFieldEstimate().
  // Convergence is determined by the coefficient of variation of the difference
  // image between the current bias field estimate and the previous estimate.

  void
  SharpenImage(const RealImageType * unsharpenedImage, RealImageType * sharpenedImage) const;

  struct FitPoints
  {
    typename PointSetType::PointsContainerPointer             points;
    typename PointSetType::PointDataContainerPointer          values;
    typename BSplineFilterType::WeightsContainerType::Pointer weights;
    std::vector<SizeValueType>                                bufferIndices;
  };

  void
  UpdateBiasFieldEstimate(RealImageType *, FitPoints &, RealImageType *);

  void
  ReconstructBiasFieldInto(const BiasFieldControlPointLatticeType *, RealImageType *);
 
  using BiasFieldReconstructerType = BSplineControlPointImageFilter<BiasFieldControlPointLatticeType, ScalarImageType>;

  typename BiasFieldReconstructerType::Pointer
  MakeBiasFieldReconstructer(const BiasFieldControlPointLatticeType *) const;

  RealType
  CalculateConvergenceMeasurement(const RealImageType *, const RealImageType *) const;
 
  MaskPixelType m_MaskLabel{};
  bool          m_UseMaskLabel{ false };
 
  // Parameters for deconvolution with Wiener filter
 
  unsigned int m_NumberOfHistogramBins{ 200 };
  RealType     m_WienerFilterNoise{ static_cast<RealType>(0.01) };
  RealType     m_BiasFieldFullWidthAtHalfMaximum{ static_cast<RealType>(0.15) };
 
  // Convergence parameters
 
  VariableSizeArrayType m_MaximumNumberOfIterations{};
  unsigned int          m_ElapsedIterations{ 0 };
  RealType              m_ConvergenceThreshold{ static_cast<RealType>(0.001) };
  RealType              m_CurrentConvergenceMeasurement{};
  unsigned int          m_CurrentLevel{ 0 };
 
  // B-spline fitting parameters
 
  typename BiasFieldControlPointLatticeType::Pointer m_LogBiasFieldControlPointLattice{};
 
  unsigned int m_SplineOrder{ 3 };
  ArrayType    m_NumberOfControlPoints{};
  ArrayType    m_NumberOfFittingLevels{};
};
 
} // end namespace itk
 
#ifndef ITK_MANUAL_INSTANTIATION
#  include "itkN4BiasFieldCorrectionImageFilter.hxx"
#endif
 
#endif