00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #ifndef __itkMRIBiasFieldCorrectionFilter_h
00018 #define __itkMRIBiasFieldCorrectionFilter_h
00019
00020 #include <time.h>
00021
00022 #include "itkImageToImageFilter.h"
00023 #include "itkImage.h"
00024 #include "itkArray2D.h"
00025 #include "itkMRASlabIdentifier.h"
00026 #include "itkCompositeValleyFunction.h"
00027 #include "itkMultivariateLegendrePolynomial.h"
00028 #include "Statistics/itkNormalVariateGenerator.h"
00029 #include "itkOnePlusOneEvolutionaryOptimizer.h"
00030 #include "itkArray.h"
00031 #include "itkImageRegionConstIterator.h"
00032 #include "itkImageRegionIterator.h"
00033
00034 namespace itk
00035 {
00047 template<class TImage, class TImageMask, class TBiasField>
00048 class MRIBiasEnergyFunction : public SingleValuedCostFunction
00049 {
00050 public:
00052 typedef MRIBiasEnergyFunction Self;
00053 typedef SingleValuedCostFunction Superclass;
00054 typedef SmartPointer<Self> Pointer;
00055 typedef SmartPointer<const Self> ConstPointer;
00056
00058 itkTypeMacro( MRIBiasEnergyFunction, SingleValuedCostFunction );
00059
00061 itkNewMacro(Self);
00062
00064 typedef TImage ImageType;
00065 typedef typename ImageType::Pointer ImagePointer;
00066 typedef typename ImageType::PixelType ImageElementType;
00067 typedef typename ImageType::IndexType ImageIndexType;
00068 typedef typename ImageType::RegionType ImageRegionType;
00069 typedef TImageMask MaskType;
00070 typedef typename MaskType::Pointer MaskPointer;
00071 typedef typename MaskType::PixelType MaskElementType;
00072
00074 typedef TBiasField BiasFieldType;
00075
00078 typedef typename Superclass::ParametersType ParametersType;
00079
00081 typedef Superclass::DerivativeType DerivativeType;
00082
00084 typedef Superclass::MeasureType MeasureType;
00085
00086 itkStaticConstMacro(SpaceDimension, unsigned int, 3);
00087
00089 typedef CompositeValleyFunction InternalEnergyFunction;
00090
00092 typedef unsigned int SamplingFactorType[SpaceDimension];
00093
00095 itkSetObjectMacro( Image, ImageType );
00096
00098 itkSetObjectMacro( Mask, MaskType );
00099
00101 itkSetMacro( Region, ImageRegionType );
00102
00104 void SetBiasField(BiasFieldType* bias)
00105 { m_BiasField = bias; }
00106
00109 void SetSamplingFactors(SamplingFactorType factor)
00110 {
00111 for (unsigned int i = 0; i < SpaceDimension; i++)
00112 {
00113 m_SamplingFactor[i] = factor[i];
00114 }
00115 }
00117
00120 double GetEnergy0(double diff)
00121 {
00122 return (*m_InternalEnergyFunction)(diff);
00123 }
00124
00127 MeasureType GetValue(const ParametersType & parameters ) const;
00128
00131 void GetDerivative( const ParametersType & itkNotUsed(parameters),
00132 DerivativeType & itkNotUsed(derivative) ) const
00133 { }
00134
00139 void InitializeDistributions( Array<double> classMeans,
00140 Array<double> classSigmas );
00141
00142 unsigned int GetNumberOfParameters(void) const;
00143
00144 protected:
00146 MRIBiasEnergyFunction();
00147
00149 virtual ~MRIBiasEnergyFunction();
00150
00151
00152 private:
00153
00155 BiasFieldType * m_BiasField;
00156
00158 ImagePointer m_Image;
00159
00161 MaskPointer m_Mask;
00162
00164 ImageRegionType m_Region;
00165
00167 InternalEnergyFunction* m_InternalEnergyFunction;
00168
00170 SamplingFactorType m_SamplingFactor;
00171
00172 MRIBiasEnergyFunction(const Self&);
00173 void operator=(const Self&);
00174
00175 };
00176
00226 template <class TInputImage, class TOutputImage, class TMaskImage>
00227 class ITK_EXPORT MRIBiasFieldCorrectionFilter :
00228 public ImageToImageFilter< TInputImage, TOutputImage >
00229 {
00230 public:
00232 typedef MRIBiasFieldCorrectionFilter Self;
00233 typedef ImageToImageFilter< TInputImage, TOutputImage > Superclass;
00234 typedef SmartPointer<Self> Pointer;
00235 typedef SmartPointer<const Self> ConstPointer;
00236
00238 itkNewMacro(Self);
00239
00241 itkTypeMacro(MRIBiasFieldCorrectionFilter, ImageToImageFilter);
00242
00244 itkStaticConstMacro(ImageDimension, unsigned int,
00245 TOutputImage::ImageDimension);
00246
00248 typedef TOutputImage OutputImageType;
00249 typedef typename TOutputImage::Pointer OutputImagePointer;
00250 typedef typename TOutputImage::IndexType OutputImageIndexType;
00251 typedef typename TOutputImage::PixelType OutputImagePixelType;
00252 typedef typename TOutputImage::SizeType OutputImageSizeType;
00253 typedef typename TOutputImage::RegionType OutputImageRegionType;
00254
00255 typedef TInputImage InputImageType;
00256 typedef typename TInputImage::Pointer InputImagePointer;
00257 typedef typename TInputImage::IndexType InputImageIndexType;
00258 typedef typename TInputImage::PixelType InputImagePixelType;
00259 typedef typename TInputImage::SizeType InputImageSizeType;
00260 typedef typename TInputImage::RegionType InputImageRegionType;
00261
00263 typedef TMaskImage ImageMaskType;
00264 typedef typename ImageMaskType::Pointer ImageMaskPointer;
00265 typedef typename ImageMaskType::RegionType ImageMaskRegionType;
00266
00268 typedef Image< float, itkGetStaticConstMacro(ImageDimension) >
00269 InternalImageType;
00270 typedef typename InternalImageType::PixelType InternalImagePixelType;
00271 typedef typename InternalImageType::Pointer InternalImagePointer;
00272 typedef typename InternalImageType::RegionType InternalImageRegionType;
00274
00276 typedef MRASlabIdentifier<InputImageType> MRASlabIdentifierType;
00277 typedef typename MRASlabIdentifierType::SlabRegionVectorType
00278 SlabRegionVectorType;
00279 typedef typename SlabRegionVectorType::iterator SlabRegionVectorIteratorType;
00280
00282 typedef MultivariateLegendrePolynomial BiasFieldType;
00283
00285 typedef MRIBiasEnergyFunction<InternalImageType,
00286 ImageMaskType,
00287 BiasFieldType> EnergyFunctionType;
00288 typedef typename EnergyFunctionType::Pointer EnergyFunctionPointer;
00289
00291 typedef Statistics::NormalVariateGenerator NormalVariateGeneratorType;
00292
00294 typedef OnePlusOneEvolutionaryOptimizer OptimizerType;
00295
00297 typedef Array2D<unsigned int> ScheduleType;
00298
00302 void SetInputMask(ImageMaskType* inputMask);
00303 itkGetObjectMacro( InputMask, ImageMaskType );
00305
00308 void SetOutputMask(ImageMaskType* outputMask);
00309
00311 itkGetObjectMacro( OutputMask, ImageMaskType );
00312
00316 void IsBiasFieldMultiplicative(bool flag)
00317 { m_BiasMultiplicative = flag; }
00318
00320 bool IsBiasFieldMultiplicative()
00321 { return m_BiasMultiplicative; }
00322
00326 itkSetMacro( UsingInterSliceIntensityCorrection, bool );
00327 itkGetConstReferenceMacro( UsingInterSliceIntensityCorrection, bool );
00329
00335 itkSetMacro( UsingSlabIdentification, bool );
00336 itkGetConstReferenceMacro( UsingSlabIdentification, bool );
00338
00339 itkSetMacro( SlabBackgroundMinimumThreshold, InputImagePixelType );
00340 itkGetConstReferenceMacro( SlabBackgroundMinimumThreshold,
00341 InputImagePixelType );
00342
00343 itkSetMacro( SlabNumberOfSamples, unsigned int );
00344 itkGetConstReferenceMacro( SlabNumberOfSamples, unsigned int );
00345
00346 itkSetMacro( SlabTolerance, double );
00347 itkGetConstReferenceMacro( SlabTolerance, double );
00348
00354 itkSetMacro( UsingBiasFieldCorrection, bool );
00355 itkGetConstReferenceMacro( UsingBiasFieldCorrection, bool );
00357
00360 itkSetMacro( GeneratingOutput, bool );
00361 itkGetConstReferenceMacro( GeneratingOutput, bool );
00363
00366 itkSetMacro( SlicingDirection , int );
00367
00369 itkSetMacro( BiasFieldDegree, int );
00370 itkGetConstMacro( BiasFieldDegree, int );
00372
00375 void SetInitialBiasFieldCoefficients(const
00376 BiasFieldType::CoefficientArrayType
00377 &coefficients)
00378 { this->Modified(); m_BiasFieldCoefficients = coefficients; }
00379
00383 BiasFieldType::CoefficientArrayType GetEstimatedBiasFieldCoefficients()
00384 { return m_EstimatedBiasFieldCoefficients; }
00385
00389 void SetTissueClassStatistics(const Array<double> & means,
00390 const Array<double> & sigmas)
00391 throw (ExceptionObject);
00392
00394 itkSetMacro( VolumeCorrectionMaximumIteration, int );
00395 itkGetConstMacro( VolumeCorrectionMaximumIteration, int );
00396 itkSetMacro( InterSliceCorrectionMaximumIteration, int );
00397 itkGetConstMacro( InterSliceCorrectionMaximumIteration, int );
00399
00401 void SetOptimizerInitialRadius(double initRadius)
00402 { m_OptimizerInitialRadius = initRadius; }
00403 double GetOptimizerInitialRadius()
00404 { return m_OptimizerInitialRadius; }
00406
00408 itkSetMacro( OptimizerGrowthFactor, double );
00409 itkGetConstMacro( OptimizerGrowthFactor, double );
00411
00414 itkSetMacro( OptimizerShrinkFactor, double );
00415 itkGetConstMacro( OptimizerShrinkFactor, double );
00416
00417
00424 void SetNumberOfLevels(unsigned int num);
00425
00427 itkGetConstMacro(NumberOfLevels, unsigned int);
00428
00435 void SetSchedule( const ScheduleType& schedule );
00436
00438 itkGetConstReferenceMacro(Schedule, ScheduleType);
00439
00444 void SetStartingShrinkFactors( unsigned int factor );
00445 void SetStartingShrinkFactors( unsigned int* factors );
00447
00449 const unsigned int * GetStartingShrinkFactors() const;
00450
00454 static bool IsScheduleDownwardDivisible( const ScheduleType& schedule );
00455
00456
00463 void Initialize() throw (ExceptionObject);
00464
00467 BiasFieldType EstimateBiasField(InputImageRegionType region,
00468 unsigned int degree,
00469 int maximumIteration);
00470
00474 void CorrectImage(BiasFieldType& bias,
00475 InputImageRegionType region);
00476
00479 void CorrectInterSliceIntensityInhomogeneity(InputImageRegionType region);
00480
00481 protected:
00482 MRIBiasFieldCorrectionFilter();
00483 virtual ~MRIBiasFieldCorrectionFilter();
00484 void PrintSelf(std::ostream& os, Indent indent) const;
00485
00488 bool CheckMaskImage(ImageMaskType* mask);
00489
00490 protected:
00494 void Log1PImage(InternalImageType* source,
00495 InternalImageType* target);
00496
00499 void ExpImage(InternalImageType* source,
00500 InternalImageType* target);
00501
00504 template<class TSource, class TTarget>
00505 void CopyAndConvertImage(const TSource * source,
00506 TTarget * target,
00507 typename TTarget::RegionType requestedRegion)
00508 {
00509 typedef ImageRegionConstIterator<TSource> SourceIterator;
00510 typedef ImageRegionIterator<TTarget> TargetIterator;
00511 typedef typename TTarget::PixelType TargetPixelType;
00512
00513 SourceIterator s_iter(source, requestedRegion);
00514 TargetIterator t_iter(target, requestedRegion);
00515
00516 s_iter.GoToBegin();
00517 t_iter.GoToBegin();
00518 while (!s_iter.IsAtEnd())
00519 {
00520 t_iter.Set(static_cast<TargetPixelType>( s_iter.Get() ) );
00521 ++s_iter;
00522 ++t_iter;
00523 }
00524 }
00525
00530 void GetBiasFieldSize(InputImageRegionType region,
00531 BiasFieldType::DomainSizeType& domainSize);
00532
00536 void AdjustSlabRegions(SlabRegionVectorType& slabs,
00537 OutputImageRegionType requestedRegion);
00538
00539 void GenerateData();
00540
00541 private:
00542 MRIBiasFieldCorrectionFilter(const Self&);
00543 void operator=(const Self&);
00544
00546 EnergyFunctionPointer m_EnergyFunction;
00547
00549 NormalVariateGeneratorType::Pointer m_NormalVariateGenerator;
00550
00552 ImageMaskPointer m_InputMask;
00553
00555 ImageMaskPointer m_OutputMask;
00556
00558 InternalImagePointer m_InternalInput;
00559
00561 SlabRegionVectorType m_Slabs;
00562
00564 int m_SlicingDirection;
00565
00567 bool m_BiasMultiplicative;
00568
00570 bool m_UsingInterSliceIntensityCorrection;
00571 bool m_UsingSlabIdentification;
00572 bool m_UsingBiasFieldCorrection;
00573 bool m_GeneratingOutput;
00574
00575 unsigned int m_SlabNumberOfSamples;
00576 InputImagePixelType m_SlabBackgroundMinimumThreshold;
00577 double m_SlabTolerance;
00579 int m_BiasFieldDegree;
00580
00582 unsigned int m_NumberOfLevels;
00583
00585 ScheduleType m_Schedule;
00586
00589 BiasFieldType::CoefficientArrayType m_BiasFieldCoefficients;
00590
00593 BiasFieldType::CoefficientArrayType m_EstimatedBiasFieldCoefficients;
00594
00596 int m_VolumeCorrectionMaximumIteration;
00597
00599 int m_InterSliceCorrectionMaximumIteration;
00600
00602 double m_OptimizerInitialRadius;
00603
00605 double m_OptimizerGrowthFactor;
00606
00608 double m_OptimizerShrinkFactor;
00609
00611 Array<double> m_TissueClassMeans;
00612
00614 Array<double> m_TissueClassSigmas;
00615 };
00616
00617
00618
00619
00620
00621 }
00622
00623 #ifndef ITK_MANUAL_INSTANTIATION
00624 #include "itkMRIBiasFieldCorrectionFilter.txx"
00625 #endif
00626
00627 #endif
00628