N4 Bias Field Correction

Overview

The N4 bias field correction algorithm is a popular method for correcting low frequency intensity non-uniformity present in MRI image data known as a bias or gain field. The method has also been successfully applied as flat-field correction in microscopy data. This method assumes a simple parametric model and does not require tissue classification.

This example demonstrates how to use the SimpleITK N4BiasFieldCorrectionImageFilter class. This filter has one required input image which is affected by a bias field we wish to correct. The primary input is required to have a “real” pixel type of either sitkFloat32 or sitkFloat64. Additionally, there is an optional “MaskImage” input which specifies which pixels are used to estimate the bias-field and suppress pixels close to zero. It is recommended that the mask image follows the common conventions for masked images, which is being of pixel type sitkUint8 and having values of 0 and 1 representing the mask. Additionally, the mask image and the main input image must occupy the same physical space to ensure pixel to pixel correspondence.

Some basic parameters for using this algorithm are parsed from the command line in this example. The shrink factor is used to reduce the size and complexity of the image. The N4 algorithm uses a multi-scale optimization approach to compute the bias field. The SetMaximumNumberOfIterations method takes an array used to set the limit of iterations per resolution level, thereby setting both the iterations and the number of scales ( from the length of the array ). The output of the filter is the bias corrected input image.

Code

#!/usr/bin/env python

from __future__ import print_function

import SimpleITK as sitk
import sys
import os

if len(sys.argv) < 2:
    print("Usage: N4BiasFieldCorrection inputImage " +
          "outputImage [shrinkFactor] [maskImage] [numberOfIterations] " +
          "[numberOfFittingLevels]")
    sys.exit(1)

inputImage = sitk.ReadImage(sys.argv[1], sitk.sitkFloat32)
image = inputImage

if len(sys.argv) > 4:
    maskImage = sitk.ReadImage(sys.argv[4], sitk.sitkUint8)
else:
    maskImage = sitk.OtsuThreshold(inputImage, 0, 1, 200)

if len(sys.argv) > 3:
    image = sitk.Shrink(inputImage,
                             [int(sys.argv[3])] * inputImage.GetDimension())
    maskImage = sitk.Shrink(maskImage,
                            [int(sys.argv[3])] * inputImage.GetDimension())

corrector = sitk.N4BiasFieldCorrectionImageFilter()

numberFittingLevels = 4

if len(sys.argv) > 6:
    numberFittingLevels = int(sys.argv[6])

if len(sys.argv) > 5:
    corrector.SetMaximumNumberOfIterations([int(sys.argv[5])]
                                           * numberFittingLevels)

output = corrector.Execute(image, maskImage)


log_bias_field = corrector.GetLogBiasFieldAsImage(inputImage)

output = inputImage / sitk.Exp( log_bias_field )

sitk.WriteImage(output, sys.argv[2])

if ("SITK_NOSHOW" not in os.environ):
    sitk.Show(output, "N4 Corrected")
# Run with:
#
# Rscript --vanilla N4BiasFieldCorrection.R inputImage, outputImage, shrinkFactor, maskImage, numberOfIterations, numberOfFittingLevels
#

library(SimpleITK)

args <- commandArgs( TRUE )

if (length(args) < 2) {
    stop("At least 2 arguments expected - inputImage, outputImage, [shrinkFactor], ",
         "[maskImage], [numberOfIterations], [numberOfFittingLevels]")
}

inputImage <- ReadImage(args[[1]], 'sitkFloat32')
image <- inputImage

if (length( args ) > 4) {
    maskImage <- ReadImage( args[[4]], 'sitkUint8' )
} else {
    maskImage <- OtsuThreshold( inputImage, 0, 1, 200 )
}

if (length( args ) > 3) {
    image <- Shrink( inputImage, rep(strtoi(args[3]), inputImage$GetDimension()) )
    maskImage <- Shrink( maskImage, rep(strtoi(args[3]), inputImage$GetDimension()) )
}

inputImage <- Cast( inputImage, 'sitkFloat32' )

corrector <- N4BiasFieldCorrectionImageFilter()

numberFittingLevels <- 4

if (length ( args ) > 6) {
    numberFittingLevels <- strtoi( args[[6]] )
}

if (length ( args ) > 5) {
    corrector$SetMaximumNumberOfIterations( rep(strtoi( args[[5]], numberFittingLevels)) )
}

output <- corrector$Execute( image, maskImage )

logBiasField <- corrector$GetLogBiasFieldAsImage(inputImage)

output <- inputImage / Exp( logBiasField )

WriteImage(output, args[[2]])