ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkMRIBiasFieldCorrectionFilter.h
Go to the documentation of this file.
00001 /*=========================================================================
00002  *
00003  *  Copyright Insight Software Consortium
00004  *
00005  *  Licensed under the Apache License, Version 2.0 (the "License");
00006  *  you may not use this file except in compliance with the License.
00007  *  You may obtain a copy of the License at
00008  *
00009  *         http://www.apache.org/licenses/LICENSE-2.0.txt
00010  *
00011  *  Unless required by applicable law or agreed to in writing, software
00012  *  distributed under the License is distributed on an "AS IS" BASIS,
00013  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00014  *  See the License for the specific language governing permissions and
00015  *  limitations under the License.
00016  *
00017  *=========================================================================*/
00018 #ifndef __itkMRIBiasFieldCorrectionFilter_h
00019 #define __itkMRIBiasFieldCorrectionFilter_h
00020 
00021 #include <time.h>
00022 
00023 #include "itkImageToImageFilter.h"
00024 #include "itkArray2D.h"
00025 #include "itkMRASlabIdentifier.h"
00026 #include "itkCompositeValleyFunction.h"
00027 #include "itkMultivariateLegendrePolynomial.h"
00028 #include "itkNormalVariateGenerator.h"
00029 #include "itkOnePlusOneEvolutionaryOptimizer.h"
00030 #include "itkImageRegionIterator.h"
00031 
00032 namespace itk
00033 {
00046 template< class TImage, class TImageMask, class TBiasField >
00047 class MRIBiasEnergyFunction:public SingleValuedCostFunction
00048 {
00049 public:
00051   typedef MRIBiasEnergyFunction      Self;
00052   typedef SingleValuedCostFunction   Superclass;
00053   typedef SmartPointer< Self >       Pointer;
00054   typedef SmartPointer< const Self > ConstPointer;
00055 
00057   itkTypeMacro(MRIBiasEnergyFunction, SingleValuedCostFunction);
00058 
00060   itkNewMacro(Self);
00061 
00063   typedef TImage                         ImageType;
00064   typedef typename ImageType::Pointer    ImagePointer;
00065   typedef typename ImageType::PixelType  ImageElementType;
00066   typedef typename ImageType::IndexType  ImageIndexType;
00067   typedef typename ImageType::RegionType ImageRegionType;
00068   typedef TImageMask                     MaskType;
00069   typedef typename MaskType::Pointer     MaskPointer;
00070   typedef typename MaskType::PixelType   MaskElementType;
00071 
00073   typedef TBiasField BiasFieldType;
00074 
00077   typedef typename Superclass::ParametersType ParametersType;
00078 
00080   typedef Superclass::DerivativeType DerivativeType;
00081 
00083   typedef Superclass::MeasureType MeasureType;
00084 
00085   itkStaticConstMacro(SpaceDimension, unsigned int, 3);
00086 
00088   typedef CompositeValleyFunction InternalEnergyFunction;
00089 
00091   typedef unsigned int SamplingFactorType[SpaceDimension];
00092 
00094   itkSetObjectMacro(Image, ImageType);
00095 
00097   itkSetObjectMacro(Mask, MaskType);
00098 
00100   itkSetMacro(Region, ImageRegionType);
00101 
00103   void SetBiasField(BiasFieldType *bias)
00104   { m_BiasField = bias; }
00105 
00108   void SetSamplingFactors(SamplingFactorType factor)
00109   {
00110     for ( unsigned int i = 0; i < SpaceDimension; i++ )
00111       {
00112       m_SamplingFactor[i] = factor[i];
00113       }
00114   }
00116 
00119   double GetEnergy0(double diff)
00120   {
00121     return ( *m_InternalEnergyFunction )( diff );
00122   }
00123 
00126   MeasureType GetValue(const ParametersType & parameters) const;
00127 
00130   void GetDerivative( const ParametersType & itkNotUsed(parameters),
00131                       DerivativeType & itkNotUsed(derivative) ) const
00132   {}
00133 
00138   void InitializeDistributions(Array< double > classMeans,
00139                                Array< double > classSigmas);
00140 
00141   unsigned int GetNumberOfParameters(void) const;
00142 
00143 protected:
00145   MRIBiasEnergyFunction();
00146 
00148   virtual ~MRIBiasEnergyFunction();
00149 private:
00150 
00152   BiasFieldType        *m_BiasField;
00153 
00155   ImagePointer m_Image;
00156 
00158   MaskPointer m_Mask;
00159 
00161   ImageRegionType m_Region;
00162 
00164   InternalEnergyFunction *m_InternalEnergyFunction;
00165 
00167   SamplingFactorType m_SamplingFactor;
00168 
00169   MRIBiasEnergyFunction(const Self &); //purposely not implemented
00170   void operator=(const Self &);        //purposely not implemented
00171 };                                     // end of class
00172 
00223 template< class TInputImage, class TOutputImage, class TMaskImage >
00224 class ITK_EXPORT MRIBiasFieldCorrectionFilter:
00225   public ImageToImageFilter< TInputImage, TOutputImage >
00226 {
00227 public:
00229   typedef MRIBiasFieldCorrectionFilter                    Self;
00230   typedef ImageToImageFilter< TInputImage, TOutputImage > Superclass;
00231   typedef SmartPointer< Self >                            Pointer;
00232   typedef SmartPointer< const Self >                      ConstPointer;
00233 
00235   itkNewMacro(Self);
00236 
00238   itkTypeMacro(MRIBiasFieldCorrectionFilter, ImageToImageFilter);
00239 
00241   itkStaticConstMacro(ImageDimension, unsigned int,
00242                       TOutputImage::ImageDimension);
00243 
00245   typedef TOutputImage                      OutputImageType;
00246   typedef typename TOutputImage::Pointer    OutputImagePointer;
00247   typedef typename TOutputImage::IndexType  OutputImageIndexType;
00248   typedef typename TOutputImage::PixelType  OutputImagePixelType;
00249   typedef typename TOutputImage::SizeType   OutputImageSizeType;
00250   typedef typename TOutputImage::RegionType OutputImageRegionType;
00251 
00252   typedef TInputImage                      InputImageType;
00253   typedef typename TInputImage::Pointer    InputImagePointer;
00254   typedef typename TInputImage::IndexType  InputImageIndexType;
00255   typedef typename TInputImage::PixelType  InputImagePixelType;
00256   typedef typename TInputImage::SizeType   InputImageSizeType;
00257   typedef typename TInputImage::RegionType InputImageRegionType;
00258 
00260   typedef TMaskImage                         ImageMaskType;
00261   typedef typename ImageMaskType::Pointer    ImageMaskPointer;
00262   typedef typename ImageMaskType::RegionType ImageMaskRegionType;
00263 
00265   typedef Image< float, itkGetStaticConstMacro(ImageDimension) > InternalImageType;
00266   typedef typename InternalImageType::PixelType                  InternalImagePixelType;
00267   typedef typename InternalImageType::Pointer                    InternalImagePointer;
00268   typedef typename InternalImageType::RegionType                 InternalImageRegionType;
00270 
00272   typedef MRASlabIdentifier< InputImageType >                  MRASlabIdentifierType;
00273   typedef typename MRASlabIdentifierType::SlabRegionVectorType SlabRegionVectorType;
00274   typedef typename SlabRegionVectorType::iterator              SlabRegionVectorIteratorType;
00275 
00277   typedef MultivariateLegendrePolynomial BiasFieldType;
00278 
00280   typedef MRIBiasEnergyFunction< InternalImageType,
00281                                  ImageMaskType,
00282                                  BiasFieldType >            EnergyFunctionType;
00283   typedef typename EnergyFunctionType::Pointer EnergyFunctionPointer;
00284 
00286   typedef Statistics::NormalVariateGenerator NormalVariateGeneratorType;
00287 
00289   typedef OnePlusOneEvolutionaryOptimizer OptimizerType;
00290 
00292   typedef Array2D< unsigned int > ScheduleType;
00293 
00297   void SetInputMask(ImageMaskType *inputMask);
00298 
00299   itkGetObjectMacro(InputMask, ImageMaskType);
00300 
00303   void SetOutputMask(ImageMaskType *outputMask);
00304 
00306   itkGetObjectMacro(OutputMask, ImageMaskType);
00307 
00311   void IsBiasFieldMultiplicative(bool flag)
00312   { m_BiasMultiplicative = flag; }
00313 
00315   bool IsBiasFieldMultiplicative()
00316   { return m_BiasMultiplicative; }
00317 
00321   itkSetMacro(UsingInterSliceIntensityCorrection, bool);
00322   itkGetConstReferenceMacro(UsingInterSliceIntensityCorrection, bool);
00324 
00330   itkSetMacro(UsingSlabIdentification, bool);
00331   itkGetConstReferenceMacro(UsingSlabIdentification, bool);
00333 
00334   itkSetMacro(SlabBackgroundMinimumThreshold, InputImagePixelType);
00335   itkGetConstReferenceMacro(SlabBackgroundMinimumThreshold,
00336                             InputImagePixelType);
00337 
00338   itkSetMacro(SlabNumberOfSamples, unsigned int);
00339   itkGetConstReferenceMacro(SlabNumberOfSamples, unsigned int);
00340 
00341   itkSetMacro(SlabTolerance, double);
00342   itkGetConstReferenceMacro(SlabTolerance, double);
00343 
00349   itkSetMacro(UsingBiasFieldCorrection, bool);
00350   itkGetConstReferenceMacro(UsingBiasFieldCorrection, bool);
00352 
00355   itkSetMacro(GeneratingOutput, bool);
00356   itkGetConstReferenceMacro(GeneratingOutput, bool);
00358 
00361   itkSetMacro(SlicingDirection, int);
00362 
00364   itkSetMacro(BiasFieldDegree, int);
00365   itkGetConstMacro(BiasFieldDegree, int);
00367 
00370   void SetInitialBiasFieldCoefficients(const
00371                                        BiasFieldType::CoefficientArrayType
00372                                        & coefficients)
00373   { this->Modified(); m_BiasFieldCoefficients = coefficients; }
00374 
00378   BiasFieldType::CoefficientArrayType GetEstimatedBiasFieldCoefficients()
00379   { return m_EstimatedBiasFieldCoefficients; }
00380 
00384   void SetTissueClassStatistics(const Array< double > & means,
00385                                 const Array< double > & sigmas)
00386   throw ( ExceptionObject );
00387 
00389   itkSetMacro(VolumeCorrectionMaximumIteration, int);
00390   itkGetConstMacro(VolumeCorrectionMaximumIteration, int);
00391   itkSetMacro(InterSliceCorrectionMaximumIteration, int);
00392   itkGetConstMacro(InterSliceCorrectionMaximumIteration, int);
00394 
00396   void SetOptimizerInitialRadius(double initRadius)
00397   { m_OptimizerInitialRadius = initRadius; }
00398   double GetOptimizerInitialRadius()
00399   { return m_OptimizerInitialRadius; }
00401 
00403   itkSetMacro(OptimizerGrowthFactor, double);
00404   itkGetConstMacro(OptimizerGrowthFactor, double);
00406 
00409   itkSetMacro(OptimizerShrinkFactor, double);
00410   itkGetConstMacro(OptimizerShrinkFactor, double);
00411 
00418   void SetNumberOfLevels(unsigned int num);
00419 
00421   itkGetConstMacro(NumberOfLevels, unsigned int);
00422 
00429   void SetSchedule(const ScheduleType & schedule);
00430 
00432   itkGetConstReferenceMacro(Schedule, ScheduleType);
00433 
00438   void SetStartingShrinkFactors(unsigned int factor);
00439 
00440   void SetStartingShrinkFactors(unsigned int *factors);
00441 
00443   const unsigned int * GetStartingShrinkFactors() const;
00444 
00448   static bool IsScheduleDownwardDivisible(const ScheduleType & schedule);
00449 
00456   void Initialize()
00457   throw ( ExceptionObject );
00458 
00461   BiasFieldType EstimateBiasField(InputImageRegionType region,
00462                                   unsigned int degree,
00463                                   int maximumIteration);
00464 
00468   void CorrectImage(BiasFieldType & bias,
00469                     InputImageRegionType region);
00470 
00473   void CorrectInterSliceIntensityInhomogeneity(InputImageRegionType region);
00474 
00475 protected:
00476   MRIBiasFieldCorrectionFilter();
00477   virtual ~MRIBiasFieldCorrectionFilter();
00478   void PrintSelf(std::ostream & os, Indent indent) const;
00479 
00482   bool CheckMaskImage(ImageMaskType *mask);
00483 
00484 protected:
00488   void Log1PImage(InternalImageType *source,
00489                   InternalImageType *target);
00490 
00493   void ExpImage(InternalImageType *source,
00494                 InternalImageType *target);
00495 
00498   template< class TSource, class TTarget >
00499   void CopyAndConvertImage(const TSource *source,
00500                            TTarget *target,
00501                            typename TTarget::RegionType requestedRegion)
00502   {
00503     typedef ImageRegionConstIterator< TSource > SourceIterator;
00504     typedef ImageRegionIterator< TTarget >      TargetIterator;
00505     typedef typename TTarget::PixelType         TargetPixelType;
00506 
00507     SourceIterator s_iter(source, requestedRegion);
00508     TargetIterator t_iter(target, requestedRegion);
00509 
00510     s_iter.GoToBegin();
00511     t_iter.GoToBegin();
00512     while ( !s_iter.IsAtEnd() )
00513       {
00514       t_iter.Set( static_cast< TargetPixelType >( s_iter.Get() ) );
00515       ++s_iter;
00516       ++t_iter;
00517       }
00518   }
00519 
00524   void GetBiasFieldSize(InputImageRegionType region,
00525                         BiasFieldType::DomainSizeType & domainSize);
00526 
00530   void AdjustSlabRegions(SlabRegionVectorType & slabs,
00531                          OutputImageRegionType requestedRegion);
00532 
00533   void GenerateData();
00534 
00535 private:
00536   MRIBiasFieldCorrectionFilter(const Self &); //purposely not implemented
00537   void operator=(const Self &);               //purposely not implemented
00538 
00540   EnergyFunctionPointer m_EnergyFunction;
00541 
00543   NormalVariateGeneratorType::Pointer m_NormalVariateGenerator;
00544 
00546   ImageMaskPointer m_InputMask;
00547 
00549   ImageMaskPointer m_OutputMask;
00550 
00552   InternalImagePointer m_InternalInput;
00553 
00555   SlabRegionVectorType m_Slabs;
00556 
00558   int m_SlicingDirection;
00559 
00561   bool m_BiasMultiplicative;
00562 
00564   bool m_UsingInterSliceIntensityCorrection;
00565   bool m_UsingSlabIdentification;
00566   bool m_UsingBiasFieldCorrection;
00567   bool m_GeneratingOutput;
00568 
00569   unsigned int        m_SlabNumberOfSamples;
00570   InputImagePixelType m_SlabBackgroundMinimumThreshold;
00571   double              m_SlabTolerance;
00573   int m_BiasFieldDegree;
00574 
00576   unsigned int m_NumberOfLevels;
00577 
00579   ScheduleType m_Schedule;
00580 
00583   BiasFieldType::CoefficientArrayType m_BiasFieldCoefficients;
00584 
00587   BiasFieldType::CoefficientArrayType m_EstimatedBiasFieldCoefficients;
00588 
00590   int m_VolumeCorrectionMaximumIteration;
00591 
00593   int m_InterSliceCorrectionMaximumIteration;
00594 
00596   double m_OptimizerInitialRadius;
00597 
00599   double m_OptimizerGrowthFactor;
00600 
00602   double m_OptimizerShrinkFactor;
00603 
00605   Array< double > m_TissueClassMeans;
00606 
00608   Array< double > m_TissueClassSigmas;
00609 };
00610 
00611 // ==================================
00612 } // end namespace itk
00613 
00614 #ifndef ITK_MANUAL_INSTANTIATION
00615 #include "itkMRIBiasFieldCorrectionFilter.hxx"
00616 #endif
00617 
00618 #endif
00619