itkMRIBiasFieldCorrectionFilter.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 __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< class TImage, class TImageMask, class TBiasField >
class MRIBiasEnergyFunction:public SingleValuedCostFunction
{
public:
  typedef MRIBiasEnergyFunction      Self;
  typedef SingleValuedCostFunction   Superclass;
  typedef SmartPointer< Self >       Pointer;
  typedef SmartPointer< const Self > ConstPointer;

  itkTypeMacro(MRIBiasEnergyFunction, SingleValuedCostFunction);

  itkNewMacro(Self);

  typedef TImage                         ImageType;
  typedef typename ImageType::Pointer    ImagePointer;
  typedef typename ImageType::PixelType  ImageElementType;
  typedef typename ImageType::IndexType  ImageIndexType;
  typedef typename ImageType::RegionType ImageRegionType;
  typedef TImageMask                     MaskType;
  typedef typename MaskType::Pointer     MaskPointer;
  typedef typename MaskType::PixelType   MaskElementType;

  typedef TBiasField BiasFieldType;

  typedef typename Superclass::ParametersType ParametersType;

  typedef Superclass::DerivativeType DerivativeType;

  typedef Superclass::MeasureType MeasureType;

  itkStaticConstMacro(SpaceDimension, unsigned int, 3);

  typedef CompositeValleyFunction InternalEnergyFunction;

  typedef unsigned int SamplingFactorType[SpaceDimension];

  itkSetObjectMacro(Image, ImageType);

  itkSetObjectMacro(Mask, MaskType);

  itkSetMacro(Region, ImageRegionType);

  void SetBiasField(BiasFieldType *bias)
  { m_BiasField = bias; }

