itkMRIBiasFieldCorrectionFilter.h

Go to the documentation of this file.

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkMRIBiasFieldCorrectionFilter.h,v $
  Language:  C++
  Date:      $Date: 2006/03/19 04:36:54 $
  Version:   $Revision: 1.36 $

  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( SingleValuedCostFunction, CostFunction );

  itkNewMacro(Self);

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

  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;

private:
  
  BiasFieldType        * m_BiasField ;

  ImagePointer           m_Image ;

  MaskPointer            m_Mask ;

  ImageRegionType        m_Region ;

  InternalEnergyFunction* m_InternalEnergyFunction ;

  SamplingFactorType m_SamplingFactor;

protected:
  MRIBiasEnergyFunction();

  virtual ~MRIBiasEnergyFunction();


private:
  
  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 TInputImage InputImageType ;
  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 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

---