ITK
4.1.0
Insight Segmentation and Registration Toolkit
|
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