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.

Example Run

Running the Python code with the following inputs:

main("BrainProtonDensitySlice.png", "N4_corrected.nrrd", 1)

produces the images below.

Input Images

Original Input

Output Images

Default Mask Image

Log Bias Image

N4 Bias Corrected Image

Code

C++C#PythonR

#include <SimpleITK.h>
#include <iostream>
#include <stdlib.h>


namespace sitk = itk::simple;

int
main(int argc, char * argv[])
{

  if (argc < 2)
  {
    std::cerr << "Usage: N4BiasFieldCorrection inputImage outputImage";
    std::cerr << " [shrinkFactor] [maskImage] [numberOfIterations]";
    std::cerr << " [numberOfFittingLevels]\n";
    return 1;
  }

  sitk::Image inputImage = sitk::ReadImage(argv[1], sitk::sitkFloat32);
  sitk::Image image = inputImage;

  sitk::Image maskImage;
  if (argc > 4)
  {
    maskImage = sitk::ReadImage(argv[4], sitk::sitkUInt8);
  }
  else
  {
    maskImage = sitk::OtsuThreshold(image, 0, 1, 200);
  }

  unsigned int shrinkFactor = 1;
  if (argc > 3)
  {
    shrinkFactor = atoi(argv[3]);
    std::vector<unsigned int> shrink(inputImage.GetDimension(), shrinkFactor);
    image = sitk::Shrink(inputImage, shrink);
    maskImage = sitk::Shrink(maskImage, shrink);
  }

  sitk::N4BiasFieldCorrectionImageFilter corrector;

  unsigned int numFittingLevels = 4;

  if (argc > 6)
  {
    numFittingLevels = atoi(argv[6]);
  }

  if (argc > 5)
  {
    unsigned int              it = atoi(argv[5]);
    std::vector<unsigned int> iterations(numFittingLevels, it);
    corrector.SetMaximumNumberOfIterations(iterations);
  }

  sitk::Image corrected_image = corrector.Execute(image, maskImage);

  sitk::Image log_bias_field = corrector.GetLogBiasFieldAsImage(inputImage);

  sitk::Image corrected_image_full_resolution = sitk::Divide(inputImage, sitk::Exp(log_bias_field));

  sitk::WriteImage(corrected_image_full_resolution, argv[2]);

  if (shrinkFactor > 1)
  {
    sitk::WriteImage(corrected_image, "CXX-Example-N4BiasFieldCorrection-shrunk.nrrd");
  }

  if (getenv("SITK_NOSHOW") == NULL)
    sitk::Show(corrected_image, "N4 Corrected");

  return 0;
}