CBICA Toolkit  1.0.0
itkN3MRIBiasFieldCorrectionImageFilter.h
1 /*=========================================================================
2 
3  Program: Insight Segmentation & Registration Toolkit
4  Module: $RCSfile: itkN3MRIBiasFieldCorrectionImageFilter.h,v $
5  Language: C++
6  Date: $Date: 2009/06/09 16:22:05 $
7  Version: $Revision: 1.6 $
8 
9  Copyright (c) Insight Software Consortium. All rights reserved.
10  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
11 
12  This software is distributed WITHOUT ANY WARRANTY; without even
13  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14  PURPOSE. See the above copyright notices for more information.
15 
16 =========================================================================*/
17 #ifndef __itkN3MRIBiasFieldCorrectionImageFilter_h
18 #define __itkN3MRIBiasFieldCorrectionImageFilter_h
19 
20 #include "itkImageToImageFilter.h"
21 
22 #include "itkBSplineScatteredDataPointSetToImageFilter.h"
23 #include "itkPointSet.h"
24 #include "itkSingleValuedCostFunction.h"
25 #include "itkVector.h"
26 
27 #include "vnl/vnl_vector.h"
28 #include "vcl_complex.h"
29 #include "vnl/vnl_complex_traits.h"
30 
31 namespace itk {
32 
33 /** \class N3MRIBiasFieldCorrectionImageFilter.h
34  * \brief Implementation of the N3 MRI bias field correction algorithm.
35  *
36  * The nonparametric nonuniform intensity normalization (N3) algorithm
37  * is a method for correcting nonuniformity associated with MR images.
38  * The algorithm assumes a simple parametric model (Gaussian) for the bias field
39  * but does not require tissue class segmentation. In addition, there are
40  * only a couple of parameters to tune with the default values performing
41  * quite well.
42  *
43  * N3 has been publicly available as a set of perl scripts
44  * (http://www.bic.mni.mcgill.ca/software/N3/) but, with this class, has been
45  * reimplemented for the ITK library with only one minor variation involving
46  * the b-spline fitting routine. We replaced the original fitting approach
47  * with the itkBSplineScatteredDataPointSetToImageFilter which is not
48  * susceptible to ill-conditioned matrix calculation as is the original proposed
49  * fitting component.
50  *
51  * Notes for the user:
52  * 0. Based on our experience with the filter, sometimes the scale of the bias
53  * field is too small particularly for large mesh element sizes. Therefore,
54  * we added an option (which is true by default) for finding the optimal
55  * scaling of the bias field which is based on minimizing the coefficient
56  * of variation of the corrected image. This cost function is given below
57  * in N3BiasFieldScaleCostFunction.
58  * 1. Since much of the image manipulation is done in the log space of the
59  * intensities, input images with negative and small values (< 1) are
60  * discouraged.
61  * 2. The original authors recommend performing the bias field correction
62  * on a downsampled version of the original image.
63  * 3. A mask and/or confidence image can be supplied.
64  * 4. The filter returns the corrected image. If the bias field is wanted, one
65  * can reconstruct it using the class itkBSplineControlPointImageFilter.
66  * See the IJ article and the test file for an example.
67  * 5. The 'Z' parameter in Sled's 1998 paper is the square root
68  * of the class variable 'm_WeinerFilterNoise'.
69  *
70  * \author Nicholas J. Tustison
71  *
72  * Contributed by Nicholas J. Tustison, James C. Gee
73  * in the Insight Journal paper:
74  *
75  * \par REFERENCE
76  * J.G. Sled, P. Zijdenbos and A.C. Evans, "A comparison of retrospective
77  * intensity non-uniformity correction methods for MRI". Information
78  * Processing Medical Imaging Proc, 15th Int. Conf. IMPI'97, vol.1230,
79  * pp 459-464,1997.
80  *
81  * J.G. Sled, A.P. Zijdenbos and A.C. Evans. "A Nonparametric Method for
82  * Automatic Correction of Intensity Nonuniformity in MRI Data"
83  * IEEE Transactions on Medical Imaging, Vol 17, No 1. Feb 1998.
84  *
85  */
86 
87 
88  /**
89  * Class definition for N3MRIBiasFieldCorrectionImageFilter
90  */
91  template< typename TInputImage, typename TMaskImage = Image<unsigned char,
92  TInputImage::ImageDimension>, typename TOutputImage = TInputImage>
94  public ImageToImageFilter<TInputImage, TOutputImage>
95  {
96  public:
97  /** Standard class typedefs. */
99  typedef ImageToImageFilter<TInputImage, TOutputImage> Superclass;
100  typedef SmartPointer<Self> Pointer;
101  typedef SmartPointer<const Self> ConstPointer;
102 
103  /** Runtime information support. */
105 
106  /** Standard New method. */
107  itkNewMacro(Self);
108 
109  /** ImageDimension constants */
110  itkStaticConstMacro(ImageDimension, unsigned int,
111  TInputImage::ImageDimension);
112 
113  /** Some convenient typedefs. */
114  typedef TInputImage InputImageType;
115  typedef TOutputImage OutputImageType;
116  typedef TMaskImage MaskImageType;
117  typedef typename MaskImageType::PixelType MaskPixelType;
118 
119  typedef float RealType;
120  typedef Image<RealType, ImageDimension> RealImageType;
121 
122  /** B-spline smoothing filter typedefs */
123  typedef Vector<RealType, 1> ScalarType;
124  typedef PointSet<ScalarType,
125  itkGetStaticConstMacro(ImageDimension)> PointSetType;
126  typedef Image<ScalarType,
127  itkGetStaticConstMacro(ImageDimension)> ScalarImageType;
128  typedef BSplineScatteredDataPointSetToImageFilter
129  <PointSetType, ScalarImageType> BSplineFilterType;
130  typedef typename
131  BSplineFilterType::PointDataImageType BiasFieldControlPointLatticeType;
132  typedef typename BSplineFilterType::ArrayType ArrayType;
133 
134  void SetMaskImage(const MaskImageType *mask)
135  {
136  this->SetNthInput(1, const_cast<MaskImageType *>(mask));
137  }
138  const MaskImageType* GetMaskImage() const
139  {
140  return static_cast<MaskImageType*>(const_cast<DataObject *>
141  (this->ProcessObject::GetInput(1)));
142  }
143  void SetConfidenceImage(const RealImageType *image)
144  {
145  this->SetNthInput(2, const_cast<RealImageType *>(image));
146  }
147  const RealImageType* GetConfidenceImage() const
148  {
149  return static_cast<RealImageType*>(const_cast<DataObject *>
150  (this->ProcessObject::GetInput(2)));
151  }
152 
153  void SetInput1(const TInputImage *input)
154  {
155  this->SetInput(input);
156  }
157  void SetInput2(const TMaskImage *mask)
158  {
159  this->SetMaskImage(mask);
160  }
161  void SetInput3(const RealImageType *image)
162  {
163  this->SetConfidenceImage(image);
164  }
165 
166  itkSetMacro(MaskLabel, MaskPixelType);
167  itkGetConstMacro(MaskLabel, MaskPixelType);
168 
169  itkSetMacro(NumberOfHistogramBins, unsigned int);
170  itkGetConstMacro(NumberOfHistogramBins, unsigned int);
171 
172  itkSetMacro(WeinerFilterNoise, RealType);
173  itkGetConstMacro(WeinerFilterNoise, RealType);
174 
175  itkSetMacro(BiasFieldFullWidthAtHalfMaximum, RealType);
176  itkGetConstMacro(BiasFieldFullWidthAtHalfMaximum, RealType);
177 
178  itkSetMacro(MaximumNumberOfIterations, unsigned int);
179  itkGetConstMacro(MaximumNumberOfIterations, unsigned int);
180 
181  itkSetMacro(ConvergenceThreshold, RealType);
182  itkGetConstMacro(ConvergenceThreshold, RealType);
183 
184  itkSetMacro(SplineOrder, unsigned int);
185  itkGetConstMacro(SplineOrder, unsigned int);
186 
187  itkSetMacro(NumberOfFittingLevels, ArrayType);
188  itkGetConstMacro(NumberOfFittingLevels, ArrayType);
189  void SetNumberOfFittingLevels(unsigned int n)
190  {
191  ArrayType nlevels;
192  nlevels.Fill(n);
193  this->SetNumberOfFittingLevels(nlevels);
194  }
195 
196  itkSetMacro(NumberOfControlPoints, ArrayType);
197  itkGetConstMacro(NumberOfControlPoints, ArrayType);
198 
199  itkGetConstMacro(LogBiasFieldControlPointLattice,
200  typename BiasFieldControlPointLatticeType::Pointer);
201 
202  itkSetMacro(UseOptimalBiasFieldScaling, bool);
203  itkGetConstMacro(UseOptimalBiasFieldScaling, bool);
204  itkBooleanMacro(UseOptimalBiasFieldScaling);
205 
206  itkGetConstMacro(BiasFieldScaling, RealType);
207 
208  itkGetConstMacro(ElapsedIterations, unsigned int);
209  itkGetConstMacro(CurrentConvergenceMeasurement, RealType);
210 
211  protected:
212  N3MRIBiasFieldCorrectionImageFilter();
213  ~N3MRIBiasFieldCorrectionImageFilter() {};
214  void PrintSelf(std::ostream& os, Indent indent) const;
215 
216  void GenerateData();
217 
218  private:
219  N3MRIBiasFieldCorrectionImageFilter(const Self&); //purposely not implemented
220  void operator=(const Self&); //purposely not implemented
221 
222  typename RealImageType::Pointer SharpenImage(
223  typename RealImageType::Pointer);
224  typename RealImageType::Pointer SmoothField(
225  typename RealImageType::Pointer);
226  RealType CalculateConvergenceMeasurement(
227  typename RealImageType::Pointer,
228  typename RealImageType::Pointer);
229  RealType CalculateOptimalBiasFieldScaling(
230  typename RealImageType::Pointer);
231 
232  MaskPixelType m_MaskLabel;
233 
234  /**
235  * Parameters for deconvolution with Weiner filter
236  */
237  unsigned int m_NumberOfHistogramBins;
238  RealType m_WeinerFilterNoise;
239  RealType m_BiasFieldFullWidthAtHalfMaximum;
240 
241  /**
242  * Convergence parameters
243  */
244  unsigned int m_MaximumNumberOfIterations;
245  unsigned int m_ElapsedIterations;
246  RealType m_ConvergenceThreshold;
247  RealType m_CurrentConvergenceMeasurement;
248 
249  /**
250  * B-spline fitting parameters
251  */
252  typename
253  BiasFieldControlPointLatticeType::Pointer m_LogBiasFieldControlPointLattice;
254  unsigned int m_SplineOrder;
255  ArrayType m_NumberOfControlPoints;
256  ArrayType m_NumberOfFittingLevels;
257 
258  /**
259  * other parameters
260  */
261  RealType m_BiasFieldScaling;
262  bool m_UseOptimalBiasFieldScaling;
263 
264  }; // end of class
265 
266 /**
267  * Class definition for N3BiasFieldScaleCostFunction
268  */
269 
270  template< typename TInputImage, typename TBiasFieldImage, typename TMaskImage, typename TConfidenceImage>
272  : public SingleValuedCostFunction
273 {
274 public:
276  typedef SingleValuedCostFunction Superclass;
277  typedef SmartPointer<Self> Pointer;
278  typedef SmartPointer<const Self> ConstPointer;
279 
280  /** Run-time type information (and related methods). */
281  itkTypeMacro( N3BiasFieldScaleCostFunction, SingleValuedCostFunction );
282 
283  /** Method for creation through the object factory. */
284  itkNewMacro( Self );
285 
286  typedef Superclass::MeasureType MeasureType;
287  typedef Superclass::DerivativeType DerivativeType;
288  typedef Superclass::ParametersType ParametersType;
289 
290  itkSetObjectMacro( InputImage, TInputImage );
291  itkSetObjectMacro( BiasFieldImage, TBiasFieldImage );
292  itkSetObjectMacro( MaskImage, TMaskImage );
293  itkSetObjectMacro( ConfidenceImage, TConfidenceImage );
294 
295  itkSetMacro( MaskLabel, typename TMaskImage::PixelType );
296  itkGetConstMacro( MaskLabel, typename TMaskImage::PixelType );
297 
298  virtual MeasureType GetValue( const ParametersType & parameters ) const;
299  virtual void GetDerivative( const ParametersType & parameters,
300  DerivativeType & derivative ) const;
301  virtual unsigned int GetNumberOfParameters() const;
302 
303 protected:
305  virtual ~N3BiasFieldScaleCostFunction();
306 
307 private:
308  N3BiasFieldScaleCostFunction(const Self&); //purposely not implemented
309  void operator=(const Self&); //purposely not implemented
310 
311  typename TInputImage::Pointer m_InputImage;
312  typename TBiasFieldImage::Pointer m_BiasFieldImage;
313  typename TMaskImage::Pointer m_MaskImage;
314  typename TConfidenceImage::Pointer m_ConfidenceImage;
315 
316  typename TMaskImage::PixelType m_MaskLabel;
317 };
318 
319 
320 } // end namespace itk
321 
322 #ifndef ITK_MANUAL_INSTANTIATION
323 #include "itkN3MRIBiasFieldCorrectionImageFilter.cpp"
324 #endif
325 
326 #endif
TInputImage InputImageType
Some convenient typedefs.
Definition: itkN3MRIBiasFieldCorrectionImageFilter.h:114
Class definition for N3MRIBiasFieldCorrectionImageFilter.
Definition: itkN3MRIBiasFieldCorrectionImageFilter.h:93
Class definition for N3BiasFieldScaleCostFunction.
Definition: itkN3MRIBiasFieldCorrectionImageFilter.h:271
itkTypeMacro(N3MRIBiasFieldCorrectionImageFilter, ImageToImageFilter)
Runtime information support.
itkNewMacro(Self)
Method for creation through the object factory.
N3MRIBiasFieldCorrectionImageFilter Self
Standard class typedefs.
Definition: itkN3MRIBiasFieldCorrectionImageFilter.h:98
Vector< RealType, 1 > ScalarType
B-spline smoothing filter typedefs.
Definition: itkN3MRIBiasFieldCorrectionImageFilter.h:123
itkNewMacro(Self)
Standard New method.
itkStaticConstMacro(ImageDimension, unsigned int, TInputImage::ImageDimension)
ImageDimension constants.
itkTypeMacro(N3BiasFieldScaleCostFunction, SingleValuedCostFunction)
Run-time type information (and related methods).