itkMRIBiasFieldCorrectionFilter.h

Go to the documentation of this file.

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkMRIBiasFieldCorrectionFilter.h,v $
  Language:  C++
  Date:      $Date: 2002/12/24 03:08:44 $
  Version:   $Revision: 1.24 $

  Copyright (c) 2002 Insight 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 "itkImageToImageFilter.h"
#include "itkImage.h"

#include "itkMRASlabIdentifier.h"
#include "itkCompositeValleyFunction.h"
#include "itkMultivariateLegendrePolynomial.h"
#include "Statistics/itkNormalVariateGenerator.h"
#include "itkOnePlusOneEvolutionaryOptimizer.h"
#include "itkArray.h"
#include "itkImageRegionConstIterator.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;
  
  itkNewMacro(Self);
  
  itkTypeMacro( SingleValuedCostFunction, CostFunction );


  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 ;

  itkSetObjectMacro( Image, ImageType );

  itkSetObjectMacro( Mask, MaskType );

  itkSetMacro( Region, ImageRegionType );

  itkSetObjectMacro( BiasField, BiasFieldType );

  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 ;

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 ;

  void SetInputMask(ImageMaskPointer inputMask);
  itkGetObjectMacro( InputMask, ImageMaskType );

  void SetOutputMask(ImageMaskPointer outputMask) ;

  itkGetObjectMacro( OutputMask, ImageMaskType );

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

  bool IsBiasFieldMultiplicative() 
    { return m_BiasMultiplicative ; }

  itkSetMacro( UsingInterSliceIntensityCorrection, bool );

   itkSetMacro( UsingSlabIdentification, bool );

  itkSetMacro( UsingBiasFieldCorrection, bool );

  itkSetMacro( GeneratingOutput, bool );

  itkSetMacro( SlicingDirection , int );

  itkSetMacro( BiasFieldDegree, int );
  itkGetMacro( BiasFieldDegree, int );

  itkGetMacro( NoOfBiasFieldCoefficients, int );

  itkGetMacro( BiasFieldDomainSize, BiasFieldType::DomainSizeType );

  void SetInitialBiasFieldCoefficients
  (const BiasFieldType::CoefficientArrayType &coefficients)
    { m_BiasFieldCoefficients = coefficients ; }

  itkGetMacro( EstimatedBiasFieldCoefficients, BiasFieldType::CoefficientArrayType );

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

  void EstimateBiasField(BiasFieldType* bias,
                         InputImageRegionType region) ;

  void CorrectImage(BiasFieldType* bias, 
                    InputImageRegionType region) ;

  void CorrectInterSliceIntensityInhomogeneity(InputImageRegionType region) ;

protected:
  MRIBiasFieldCorrectionFilter() ;
  virtual ~MRIBiasFieldCorrectionFilter() ;

  bool CheckMaskImage(ImageMaskPointer mask) ;

protected:
  void Log1PImage(InternalImagePointer source, 
                InternalImagePointer target) ;

  void ExpImage(InternalImagePointer source, 
                InternalImagePointer target) ;

  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 ;

  OptimizerType::Pointer m_Optimizer ;

  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 ;

  int m_BiasFieldDegree ;

  int m_BiasFieldDimension ;

  int m_NoOfBiasFieldCoefficients ;

  BiasFieldType::DomainSizeType m_BiasFieldDomainSize ;

  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 ;
};

  
// ==================================

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

  typename TSource::ConstPointer source( sourceInp );
  typename TTarget::Pointer      target( targetInp );

  target->SetRequestedRegion(requestedRegion) ;
  target->SetBufferedRegion(requestedRegion) ;
  target->Allocate() ;

  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 ;
    }
}

} // end namespace itk

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

#endif

---