  void SetSamplingFactors(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;

  void GetDerivative( const ParametersType & itkNotUsed(parameters),
                      DerivativeType & itkNotUsed(derivative) ) const
  {}

  void InitializeDistributions(Array< double > classMeans,
                               Array< double > classSigmas);

  unsigned int GetNumberOfParameters(void) const;

protected:
  MRIBiasEnergyFunction();

  virtual ~MRIBiasEnergyFunction();

private:

  BiasFieldType        *m_BiasField;

  ImagePointer m_Image;

  MaskPointer m_Mask;

  ImageRegionType m_Region;

  InternalEnergyFunction *m_InternalEnergyFunction;

  SamplingFactorType m_SamplingFactor;

  MRIBiasEnergyFunction(const Self &); //purposely not implemented
  void operator=(const Self &);        //purposely not implemented
};                                     // end of class

template< class TInputImage, class TOutputImage, class TMaskImage >
class ITK_EXPORT MRIBiasFieldCorrectionFilter:
  public ImageToImageFilter< TInputImage, TOutputImage >
{
public:
  typedef MRIBiasFieldCorrectionFilter                    Self;
  typedef ImageToImageFilter< TInputImage, TOutputImage > Superclass;
  typedef SmartPointer< Self >                            Pointer;
  typedef SmartPointer< const Self >                      ConstPointer;

  itkNewMacro(Self);

  itkTypeMacro(MRIBiasFieldCorrectionFilter, ImageToImageFilter);

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

  typedef TOutputImage                      OutputImageType;
  typedef typename TOutputImage::Pointer    OutputImagePointer;
  typedef typename TOutputImage::IndexType  OutputImageIndexType;
  typedef typename TOutputImage::PixelType  OutputImagePixelType;
  typedef typename TOutputImage::SizeType   OutputImageSizeType;
  typedef typename TOutputImage::RegionType OutputImageRegionType;

  typedef TInputImage                      InputImageType;
  typedef typename TInputImage::Pointer    InputImagePointer;
  typedef typename TInputImage::IndexType  InputImageIndexType;
  typedef typename TInputImage::PixelType  InputImagePixelType;
  typedef typename TInputImage::SizeType   InputImageSizeType;
  typedef typename TInputImage::RegionType InputImageRegionType;

  typedef TMaskImage                         ImageMaskType;
  typedef typename ImageMaskType::Pointer    ImageMaskPointer;
  typedef typename ImageMaskType::RegionType ImageMaskRegionType;

  typedef Image< float, itkGetStaticConstMacro(ImageDimension) > InternalImageType;
  typedef typename InternalImageType::PixelType                  InternalImagePixelType;
  typedef typename InternalImageType::Pointer                    InternalImagePointer;
  typedef typename InternalImageType::RegionType                 InternalImageRegionType;

  typedef MRASlabIdentifier< InputImageType >                  MRASlabIdentifierType;
  typedef typename MRASlabIdentifierType::SlabRegionVectorType SlabRegionVectorType;
  typedef typename SlabRegionVectorType::iterator              SlabRegionVectorIteratorType;

  typedef MultivariateLegendrePolynomial BiasFieldType;

  typedef MRIBiasEnergyFunction< InternalImageType,
                                 ImageMaskType,
                                 BiasFieldType >            EnergyFunctionType;
  typedef typename EnergyFunctionType::Pointer EnergyFunctionPointer;

  typedef Statistics::NormalVariateGenerator NormalVariateGeneratorType;

  typedef OnePlusOneEvolutionaryOptimizer OptimizerType;

  typedef Array2D< unsigned int > ScheduleType;

  void SetInputMask(ImageMaskType *inputMask);
  itkGetModifiableObjectMacro(InputMask, ImageMaskType);

  void SetOutputMask(ImageMaskType *outputMask);
  itkGetModifiableObjectMacro(OutputMask, ImageMaskType);

  void IsBiasFieldMultiplicative(bool flag)
  { m_BiasMultiplicative = flag; }

  bool IsBiasFieldMultiplicative()
  { return m_BiasMultiplicative; }

  itkSetMacro(UsingInterSliceIntensityCorrection, bool);
  itkGetConstReferenceMacro(UsingInterSliceIntensityCorrection, bool);

  itkSetMacro(UsingSlabIdentification, bool);
  itkGetConstReferenceMacro(UsingSlabIdentification, bool);

  itkSetMacro(SlabBackgroundMinimumThreshold, InputImagePixelType);
  itkGetConstReferenceMacro(SlabBackgroundMinimumThreshold,
                            InputImagePixelType);

  itkSetMacro(SlabNumberOfSamples, unsigned int);
  itkGetConstReferenceMacro(SlabNumberOfSamples, unsigned int);

  itkSetMacro(SlabTolerance, double);
  itkGetConstReferenceMacro(SlabTolerance, double);

  itkSetMacro(UsingBiasFieldCorrection, bool);
  itkGetConstReferenceMacro(UsingBiasFieldCorrection, bool);

  itkSetMacro(GeneratingOutput, bool);
  itkGetConstReferenceMacro(GeneratingOutput, bool);

  itkSetMacro(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)
  throw ( ExceptionObject );

  itkSetMacro(VolumeCorrectionMaximumIteration, int);
  itkGetConstMacro(VolumeCorrectionMaximumIteration, int);
  itkSetMacro(InterSliceCorrectionMaximumIteration, int);
  itkGetConstMacro(InterSliceCorrectionMaximumIteration, int);

  void SetOptimizerInitialRadius(double initRadius)
  { m_OptimizerInitialRadius = initRadius; }
  double GetOptimizerInitialRadius()
  { return m_OptimizerInitialRadius; }

  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(unsigned int *factors);

  const unsigned int * GetStartingShrinkFactors() const;

  static bool IsScheduleDownwardDivisible(const ScheduleType & schedule);

  void Initialize()
  throw ( ExceptionObject );

  BiasFieldType EstimateBiasField(InputImageRegionType region,
                                  unsigned int degree,
                                  int maximumIteration);

  void CorrectImage(BiasFieldType & bias,
                    InputImageRegionType region);

  void CorrectInterSliceIntensityInhomogeneity(InputImageRegionType region);

protected:
  MRIBiasFieldCorrectionFilter();
  virtual ~MRIBiasFieldCorrectionFilter();
  void PrintSelf(std::ostream & os, Indent indent) const;

  bool CheckMaskImage(ImageMaskType *mask);

protected:
  void Log1PImage(InternalImageType *source,
                  InternalImageType *target);

  void ExpImage(InternalImageType *source,
                InternalImageType *target);

  template< class TSource, class TTarget >
  void CopyAndConvertImage(const TSource *source,
                           TTarget *target,
                           typename TTarget::RegionType requestedRegion)
  {
    typedef ImageRegionConstIterator< TSource > SourceIterator;
    typedef ImageRegionIterator< TTarget >      TargetIterator;
    typedef typename TTarget::PixelType         TargetPixelType;

    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 & domainSize);

  void AdjustSlabRegions(SlabRegionVectorType & slabs,
                         OutputImageRegionType requestedRegion);

  void GenerateData();

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

  EnergyFunctionPointer m_EnergyFunction;

  NormalVariateGeneratorType::Pointer m_NormalVariateGenerator;

  ImageMaskPointer m_InputMask;

  ImageMaskPointer m_OutputMask;

  InternalImagePointer m_InternalInput;

  SlabRegionVectorType m_Slabs;

  int m_SlicingDirection;

  bool m_BiasMultiplicative;

  bool m_UsingInterSliceIntensityCorrection;
  bool m_UsingSlabIdentification;
  bool m_UsingBiasFieldCorrection;
  bool m_GeneratingOutput;

  unsigned int        m_SlabNumberOfSamples;
  InputImagePixelType m_SlabBackgroundMinimumThreshold;
  double              m_SlabTolerance;
  int m_BiasFieldDegree;

  unsigned int m_NumberOfLevels;

  ScheduleType m_Schedule;

  BiasFieldType::CoefficientArrayType m_BiasFieldCoefficients;

  BiasFieldType::CoefficientArrayType m_EstimatedBiasFieldCoefficients;

  int m_VolumeCorrectionMaximumIteration;

  int m_InterSliceCorrectionMaximumIteration;

  double m_OptimizerInitialRadius;

  double m_OptimizerGrowthFactor;

  double m_OptimizerShrinkFactor;

  Array< double > m_TissueClassMeans;

  Array< double > m_TissueClassSigmas;
};

// ==================================
} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkMRIBiasFieldCorrectionFilter.hxx"
#endif

#endif