itkMRIBiasFieldCorrectionFilter.h

Go to the documentation of this file.

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkMRIBiasFieldCorrectionFilter.h,v $
  Language:  C++
  Date:      $Date: 2009-01-26 21:45:51 $
  Version:   $Revision: 1.38 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#ifndef __itkMRIBiasFieldCorrectionFilter_h
#define __itkMRIBiasFieldCorrectionFilter_h

#include <time.h>

#include "itkImageToImageFilter.h"
#include "itkImage.h"
#include "itkArray2D.h"
#include "itkMRASlabIdentifier.h"
#include "itkCompositeValleyFunction.h"
#include "itkMultivariateLegendrePolynomial.h"
#include "Statistics/itkNormalVariateGenerator.h"
#include "itkOnePlusOneEvolutionaryOptimizer.h"
#include "itkArray.h"
#include "itkImageRegionConstIterator.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);
  itkGetObjectMacro( InputMask, ImageMaskType );

  void SetOutputMask(ImageMaskType* outputMask);

  itkGetObjectMacro( 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 );
  itkGetMacro( 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 );
  itkGetMacro( VolumeCorrectionMaximumIteration, int );
  itkSetMacro( InterSliceCorrectionMaximumIteration, int );
  itkGetMacro( InterSliceCorrectionMaximumIteration, int );

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

  itkSetMacro( OptimizerGrowthFactor, double );
  itkGetMacro( OptimizerGrowthFactor, double );

  itkSetMacro( OptimizerShrinkFactor, double );
  itkGetMacro( OptimizerShrinkFactor, double );


  void SetNumberOfLevels(unsigned int num);

  itkGetMacro(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.txx"
#endif

#endif