N3 ITK Implementation for MRI Bias Field Correction
We forego theoretical discussions of MRI bias field correction and defer to relevant references [2, 7, 1]. Instead, we discuss our implementation and how it relates to both Sled’s paper [3] and the original N3 public offering. For notational purposes in this article, we denote the MNI N3 implementation as ‘N3MNI’ and the ITK implementation we offer as ‘N4ITK’.

Implemention
The N4ITK implementation is given as a single class itk::N3MRIBiasFieldCorrectionImageFilter. It is derived from the itk::ImageToImageFilter class (as is the related class itk::MRIBiasFieldCorrectionFilter) since its operation takes the MR image (with an associated mask) corrupted by a bias field as input, and outputs the corrected image. For the user that wants to reconstruct the bias field once the algorithm terminates, we demonstrate how that can be accomplished with the additional class itk::BSplineControlPointImageFilter which is included with this submission. Note that it is only needed if the bias field is to be reconstructed after the N3 algorithm terminates.

Algorithmic Overview
The steps for the N3 algorithm are illustrated in Figure 4 of [3]. Initially, the intensities of the input image are transformed into the log space and an initial log bias field of all zeros is instantiated. In N3MNI, an option is given whereby the user can provide an initial bias field estimate but, to keep the options to a minimum, we decided to omit that option. However, given the open-source nature of the code, the ITK user can modify the code according to preference. After initialization, we iterate by alternating between estimating the unbiased log image and estimating the log of the bias field.

Parameters
One of the attractive aspects of the N3 algorithm is the minimal number of parameters available to tune and the relatively good performance achieved with the default parameters which we tried to maintain, where we could, for both N3MNI and [3]. The available parameters are:

• m_MaskLabel (default = 1): The algorithm requires a mask be supplied by the user with the corresponding mask label. According to Sled, mask generation is not crucial and good results can be achieved with a simple scheme like Otsu thresholding.
• m_NumberOfHistogramBins (default = 200): One of the steps of N3 requires intensity profile construction from the intensities of the uncorrected input image and a triangular parzen windowing scheme. The default value is the same as in N3MNI.
• m_WeinerFilterNoise (default = 0.01): Field estimation is performed by deconvolution using a Wiener filter which has an additive noise term to prevent division by zero (see Equation (12) of [3]). This is identical to the noise variable in N3MNI and equal to Z2 in [3].
• m_BiasFieldFullWidthAtHalfMaximum (default = 0.15): A key contribution to N3 is the usage of a simple Gaussian to model the bias field. This variable characterizes that Gaussian and is the same as the FWHM variable in both N3MNI and [3].
• m_MaximumNumberOfIterations (default = 50): Optimization occurs iteratively until the number of iterations exceeds the maximum specified by this variable.
• m_ConvergenceThreshold (default = 0.001): In [3], the authors propose the coefficient of variation between the ratio of subsequent field estimates as the convergence criterion. However, in both N3MNI and N4ITK, the standard deviation of the ratio between subsequent field estimates is used.
• m_SplineOrder (default = 3): A smooth field estimate is produced after each iterative correction using B-splines. In both N3MNI and [3], cubic splines are used. Although any feasible order of spline is available, the default in N4ITK is also cubic.
• m_NumberOfControlPoints Since the field is usually low frequency, by default we set the number of control points to the minimum m_SplineOrder+1.
• m_NumberOfFittingLevels (default = 4): The B-spline fitting algorithm [6] is different from what is used in N3MNI and proposed in [3]. The version we use was already available in ITK as one of our earlier contributions [5] and is not susceptible to ill-conditioned fitting matrices. One of the parameters for that fitting is the number of hierarchical levels to fit where each successive level doubles the B-spline mesh resolution.

Figure 1: (a) Uncorrected image. (b) Mask image. (c) Bias field corrected image. (d) Uncorrected image with the calculated bias field superimposed.

