itkN4MRIBiasFieldCorrectionImageFilter.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 __itkN4MRIBiasFieldCorrectionImageFilter_h
#define __itkN4MRIBiasFieldCorrectionImageFilter_h

#include "itkImageToImageFilter.h"

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

#include "vnl/vnl_vector.h"

namespace itk {

template<class TInputImage, class TMaskImage =
  Image<unsigned char, ::itk::GetImageDimension<TInputImage>::ImageDimension>,
  class TOutputImage = TInputImage>
class ITK_EXPORT N4MRIBiasFieldCorrectionImageFilter :
  public ImageToImageFilter<TInputImage, TOutputImage>
{
public:
  typedef N4MRIBiasFieldCorrectionImageFilter           Self;
  typedef ImageToImageFilter<TInputImage, TOutputImage> Superclass;
  typedef SmartPointer<Self>                            Pointer;
  typedef SmartPointer<const Self>                      ConstPointer;

  itkTypeMacro( N4MRIBiasFieldCorrectionImageFilter, ImageToImageFilter );

  itkNewMacro( Self );

  itkStaticConstMacro( ImageDimension, unsigned int,
                       TInputImage::ImageDimension );

  typedef TInputImage                       InputImageType;
  typedef TOutputImage                      OutputImageType;
  typedef TMaskImage                        MaskImageType;
  typedef typename MaskImageType::PixelType MaskPixelType;

  typedef float                           RealType;
  typedef Image<RealType, ImageDimension> RealImageType;
  typedef typename RealImageType::Pointer RealImagePointer;
  typedef Array<unsigned int>             VariableSizeArrayType;

  typedef Vector<RealType, 1>                                             ScalarType;
  typedef PointSet<ScalarType, itkGetStaticConstMacro( ImageDimension )>  PointSetType;
  typedef Image<ScalarType, itkGetStaticConstMacro( ImageDimension )>     ScalarImageType;
  typedef typename PointSetType::Pointer                                  PointSetPointer;
  typedef typename PointSetType::PointType                                PointType;

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

  void SetMaskImage( const MaskImageType *mask )
    {
    this->SetNthInput( 1, const_cast<MaskImageType *>( mask ) );
    }

  const MaskImageType* GetMaskImage() const
    {
    return static_cast<const MaskImageType*>( this->ProcessObject::GetInput( 1 ) );
    }

  void SetConfidenceImage( const RealImageType *image )
    {
    this->SetNthInput( 2, const_cast<RealImageType *>( image ) );
    }

  const RealImageType* GetConfidenceImage() const
    {
    return static_cast<const RealImageType*>( this->ProcessObject::GetInput( 2 ) );
    }

  itkSetMacro( MaskLabel, MaskPixelType );

  itkGetConstMacro( MaskLabel, MaskPixelType );

  // 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 );

  itkGetConstMacro( LogBiasFieldControlPointLattice,
                    typename BiasFieldControlPointLatticeType::Pointer );

  itkGetConstMacro( ElapsedIterations, unsigned int );

  itkGetConstMacro( CurrentConvergenceMeasurement, RealType );

  itkGetConstMacro( CurrentLevel, unsigned int );

protected:
  N4MRIBiasFieldCorrectionImageFilter();
  ~N4MRIBiasFieldCorrectionImageFilter() {}
  void PrintSelf( std::ostream& os, Indent indent ) const;

  void GenerateData();

private:
  N4MRIBiasFieldCorrectionImageFilter( const Self& ); //purposely not
                                                      // implemented
  void operator=( const Self& );                      //purposely not
                                                      // implemented

  // 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 * );

  RealType CalculateConvergenceMeasurement( const RealImageType *, const RealImageType * ) const;


  MaskPixelType m_MaskLabel;

  // Parameters for deconvolution with Wiener filter

  unsigned int m_NumberOfHistogramBins;
  RealType     m_WienerFilterNoise;
  RealType     m_BiasFieldFullWidthAtHalfMaximum;

  // Convergence parameters

  VariableSizeArrayType m_MaximumNumberOfIterations;
  unsigned int          m_ElapsedIterations;
  RealType              m_ConvergenceThreshold;
  RealType              m_CurrentConvergenceMeasurement;
  unsigned int          m_CurrentLevel;

  // B-spline fitting parameters

  typename
  BiasFieldControlPointLatticeType::Pointer m_LogBiasFieldControlPointLattice;

  unsigned int m_SplineOrder;
  ArrayType    m_NumberOfControlPoints;
  ArrayType    m_NumberOfFittingLevels;

};

} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkN4MRIBiasFieldCorrectionImageFilter.txx"
#endif

#endif