itkN4BiasFieldCorrectionImageFilter.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright Insight Software 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
 *
 *         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 itkN4BiasFieldCorrectionImageFilter_h
#define itkN4BiasFieldCorrectionImageFilter_h

#include "itkImageToImageFilter.h"

#include "itkArray.h"
#include "itkBSplineScatteredDataPointSetToImageFilter.h"
#include "itkPointSet.h"
#include "itkVector.h"

#include "vnl/vnl_vector.h"

namespace itk {

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

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

  itkTypeMacro( N4BiasFieldCorrectionImageFilter, ImageToImageFilter );

  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 )
    {
    ArrayType nlevels;

    nlevels.Fill( 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 );

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.

  RealImagePointer SharpenImage( const RealImageType * ) const;

  RealImagePointer UpdateBiasFieldEstimate( RealImageType *, std::size_t );

  RealImagePointer ReconstructBiasField( BiasFieldControlPointLatticeType * );

  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