Bias Field Generation
Oftentimes, the user would like to see the calculated bias field. One of the more obvious reasons for this would be when the bias field is calculated on a downsampled image (suggested in [3] and given as an option in N3MNI and included in the testing code). One would then like to reconstruct the bias field to estimate the corrected image in full resolution. Since the B-spline bias field is a continuous object defined by the control point values and spline order, we can reconstruct the bias field of the full resolution image without loss of accuracy. We demonstrate how this is to be done in the test code itkN3MRIBiasFieldCorrectionImageFilterTest.cxx. Note that the control points describe a B-spline scalar field in log space so the itk::ExpImageFilter has to be used after reconstruction.

Test Code
The test code included with this submission, itkN3MRIBiasFieldCorrectionImageFilterTest.cxx, is designed to allow the user to immediately apply the N4ITK classes to their own images. Usage is given as follows:

This class takes the input image, subsamples it according to the optional shrinkFactor option, and creates the bias field corrected output image. Other optional parameters are the maskImage (if not available, one is created using the itk::OtsuThresholdImageFilter), the number of iterations (default = 50), the number of fitting levels (default = 4), and a file name for writing out the resulting bias field.

Sample Results
We demonstrate usage with two MR images—a 2D brain slice and a volume from a hyperpolarized helium-3 image. We use ITK-SNAP to visualize the results.

2D Brain Slice
Figure 1(a) is the uncorrected image used in our 2D brain test. Close inspection demonstrates a darkening in the white matter toward the upper right of the image. This darkening is corrected in Figure 1(c).

3D Hyperpolarized Helium-3 Lung MRI
Figure 2(a) is the uncorrected image used in our 3D helium-3 MR image volume. Close inspection demonstrates a darkening in the white matter toward the upper portion of the given axial slice. This darkening is corrected in Figure 2(c).

Figure 2: (a) Uncorrected image. (b) Mask image. (c) Bias field corrected image. (d) Uncorrected image with the calculated bias field superimposed.

References
[1]  R. G. Boyes, J. L. Gunter, C. Frost, A. L. Janke, T. Yeatman, D. L.G. Hill, M. A. Bernstein, P. M. Thompson, M. W. Weiner, N. Schuff, G. E.       Alexander, R. J. Killiany, C. DeCarli, C. R. Jack, N. C. Fox, and A. D. N. I. Study. Intensity non-uniformity correction using n3 on 3-t       scanners with multichannel phased array coils. Neuroimage, 39(4):1752–1762, Feb 2008.
[2]  Zujun Hou. A review on mr image intensity inhomogeneity correction. Internation Journal of Biomedical Imaging, 2006:1–11, 2006.
[3]  J. G. Sled, A. P. Zijdenbos, and A. C. Evans. A nonparametric method for automatic correction of intensity nonuniformity in MRI data.       IEEE Transactions on Medical Imaging, 17(1):87–97, Feb 1998.
[4]  M. Styner, C. Brechbuhler, G. Szckely, and G. Gerig. Parametric Estimate of Intensity Inhomogeneities Applied to MRI. IEEE Transactions       on Medical Imaging, 19(3):153–165, March 2000.
[5]  N. J. Tustison and J. C. Gee. N-d Ck B-spline scattered data approximation. The Insight Journal, 2005.
[6]  N. J. Tustison and J. C. Gee. Generalized n-d Ck B-spline scattered data approximation with confidence values. In Proc. Third       International Workshop Medical Imaging and Augmented Reality, pages 76–83, 2006.
[7]  U. Vovk, F. Pernus, and B. Likar. A review of methods for correction of intensity inhomogeneity in mri. IEEE Transactions on Medical       Imaging, 26(3):405–421, March 2007.
[8]  A listing of several relevant algorithms compiled by Finn A. Nielsen at the Technical University of Denmark is provided at       http://neuro.imm.dtu.dk/staff/fnielsen/bib/ and clicking on the folder "Nielsen2001BibSegmentation/".
[9]  http://www.bic.mni.mcgill.ca/software/N3/
[10]  A complete discussion of ‘N4’ is provided at hdl.handle.net/10380/3053.

  Nick Tustison borders his moments of software-writing serenity at the Penn Image Computing and Science Lab (PICSL) with attempts at integrating Heidegger’s notion of “Dasein” with the open source software paradigm---oh, and battling ninjas, too. .