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 "itkImageToImageFilter.h"
00021 #include "itkImage.h"
00022
00023 #include "itkMRASlabIdentifier.h"
00024 #include "itkCompositeValleyFunction.h"
00025 #include "itkMultivariateLegendrePolynomial.h"
00026 #include "Statistics/itkNormalVariateGenerator.h"
00027 #include "itkOnePlusOneEvolutionaryOptimizer.h"
00028 #include "itkArray.h"
00029 #include "itkImageRegionConstIterator.h"
00030
00031
00032 namespace itk
00033 {
00045 template<class TImage, class TImageMask, class TBiasField>
00046 class MRIBiasEnergyFunction : public SingleValuedCostFunction
00047 {
00048 public:
00050 typedef MRIBiasEnergyFunction Self;
00051 typedef SingleValuedCostFunction Superclass;
00052 typedef SmartPointer<Self> Pointer;
00053 typedef SmartPointer<const Self> ConstPointer;
00054
00056 itkTypeMacro( SingleValuedCostFunction, CostFunction );
00057
00059 itkNewMacro(Self);
00060
00062 typedef TImage ImageType ;
00063 typedef TImageMask MaskType ;
00064 typedef typename ImageType::Pointer ImagePointer ;
00065 typedef typename MaskType::Pointer MaskPointer ;
00066 typedef typename ImageType::PixelType ImageElementType ;
00067 typedef typename MaskType::PixelType MaskElementType ;
00068 typedef typename ImageType::IndexType ImageIndexType ;
00069 typedef typename ImageType::RegionType ImageRegionType ;
00070
00072 typedef TBiasField BiasFieldType;
00073
00076 typedef typename Superclass::ParametersType ParametersType ;
00077
00079 typedef Superclass::DerivativeType DerivativeType;
00080
00082 typedef Superclass::MeasureType MeasureType;
00083
00085 itkStaticConstMacro(SpaceDimension, unsigned int, 3);
00086
00088 typedef CompositeValleyFunction InternalEnergyFunction ;
00089
00091 itkSetObjectMacro( Image, ImageType );
00092
00094 itkSetObjectMacro( Mask, MaskType );
00095
00097 itkSetMacro( Region, ImageRegionType );
00098
00100 void SetBiasField(BiasFieldType* bias)
00101 { m_BiasField = bias ; }
00102
00105 double GetEnergy0(double diff)
00106 { return (*m_InternalEnergyFunction)(diff); }
00107
00110 MeasureType GetValue(const ParametersType & parameters ) const ;
00111
00114 void GetDerivative( const ParametersType & itkNotUsed(parameters),
00115 DerivativeType & itkNotUsed(derivative) ) const
00116 { }
00117
00122 void InitializeDistributions( Array<double> classMeans,
00123 Array<double> classSigmas );
00124
00125 unsigned int GetNumberOfParameters(void) const;
00126
00127 private:
00129 BiasFieldType * m_BiasField ;
00130
00132 ImagePointer m_Image ;
00133
00135 MaskPointer m_Mask ;
00136
00138 ImageRegionType m_Region ;
00139
00141 InternalEnergyFunction* m_InternalEnergyFunction ;
00142
00143 protected:
00145 MRIBiasEnergyFunction();
00146
00148 virtual ~MRIBiasEnergyFunction();
00149
00150
00151 private:
00152
00153 MRIBiasEnergyFunction(const Self&);
00154 void operator=(const Self&);
00155
00156 } ;
00157
00158
00159
00197 template <class TInputImage, class TOutputImage, class TMaskImage>
00198 class ITK_EXPORT MRIBiasFieldCorrectionFilter :
00199 public ImageToImageFilter< TInputImage, TOutputImage >
00200 {
00201 public:
00203 typedef MRIBiasFieldCorrectionFilter Self;
00204 typedef ImageToImageFilter< TInputImage, TOutputImage > Superclass;
00205 typedef SmartPointer<Self> Pointer;
00206 typedef SmartPointer<const Self> ConstPointer;
00207
00209 itkNewMacro(Self);
00210
00212 itkTypeMacro(MRIBiasFieldCorrectionFilter, ImageToImageFilter);
00213
00215 itkStaticConstMacro(ImageDimension, unsigned int,
00216 TOutputImage::ImageDimension);
00217
00219 typedef TOutputImage OutputImageType ;
00220 typedef TInputImage InputImageType ;
00221 typedef typename TOutputImage::Pointer OutputImagePointer ;
00222 typedef typename TOutputImage::IndexType OutputImageIndexType ;
00223 typedef typename TOutputImage::PixelType OutputImagePixelType ;
00224 typedef typename TOutputImage::SizeType OutputImageSizeType;
00225 typedef typename TOutputImage::RegionType OutputImageRegionType;
00226 typedef typename TInputImage::Pointer InputImagePointer ;
00227 typedef typename TInputImage::IndexType InputImageIndexType;
00228 typedef typename TInputImage::PixelType InputImagePixelType;
00229 typedef typename TInputImage::SizeType InputImageSizeType;
00230 typedef typename TInputImage::RegionType InputImageRegionType;
00231
00233 typedef TMaskImage ImageMaskType ;
00234 typedef typename ImageMaskType::Pointer ImageMaskPointer ;
00235 typedef typename ImageMaskType::RegionType ImageMaskRegionType ;
00236
00238 typedef Image< float, itkGetStaticConstMacro(ImageDimension) > InternalImageType ;
00239 typedef typename InternalImageType::PixelType InternalImagePixelType ;
00240 typedef typename InternalImageType::Pointer InternalImagePointer ;
00241 typedef typename InternalImageType::RegionType InternalImageRegionType ;
00242
00244 typedef MRASlabIdentifier<InputImageType> MRASlabIdentifierType;
00245 typedef typename MRASlabIdentifierType::SlabRegionVectorType
00246 SlabRegionVectorType;
00247 typedef typename SlabRegionVectorType::iterator SlabRegionVectorIteratorType;
00248
00250 typedef MultivariateLegendrePolynomial BiasFieldType;
00251
00253 typedef MRIBiasEnergyFunction<InternalImageType, ImageMaskType, BiasFieldType>
00254 EnergyFunctionType;
00255 typedef typename EnergyFunctionType::Pointer EnergyFunctionPointer;
00256
00258 typedef Statistics::NormalVariateGenerator NormalVariateGeneratorType ;
00259
00261 typedef OnePlusOneEvolutionaryOptimizer OptimizerType ;
00262
00266 void SetInputMask(ImageMaskType* inputMask);
00267 itkGetObjectMacro( InputMask, ImageMaskType );
00268
00271 void SetOutputMask(ImageMaskType* outputMask) ;
00272
00274 itkGetObjectMacro( OutputMask, ImageMaskType );
00275
00279 void IsBiasFieldMultiplicative(bool flag)
00280 { m_BiasMultiplicative = flag ; }
00281
00283 bool IsBiasFieldMultiplicative()
00284 { return m_BiasMultiplicative ; }
00285
00289 itkSetMacro( UsingInterSliceIntensityCorrection, bool );
00290 itkGetConstReferenceMacro( UsingInterSliceIntensityCorrection, bool );
00291
00297 itkSetMacro( UsingSlabIdentification, bool );
00298 itkGetConstReferenceMacro( UsingSlabIdentification, bool );
00299
00300 itkSetMacro( SlabBackgroundMinimumThreshold, InputImagePixelType );
00301 itkGetConstReferenceMacro( SlabBackgroundMinimumThreshold, InputImagePixelType );
00302
00303 itkSetMacro( SlabNumberOfSamples, unsigned int );
00304 itkGetConstReferenceMacro( SlabNumberOfSamples, unsigned int );
00305
00306 itkSetMacro( SlabTolerance, double );
00307 itkGetConstReferenceMacro( SlabTolerance, double );
00308
00314 itkSetMacro( UsingBiasFieldCorrection, bool );
00315 itkGetConstReferenceMacro( UsingBiasFieldCorrection, bool );
00316
00319 itkSetMacro( GeneratingOutput, bool );
00320 itkGetConstReferenceMacro( GeneratingOutput, bool );
00321
00324 itkSetMacro( SlicingDirection , int );
00325
00327 itkSetMacro( BiasFieldDegree, int );
00328 itkGetMacro( BiasFieldDegree, int );
00329
00332 void SetInitialBiasFieldCoefficients
00333 (const BiasFieldType::CoefficientArrayType &coefficients)
00334 { m_BiasFieldCoefficients = coefficients ; }
00335
00339 BiasFieldType::CoefficientArrayType GetEstimatedBiasFieldCoefficients()
00340 { return m_EstimatedBiasFieldCoefficients ; }
00341
00345 void SetTissueClassStatistics(const Array<double> & means,
00346 const Array<double> & sigmas)
00347 throw (ExceptionObject) ;
00348
00350 itkSetMacro( VolumeCorrectionMaximumIteration, int );
00351 itkGetMacro( VolumeCorrectionMaximumIteration, int );
00352 itkSetMacro( InterSliceCorrectionMaximumIteration, int );
00353 itkGetMacro( InterSliceCorrectionMaximumIteration, int );
00354
00356 void SetOptimizerInitialRadius(double initRadius)
00357 { m_OptimizerInitialRadius = initRadius ; }
00358 double GetOptimizerInitialRadius()
00359 { return m_OptimizerInitialRadius ; }
00360
00362 itkSetMacro( OptimizerGrowthFactor, double );
00363 itkGetMacro( OptimizerGrowthFactor, double );
00364
00367 itkSetMacro( OptimizerShrinkFactor, double );
00368 itkGetMacro( OptimizerShrinkFactor, double );
00369
00376 void Initialize() throw (ExceptionObject) ;
00377
00380 BiasFieldType EstimateBiasField(InputImageRegionType region,
00381 unsigned int degree,
00382 int maximumIteration) ;
00383
00387 void CorrectImage(BiasFieldType& bias,
00388 InputImageRegionType region) ;
00389
00392 void CorrectInterSliceIntensityInhomogeneity(InputImageRegionType region) ;
00393
00394 protected:
00395 MRIBiasFieldCorrectionFilter() ;
00396 virtual ~MRIBiasFieldCorrectionFilter() ;
00397 void PrintSelf(std::ostream& os, Indent indent) const;
00398
00401 bool CheckMaskImage(ImageMaskType* mask) ;
00402
00403 protected:
00407 void Log1PImage(InternalImageType* source,
00408 InternalImageType* target) ;
00409
00412 void ExpImage(InternalImageType* source,
00413 InternalImageType* target) ;
00414
00417 template<class TSource, class TTarget>
00418 void CopyAndConvertImage(const TSource * source,
00419 TTarget * target,
00420 typename TTarget::RegionType requestedRegion)
00421 {
00422 typedef ImageRegionConstIterator<TSource> SourceIterator ;
00423 typedef ImageRegionIterator<TTarget> TargetIterator ;
00424 typedef typename TTarget::PixelType TargetPixelType ;
00425
00426 SourceIterator s_iter(source, requestedRegion) ;
00427 TargetIterator t_iter(target, requestedRegion) ;
00428
00429 s_iter.GoToBegin() ;
00430 t_iter.GoToBegin() ;
00431 while (!s_iter.IsAtEnd())
00432 {
00433 t_iter.Set(static_cast<TargetPixelType>( s_iter.Get() ) ) ;
00434 ++s_iter ;
00435 ++t_iter ;
00436 }
00437 }
00438
00443 void GetBiasFieldSize(InputImageRegionType region,
00444 BiasFieldType::DomainSizeType& domainSize) ;
00445
00449 void AdjustSlabRegions(SlabRegionVectorType& slabs,
00450 OutputImageRegionType requestedRegion) ;
00451
00452 void GenerateData() ;
00453
00454 private:
00455 MRIBiasFieldCorrectionFilter(const Self&);
00456 void operator=(const Self&);
00457
00459 EnergyFunctionPointer m_EnergyFunction ;
00460
00462 NormalVariateGeneratorType::Pointer m_NormalVariateGenerator ;
00463
00465 ImageMaskPointer m_InputMask ;
00466
00468 ImageMaskPointer m_OutputMask ;
00469
00471 InternalImagePointer m_InternalInput ;
00472
00474 SlabRegionVectorType m_Slabs ;
00475
00477 int m_SlicingDirection ;
00478
00480 bool m_BiasMultiplicative ;
00481
00483 bool m_UsingInterSliceIntensityCorrection ;
00484 bool m_UsingSlabIdentification ;
00485 bool m_UsingBiasFieldCorrection ;
00486 bool m_GeneratingOutput ;
00487
00488 unsigned int m_SlabNumberOfSamples ;
00489 InputImagePixelType m_SlabBackgroundMinimumThreshold ;
00490 double m_SlabTolerance ;
00492 int m_BiasFieldDegree ;
00493
00496 BiasFieldType::CoefficientArrayType m_BiasFieldCoefficients ;
00497
00500 BiasFieldType::CoefficientArrayType m_EstimatedBiasFieldCoefficients ;
00501
00503 int m_VolumeCorrectionMaximumIteration ;
00504
00506 int m_InterSliceCorrectionMaximumIteration ;
00507
00509 double m_OptimizerInitialRadius ;
00510
00512 double m_OptimizerGrowthFactor ;
00513
00515 double m_OptimizerShrinkFactor ;
00516
00518 Array<double> m_TissueClassMeans ;
00519
00521 Array<double> m_TissueClassSigmas ;
00522 };
00523
00524
00525
00526
00527
00528 }
00529
00530 #ifndef ITK_MANUAL_INSTANTIATION
00531 #include "itkMRIBiasFieldCorrectionFilter.txx"
00532 #endif
00533
00534 #endif