Examples/N4BiasFieldCorrection.cxx

#include "antsUtilities.h"
#include "antsAllocImage.h"
#include "antsCommandLineParser.h"

#include "ReadWriteData.h"

#include "itkBinaryThresholdImageFilter.h"
#include "itkBSplineControlPointImageFilter.h"
#include "itkConstantPadImageFilter.h"
#include "itkExpImageFilter.h"
#include "itkExtractImageFilter.h"
#include "itkImageRegionIterator.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkLabelStatisticsImageFilter.h"
#include "itkN4BiasFieldCorrectionImageFilter.h"
#include "itkShrinkImageFilter.h"
#include "itkTimeProbe.h"

#include <string>
#include <algorithm>
#include <vector>

#include "ANTsVersion.h"

namespace ants
{
template <typename TFilter>
class CommandIterationUpdate final : public itk::Command
{
public:
  using Self = CommandIterationUpdate<TFilter>;
  using Superclass = itk::Command;
  using Pointer = itk::SmartPointer<Self>;
  itkNewMacro(Self);

protected:
  CommandIterationUpdate() = default;

public:
  void
  Execute(itk::Object * caller, const itk::EventObject & event) override
  {
    Execute((const itk::Object *)caller, event);
  }

  void
  Execute(const itk::Object * object, const itk::EventObject & event) override
  {
    const auto * filter = dynamic_cast<const TFilter *>(object);

    if (typeid(event) != typeid(itk::IterationEvent))
    {
      return;
    }
    if (filter->GetElapsedIterations() == 1)
    {
      std::cout << "Current level = " << filter->GetCurrentLevel() + 1 << std::endl;
    }
    std::cout << "  Iteration " << filter->GetElapsedIterations() << " (of "
              << filter->GetMaximumNumberOfIterations()[filter->GetCurrentLevel()] << ").  ";
    std::cout << " Current convergence value = " << filter->GetCurrentConvergenceMeasurement()
              << " (threshold = " << filter->GetConvergenceThreshold() << ")" << std::endl;
  }
};

template <unsigned int ImageDimension>
int
N4(itk::ants::CommandLineParser * parser)
{
  using RealType = float;

  using ImageType = itk::Image<RealType, ImageDimension>;
  typename ImageType::Pointer inputImage = nullptr;

  using MaskImageType = itk::Image<RealType, ImageDimension>;
  typename MaskImageType::Pointer maskImage = nullptr;

  bool                                                       verbose = false;
  typename itk::ants::CommandLineParser::OptionType::Pointer verboseOption = parser->GetOption("verbose");
  if (verboseOption && verboseOption->GetNumberOfFunctions())
  {
    verbose = parser->Convert<bool>(verboseOption->GetFunction(0)->GetName());
  }
  if (verbose)
  {
    std::cout << std::endl << "Running N4 for " << ImageDimension << "-dimensional images." << std::endl << std::endl;
  }

  using CorrecterType = itk::N4BiasFieldCorrectionImageFilter<ImageType, MaskImageType, ImageType>;
  typename CorrecterType::Pointer correcter = CorrecterType::New();

  typename itk::ants::CommandLineParser::OptionType::Pointer inputImageOption = parser->GetOption("input-image");
  if (inputImageOption && inputImageOption->GetNumberOfFunctions())
  {
    std::string inputFile = inputImageOption->GetFunction(0)->GetName();
    ReadImage<ImageType>(inputImage, inputFile.c_str());
  }
  else
  {
    if (verbose)
    {
      std::cerr << "Input image not specified." << std::endl;
    }
    return EXIT_FAILURE;
  }

  /**
   * handle the mask image
   */

  bool isMaskImageSpecified = false;

  typename itk::ants::CommandLineParser::OptionType::Pointer maskImageOption = parser->GetOption("mask-image");
  if (maskImageOption && maskImageOption->GetNumberOfFunctions())
  {
    std::string inputMaskFile = maskImageOption->GetFunction(0)->GetName();
    ReadImage<MaskImageType>(maskImage, inputMaskFile.c_str());

    isMaskImageSpecified = true;
  }
  if (maskImage.IsNull())
  {
    if (verbose)
    {
      std::cout << "Mask not read.  Using the entire image as the mask." << std::endl << std::endl;
    }
    maskImage = MaskImageType::New();
    maskImage->CopyInformation(inputImage);
    maskImage->SetRegions(inputImage->GetRequestedRegion());
    maskImage->Allocate();
    maskImage->FillBuffer(itk::NumericTraits<typename MaskImageType::PixelType>::OneValue());
  }

  /**
   * check for negative values in the masked region
   */
  using ThresholderType = itk::BinaryThresholdImageFilter<MaskImageType, MaskImageType>;
  typename ThresholderType::Pointer thresholder = ThresholderType::New();
  thresholder->SetInsideValue(itk::NumericTraits<typename MaskImageType::PixelType>::ZeroValue());
  thresholder->SetOutsideValue(itk::NumericTraits<typename MaskImageType::PixelType>::OneValue());
  thresholder->SetLowerThreshold(itk::NumericTraits<typename MaskImageType::PixelType>::ZeroValue());
  thresholder->SetUpperThreshold(itk::NumericTraits<typename MaskImageType::PixelType>::ZeroValue());
  thresholder->SetInput(maskImage);

  using StatsType = itk::LabelStatisticsImageFilter<ImageType, MaskImageType>;
  typename StatsType::Pointer statsOriginal = StatsType::New();
  statsOriginal->SetInput(inputImage);
  statsOriginal->SetLabelInput(thresholder->GetOutput());
  statsOriginal->UseHistogramsOff();
  statsOriginal->Update();

  using StatsLabelType = typename StatsType::LabelPixelType;
  StatsLabelType maskLabel = itk::NumericTraits<StatsLabelType>::OneValue();

  RealType minOriginal = statsOriginal->GetMinimum(maskLabel);
  RealType maxOriginal = statsOriginal->GetMaximum(maskLabel);

  if (verbose)
  {
    std::cout << "Original intensity range:  [" << minOriginal << ", " << maxOriginal << "]" << std::endl;
  }

  if (minOriginal <= 0)
  {
    if (verbose)
    {
      std::cout << std::endl;
      std::cout << "***********************************************************" << std::endl;
      std::cout << "Warning:  Your input image contains nonpositive values" << std::endl;
      std::cout << "which could cause failure or problematic results.  A" << std::endl;
      std::cout << "possible workaround would be to:" << std::endl;
      std::cout << "   1. rescale your image to positive values e.g., [10,100]." << std::endl;
      std::cout << "   2. run N4 on your rescaled image." << std::endl;
      std::cout << "   3. (optional) rescale the N4 output to the original" << std::endl;
      std::cout << "      intensity range." << std::endl;
      std::cout << "***********************************************************" << std::endl;
      std::cout << std::endl;
    }
  }

  /**
   * handle the weight image
   */

  typename ImageType::Pointer weightImage = nullptr;

  typename itk::ants::CommandLineParser::OptionType::Pointer weightImageOption = parser->GetOption("weight-image");
  if (weightImageOption && weightImageOption->GetNumberOfFunctions())
  {
    std::string inputFile = weightImageOption->GetFunction(0)->GetName();
    ReadImage<ImageType>(weightImage, inputFile.c_str());
  }

  /**
   * convergence options
   */
  typename itk::ants::CommandLineParser::OptionType::Pointer convergenceOption = parser->GetOption("convergence");
  if (convergenceOption && convergenceOption->GetNumberOfFunctions())
  {
    if (convergenceOption->GetFunction(0)->GetNumberOfParameters() > 0)
    {
      std::vector<unsigned int> numIters =
        parser->ConvertVector<unsigned int>(convergenceOption->GetFunction(0)->GetParameter(0));
      typename CorrecterType::VariableSizeArrayType maximumNumberOfIterations(numIters.size());
      for (unsigned int d = 0; d < numIters.size(); d++)
      {
        maximumNumberOfIterations[d] = numIters[d];
      }
      correcter->SetMaximumNumberOfIterations(maximumNumberOfIterations);

      typename CorrecterType::ArrayType numberOfFittingLevels;
      numberOfFittingLevels.Fill(numIters.size());
      correcter->SetNumberOfFittingLevels(numberOfFittingLevels);
      correcter->SetConvergenceThreshold(0.0);
    }
    if (convergenceOption->GetFunction(0)->GetNumberOfParameters() > 1)
    {
      correcter->SetConvergenceThreshold(parser->Convert<float>(convergenceOption->GetFunction(0)->GetParameter(1)));
    }
  }
  else // set default values
  {
    typename CorrecterType::VariableSizeArrayType maximumNumberOfIterations(4);
    maximumNumberOfIterations.Fill(50);
    correcter->SetMaximumNumberOfIterations(maximumNumberOfIterations);
    correcter->SetNumberOfFittingLevels(4);
    correcter->SetConvergenceThreshold(0.0);
  }

  /**
   * B-spline options -- we place this here to take care of the case where
   * the user wants to specify things in terms of the spline distance.
   */

  typename ImageType::IndexType inputImageIndex = inputImage->GetLargestPossibleRegion().GetIndex();
  typename ImageType::SizeType  inputImageSize = inputImage->GetLargestPossibleRegion().GetSize();

  typename ImageType::PointType newOrigin = inputImage->GetOrigin();

  typename itk::ants::CommandLineParser::OptionType::Pointer bsplineOption = parser->GetOption("bspline-fitting");
  if (bsplineOption && bsplineOption->GetNumberOfFunctions())
  {
    if (bsplineOption->GetFunction(0)->GetNumberOfParameters() > 1)
    {
      correcter->SetSplineOrder(parser->Convert<unsigned int>(bsplineOption->GetFunction(0)->GetParameter(1)));
    }
    if (bsplineOption->GetFunction(0)->GetNumberOfParameters() > 0)
    {
      std::vector<float> array = parser->ConvertVector<float>(bsplineOption->GetFunction(0)->GetParameter(0));
      typename CorrecterType::ArrayType numberOfControlPoints;
      if (array.size() == 1)
      {
        // the user wants to specify things in terms of spline distance.
        //  1. need to pad the images to get as close to possible to the
        //     requested domain size.
        float splineDistance = array[0];

        typename ImageType::SizeType originalImageSize = inputImage->GetLargestPossibleRegion().GetSize();

        itk::Size<ImageDimension> lowerBound;
        itk::Size<ImageDimension> upperBound;
        for (unsigned int d = 0; d < ImageDimension; d++)
        {
          float domain = static_cast<float>(originalImageSize[d] - 1) * static_cast<float>(inputImage->GetSpacing()[d]);
          auto  numberOfSpans = static_cast<unsigned int>(std::ceil(domain / splineDistance));
          auto  extraPadding = static_cast<unsigned long>((numberOfSpans * splineDistance - domain) /
                                                           static_cast<float>(inputImage->GetSpacing()[d]) +
                                                         static_cast<float>(0.5));
          lowerBound[d] = static_cast<unsigned long>(0.5 * extraPadding);
          upperBound[d] = extraPadding - lowerBound[d];
          newOrigin[d] -= (static_cast<double>(lowerBound[d]) * inputImage->GetSpacing()[d]);
          numberOfControlPoints[d] = numberOfSpans + correcter->GetSplineOrder();
        }

        using PadderType = itk::ConstantPadImageFilter<ImageType, ImageType>;
        typename PadderType::Pointer padder = PadderType::New();
        padder->SetInput(inputImage);
        padder->SetPadLowerBound(lowerBound);
        padder->SetPadUpperBound(upperBound);
        padder->SetConstant(itk::NumericTraits<typename ImageType::PixelType>::ZeroValue());
        padder->Update();

        inputImage = padder->GetOutput();
        inputImage->DisconnectPipeline();

        using MaskPadderType = itk::ConstantPadImageFilter<MaskImageType, MaskImageType>;
        typename MaskPadderType::Pointer maskPadder = MaskPadderType::New();
        maskPadder->SetInput(maskImage);
        maskPadder->SetPadLowerBound(lowerBound);
        maskPadder->SetPadUpperBound(upperBound);
        maskPadder->SetConstant(0);
        maskPadder->Update();

        maskImage = maskPadder->GetOutput();
        maskImage->DisconnectPipeline();

        if (weightImage)
        {
          typename PadderType::Pointer weightPadder = PadderType::New();
          weightPadder->SetInput(weightImage);
          weightPadder->SetPadLowerBound(lowerBound);
          weightPadder->SetPadUpperBound(upperBound);
          weightPadder->SetConstant(0);
          weightPadder->Update();

          weightImage = weightPadder->GetOutput();
          weightImage->DisconnectPipeline();
        }

        if (verbose)
        {
          std::cout << "Specified spline distance: " << splineDistance << "mm" << std::endl;
          std::cout << "  original image size:  " << originalImageSize << std::endl;
          std::cout << "  padded image size:  " << inputImage->GetLargestPossibleRegion().GetSize() << std::endl;
          std::cout << "  number of control points:  " << numberOfControlPoints << std::endl;
          std::cout << std::endl;
        }
      }
      else if (array.size() == ImageDimension)
      {
        for (unsigned int d = 0; d < ImageDimension; d++)
        {
          numberOfControlPoints[d] = static_cast<unsigned int>(array[d]) + correcter->GetSplineOrder();
        }
      }
      else
      {
        if (verbose)
        {
          std::cerr << "Incorrect mesh resolution" << std::endl;
        }
        return EXIT_FAILURE;
      }
      correcter->SetNumberOfControlPoints(numberOfControlPoints);
    }
  }

  using ShrinkerType = itk::ShrinkImageFilter<ImageType, ImageType>;
  typename ShrinkerType::Pointer shrinker = ShrinkerType::New();
  shrinker->SetInput(inputImage);
  shrinker->SetShrinkFactors(1);

  using MaskShrinkerType = itk::ShrinkImageFilter<MaskImageType, MaskImageType>;
  typename MaskShrinkerType::Pointer maskshrinker = MaskShrinkerType::New();
  maskshrinker->SetInput(maskImage);
  maskshrinker->SetShrinkFactors(1);

  typename itk::ants::CommandLineParser::OptionType::Pointer shrinkFactorOption = parser->GetOption("shrink-factor");
  int                                                        shrinkFactor = 4;
  if (shrinkFactorOption && shrinkFactorOption->GetNumberOfFunctions())
  {
    shrinkFactor = parser->Convert<int>(shrinkFactorOption->GetFunction(0)->GetName());
  }
  shrinker->SetShrinkFactors(shrinkFactor);
  maskshrinker->SetShrinkFactors(shrinkFactor);
  if (ImageDimension == 4)
  {
    shrinker->SetShrinkFactor(3, 1);
    maskshrinker->SetShrinkFactor(3, 1);
  }
  shrinker->Update();
  maskshrinker->Update();

  itk::TimeProbe timer;
  timer.Start();

  correcter->SetInput(shrinker->GetOutput());
  correcter->SetMaskImage(maskshrinker->GetOutput());

  using WeightShrinkerType = itk::ShrinkImageFilter<ImageType, ImageType>;
  typename WeightShrinkerType::Pointer weightshrinker = WeightShrinkerType::New();
  if (weightImage)
  {
    weightshrinker->SetInput(weightImage);
    weightshrinker->SetShrinkFactors(shrinker->GetShrinkFactors());
    weightshrinker->Update();

    correcter->SetConfidenceImage(weightshrinker->GetOutput());
  }

  if (verbose)
  {
    using CommandType = CommandIterationUpdate<CorrecterType>;
    typename CommandType::Pointer observer = CommandType::New();
    correcter->AddObserver(itk::IterationEvent(), observer);
  }

  /**
   * histogram sharpening options
   */
  typename itk::ants::CommandLineParser::OptionType::Pointer histOption = parser->GetOption("histogram-sharpening");
  if (histOption && histOption->GetNumberOfFunctions())
  {
    if (histOption->GetFunction(0)->GetNumberOfParameters() > 0)
    {
      correcter->SetBiasFieldFullWidthAtHalfMaximum(
        parser->Convert<float>(histOption->GetFunction(0)->GetParameter(0)));
    }
    if (histOption->GetFunction(0)->GetNumberOfParameters() > 1)
    {
      correcter->SetWienerFilterNoise(parser->Convert<float>(histOption->GetFunction(0)->GetParameter(1)));
    }
    if (histOption->GetFunction(0)->GetNumberOfParameters() > 2)
    {
      correcter->SetNumberOfHistogramBins(parser->Convert<unsigned int>(histOption->GetFunction(0)->GetParameter(2)));
    }
  }

  try
  {
    // correcter->DebugOn();
    correcter->Update();
  }
  catch (const itk::ExceptionObject & e)
  {
    if (verbose)
    {
      std::cerr << "Exception caught: " << e << std::endl;
    }
    return EXIT_FAILURE;
  }

  if (verbose)
  {
    correcter->Print(std::cout, 3);
  }

  timer.Stop();
  if (verbose)
  {
    std::cout << "Elapsed time: " << timer.GetMean() << std::endl;
  }

  /**
   * output
   */

  typename itk::ants::CommandLineParser::OptionType::Pointer outputOption = parser->GetOption("output");
  if (outputOption && outputOption->GetNumberOfFunctions())
  {
    /**
     * Reconstruct the bias field at full image resolution.  Divide
     * the original input image by the bias field to get the final
     * corrected image.
     */
    using BSplinerType = itk::BSplineControlPointImageFilter<typename CorrecterType::BiasFieldControlPointLatticeType,
                                                             typename CorrecterType::ScalarImageType>;
    typename BSplinerType::Pointer bspliner = BSplinerType::New();
    bspliner->SetInput(correcter->GetLogBiasFieldControlPointLattice());
    bspliner->SetSplineOrder(correcter->GetSplineOrder());
    bspliner->SetSize(inputImage->GetLargestPossibleRegion().GetSize());
    bspliner->SetOrigin(newOrigin);
    bspliner->SetDirection(inputImage->GetDirection());
    bspliner->SetSpacing(inputImage->GetSpacing());
    bspliner->Update();

    typename ImageType::Pointer logField = AllocImage<ImageType>(inputImage);

    itk::ImageRegionIterator<typename CorrecterType::ScalarImageType> ItB(
      bspliner->GetOutput(), bspliner->GetOutput()->GetLargestPossibleRegion());
    itk::ImageRegionIterator<ImageType> ItF(logField, logField->GetLargestPossibleRegion());
    for (ItB.GoToBegin(), ItF.GoToBegin(); !ItB.IsAtEnd(); ++ItB, ++ItF)
    {
      ItF.Set(ItB.Get()[0]);
    }

    using ExpFilterType = itk::ExpImageFilter<ImageType, ImageType>;
    typename ExpFilterType::Pointer expFilter = ExpFilterType::New();
    expFilter->SetInput(logField);
    expFilter->Update();

    using DividerType = itk::DivideImageFilter<ImageType, ImageType, ImageType>;
    typename DividerType::Pointer divider = DividerType::New();
    divider->SetInput1(inputImage);
    divider->SetInput2(expFilter->GetOutput());

    typename ImageType::Pointer dividedImage = divider->GetOutput();
    dividedImage->Update();
    dividedImage->DisconnectPipeline();

    if (maskImageOption && maskImageOption->GetNumberOfFunctions() > 0)
    {
      itk::ImageRegionIteratorWithIndex<ImageType> ItD(dividedImage, dividedImage->GetLargestPossibleRegion());
      itk::ImageRegionIterator<ImageType>          ItI(inputImage, inputImage->GetLargestPossibleRegion());
      for (ItD.GoToBegin(), ItI.GoToBegin(); !ItD.IsAtEnd(); ++ItD, ++ItI)
      {
        if (itk::Math::FloatAlmostEqual(maskImage->GetPixel(ItD.GetIndex()),
                                        itk::NumericTraits<typename MaskImageType::PixelType>::ZeroValue()))
        {
          ItD.Set(ItI.Get());
        }
      }
    }

    bool doRescale = true;

    typename itk::ants::CommandLineParser::OptionType::Pointer rescaleOption = parser->GetOption("rescale-intensities");
    if (!isMaskImageSpecified || (rescaleOption && rescaleOption->GetNumberOfFunctions() &&
                                  !parser->Convert<bool>(rescaleOption->GetFunction()->GetName())))
    {
      doRescale = false;
    }

    if (doRescale)
    {
      typename ThresholderType::Pointer thresholder2 = ThresholderType::New();
      thresholder2->SetInsideValue(itk::NumericTraits<typename MaskImageType::PixelType>::ZeroValue());
      thresholder2->SetOutsideValue(itk::NumericTraits<typename MaskImageType::PixelType>::OneValue());
      thresholder2->SetLowerThreshold(itk::NumericTraits<typename MaskImageType::PixelType>::ZeroValue());
      thresholder2->SetUpperThreshold(itk::NumericTraits<typename MaskImageType::PixelType>::ZeroValue());
      thresholder2->SetInput(maskImage);

      typename StatsType::Pointer statsBiasCorrected = StatsType::New();
      statsBiasCorrected->SetInput(dividedImage);
      statsBiasCorrected->SetLabelInput(thresholder2->GetOutput());
      statsBiasCorrected->UseHistogramsOff();
      statsBiasCorrected->Update();

      RealType minBiasCorrected = statsBiasCorrected->GetMinimum(maskLabel);
      RealType maxBiasCorrected = statsBiasCorrected->GetMaximum(maskLabel);

      RealType slope = (maxOriginal - minOriginal) / (maxBiasCorrected - minBiasCorrected);

      itk::ImageRegionIteratorWithIndex<ImageType> ItD(dividedImage, dividedImage->GetLargestPossibleRegion());
      for (ItD.GoToBegin(); !ItD.IsAtEnd(); ++ItD)
      {
        if (itk::Math::FloatAlmostEqual(maskImage->GetPixel(ItD.GetIndex()), static_cast<RealType>(maskLabel)))
        {
          RealType originalIntensity = ItD.Get();
          RealType rescaledIntensity = maxOriginal - slope * (maxBiasCorrected - originalIntensity);
          ItD.Set(rescaledIntensity);
        }
      }
    }

    typename ImageType::RegionType inputRegion;
    inputRegion.SetIndex(inputImageIndex);
    inputRegion.SetSize(inputImageSize);

    using CropperType = itk::ExtractImageFilter<ImageType, ImageType>;
    typename CropperType::Pointer cropper = CropperType::New();
    cropper->SetInput(dividedImage);
    cropper->SetExtractionRegion(inputRegion);
    cropper->SetDirectionCollapseToSubmatrix();
    cropper->Update();

    typename CropperType::Pointer biasFieldCropper = CropperType::New();
    biasFieldCropper->SetInput(expFilter->GetOutput());
    biasFieldCropper->SetExtractionRegion(inputRegion);
    biasFieldCropper->SetDirectionCollapseToSubmatrix();
    biasFieldCropper->Update();

    if (outputOption->GetFunction(0)->GetNumberOfParameters() == 0)
    {
      ANTs::WriteImage<ImageType>(cropper->GetOutput(), (outputOption->GetFunction(0)->GetName()).c_str());
    }
    else if (outputOption->GetFunction(0)->GetNumberOfParameters() > 0)
    {
      ANTs::WriteImage<ImageType>(cropper->GetOutput(), (outputOption->GetFunction(0)->GetParameter(0)).c_str());
      if (outputOption->GetFunction(0)->GetNumberOfParameters() > 1)
      {
        ANTs::WriteImage<ImageType>(biasFieldCropper->GetOutput(),
                                    (outputOption->GetFunction(0)->GetParameter(1)).c_str());
      }
    }
  }

  return EXIT_SUCCESS;
}

void
N4InitializeCommandLineOptions(itk::ants::CommandLineParser * parser)
{
  using OptionType = itk::ants::CommandLineParser::OptionType;

  {
    std::string description = std::string("This option forces the image to be treated as a specified-") +
                              std::string("dimensional image.  If not specified, N4 tries to ") +
                              std::string("infer the dimensionality from the input image.");
    OptionType::Pointer option = OptionType::New();
    option->SetLongName("image-dimensionality");
    option->SetShortName('d');
    option->SetUsageOption(0, "2/3/4");
    option->SetDescription(description);
    parser->AddOption(option);
  }

  {
    std::string description = std::string("A scalar image is expected as input for bias correction.  ") +
                              std::string("Since N4 log transforms the intensities, negative values ") +
                              std::string("or values close to zero should be processed prior to ") +
                              std::string("correction.");

    OptionType::Pointer option = OptionType::New();
    option->SetLongName("input-image");
    option->SetShortName('i');
    option->SetUsageOption(0, "inputImageFilename");
    option->SetDescription(description);
    parser->AddOption(option);
  }

  {
    std::string description =
      std::string("If a mask image is specified, the final bias correction is ") +
      std::string("only performed in the mask region.  If a weight image is not ") +
      std::string("specified, only intensity values inside the masked region are ") +
      std::string("used during the execution of the algorithm.  If a weight ") +
      std::string("image is specified, only the non-zero weights are used in the ") +
      std::string("execution of the algorithm although the mask region defines ") +
      std::string("where bias correction is performed in the final output. ") +
      std::string("Otherwise bias correction occurs over the entire image domain. ") +
      std::string("See also the option description for the weight image. ") +
      std::string("If a mask image is *not* specified then the entire image region ") +
      std::string("will be used as the mask region.  Note that this is different than ") +
      std::string("the N3 implementation which uses the results of Otsu thresholding ") +
      std::string("to define a mask.  However, this leads to unknown anatomical regions being ") +
      std::string("included and excluded during the bias correction.");

    OptionType::Pointer option = OptionType::New();
    option->SetLongName("mask-image");
    option->SetShortName('x');
    option->SetUsageOption(0, "maskImageFilename");
    option->SetDescription(description);
    parser->AddOption(option);
  }

  {
    std::string description = std::string("At each iteration, a new intensity mapping is calculated ") +
                              std::string("and applied but there is nothing which constrains the ") +
                              std::string("new intensity range to be within certain values.  The ") +
                              std::string("result is that the range can \"drift\" from the original ") +
                              std::string("at each iteration.  This option rescales to the [min,max] ") +
                              std::string("range of the original image intensities within the user-specified mask. ") +
                              std::string("A mask is required to perform rescaling.");

    OptionType::Pointer option = OptionType::New();
    option->SetLongName("rescale-intensities");
    option->SetShortName('r');
    option->SetUsageOption(0, "0/(1)");
    option->SetDescription(description);
    parser->AddOption(option);
  }

  {
    std::string description = std::string("The weight image allows the user to perform a relative ") +
                              std::string("weighting of specific voxels during the B-spline fitting. ") +
                              std::string("For example, some studies have shown that N3 performed on ") +
                              std::string("white matter segmentations improves performance.  If one ") +
                              std::string("has a spatial probability map of the white matter, one can ") +
                              std::string("use this map to weight the b-spline fitting towards those ") +
                              std::string("voxels which are more probabilistically classified as white ") +
                              std::string("matter.  See also the option description for the mask image.");

    OptionType::Pointer option = OptionType::New();
    option->SetLongName("weight-image");
    option->SetUsageOption(0, "weightImageFilename");
    option->SetShortName('w');
    option->SetDescription(description);
    parser->AddOption(option);
  }

  {
    std::string description = std::string("Running N4 on large images can be time consuming. ") +
                              std::string("To lessen computation time, the input image can be resampled. ") +
                              std::string("The shrink factor, specified as a single integer, describes ") +
                              std::string("this resampling.  Shrink factors <= 4 are commonly used.") +
                              std::string("Note that the shrink factor is only applied to the first two or ") +
                              std::string("three dimensions which we assume are spatial.  ");

    OptionType::Pointer option = OptionType::New();
    option->SetLongName("shrink-factor");
    option->SetShortName('s');
    option->SetUsageOption(0, "1/2/3/(4)/...");
    option->SetDescription(description);
    parser->AddOption(option);
  }

  {
    std::string description = std::string("Convergence is determined by calculating the coefficient of ") +
                              std::string("variation between subsequent iterations. When this value ") +
                              std::string("is less than the specified threshold ") +
                              std::string("from the previous iteration or the maximum number of ") +
                              std::string("iterations is exceeded the program terminates.  Multiple ") +
                              std::string("resolutions can be specified by using 'x' between the number ") +
                              std::string("of iterations at each resolution, e.g. 100x50x50.");

    OptionType::Pointer option = OptionType::New();
    option->SetLongName("convergence");
    option->SetShortName('c');
    option->SetUsageOption(0, "[<numberOfIterations=50x50x50x50>,<convergenceThreshold=0.0>]");
    option->SetDescription(description);
    parser->AddOption(option);
  }

  {
    std::string description = std::string("These options describe the b-spline fitting parameters. ") +
                              std::string("The initial b-spline mesh at the coarsest resolution is ") +
                              std::string("specified either as the number of elements in each dimension, ") +
                              std::string("e.g. 2x2x3 for 3-D images, or it can be specified as a ") +
                              std::string("single scalar parameter which describes the isotropic sizing ") +
                              std::string("of the mesh elements. The latter option is typically preferred. ") +
                              std::string("For each subsequent level, the spline distance decreases in ") +
                              std::string("half, or equivalently, the number of mesh elements doubles ") +
                              std::string("Cubic splines (order = 3) are typically used.  The default setting ") +
                              std::string("is to employ a single mesh element over the entire domain, i.e., ") +
                              std::string("-b [1x1x1,3].");

    OptionType::Pointer option = OptionType::New();
    option->SetLongName("bspline-fitting");
    option->SetShortName('b');
    option->SetUsageOption(0, "[splineDistance,<splineOrder=3>]");
    option->SetUsageOption(1, "[initialMeshResolution,<splineOrder=3>]");
    option->SetDescription(description);
    parser->AddOption(option);
  }

  {
    std::string description = std::string("These options describe the histogram sharpening parameters, ") +
                              std::string("i.e. the deconvolution step parameters described in the ") +
                              std::string("original N3 algorithm.  The default values have been shown ") +
                              std::string("to work fairly well.");

    OptionType::Pointer option = OptionType::New();
    option->SetLongName("histogram-sharpening");
    option->SetShortName('t');
    option->SetUsageOption(0, "[<FWHM=0.15>,<wienerNoise=0.01>,<numberOfHistogramBins=200>]");
    option->SetDescription(description);
    parser->AddOption(option);
  }

  {
    std::string description = std::string("The output consists of the bias corrected version of the ") +
                              std::string("input image.  Optionally, one can also output the estimated ") +
                              std::string("bias field.");

    OptionType::Pointer option = OptionType::New();
    option->SetLongName("output");
    option->SetShortName('o');
    option->SetUsageOption(0, "correctedImage");
    option->SetUsageOption(1, "[correctedImage,<biasField>]");
    option->SetDescription(description);
    parser->AddOption(option);
  }

  {
    std::string         description = std::string("Get Version Information.");
    OptionType::Pointer option = OptionType::New();
    option->SetLongName("version");
    option->SetDescription(description);
    parser->AddOption(option);
  }

  {
    std::string description = std::string("Verbose output.");

    OptionType::Pointer option = OptionType::New();
    option->SetShortName('v');
    option->SetLongName("verbose");
    option->SetUsageOption(0, "(0)/1");
    option->SetDescription(description);
    parser->AddOption(option);
  }

  {
    std::string description = std::string("Print the help menu (short version).");

    OptionType::Pointer option = OptionType::New();
    option->SetShortName('h');
    option->SetDescription(description);
    parser->AddOption(option);
  }

  {
    std::string description = std::string("Print the help menu.");

    OptionType::Pointer option = OptionType::New();
    option->SetLongName("help");
    option->SetDescription(description);
    parser->AddOption(option);
  }
}

// entry point for the library; parameter 'args' is equivalent to 'argv' in (argc,argv) of commandline parameters to
// 'main()'
int
N4BiasFieldCorrection(std::vector<std::string> args, std::ostream * /*out_stream = nullptr */)
{
  // put the arguments coming in as 'args' into standard (argc,argv) format;
  // 'args' doesn't have the command name as first, argument, so add it manually;
  // 'args' may have adjacent arguments concatenated into one argument,
  // which the parser should handle
  args.insert(args.begin(), "N4BiasFieldCorrection");

  int     argc = args.size();
  char ** argv = new char *[args.size() + 1];
  for (unsigned int i = 0; i < args.size(); ++i)
  {
    // allocate space for the string plus a null character
    argv[i] = new char[args[i].length() + 1];
    std::strncpy(argv[i], args[i].c_str(), args[i].length());
    // place the null character in the end
    argv[i][args[i].length()] = '\0';
  }
  argv[argc] = nullptr;
  // class to automatically cleanup argv upon destruction
  class Cleanup_argv
  {
  public:
    Cleanup_argv(char ** argv_, int argc_plus_one_)
      : argv(argv_)
      , argc_plus_one(argc_plus_one_)
    {}

    ~Cleanup_argv()
    {
      for (unsigned int i = 0; i < argc_plus_one; ++i)
      {
        delete[] argv[i];
      }
      delete[] argv;
    }

  private:
    char **      argv;
    unsigned int argc_plus_one;
  };
  Cleanup_argv cleanup_argv(argv, argc + 1);

  // antscout->set_stream( out_stream );

  itk::ants::CommandLineParser::Pointer parser = itk::ants::CommandLineParser::New();

  parser->SetCommand(argv[0]);

  std::string commandDescription = std::string("N4 is a variant of the popular N3 (nonparameteric nonuniform ") +
                                   std::string("normalization) retrospective bias correction algorithm. Based ") +
                                   std::string("on the assumption that the corruption of the low frequency bias ") +
                                   std::string("field can be modeled as a convolution of the intensity histogram ") +
                                   std::string("by a Gaussian, the basic algorithmic protocol is to iterate ") +
                                   std::string("between deconvolving the intensity histogram by a Gaussian, ") +
                                   std::string("remapping the intensities, and then spatially smoothing this ") +
                                   std::string("result by a B-spline modeling of the bias field itself. ") +
                                   std::string("The modifications from and improvements obtained over ") +
                                   std::string("the original N3 algorithm are described in the following paper:  ") +
                                   std::string("N. Tustison et al., N4ITK:  Improved N3 Bias Correction, ") +
                                   std::string("IEEE Transactions on Medical Imaging, 29(6):1310-1320, June 2010.");

  parser->SetCommandDescription(commandDescription);
  N4InitializeCommandLineOptions(parser);

  if (parser->Parse(argc, argv) == EXIT_FAILURE)
  {
    return EXIT_FAILURE;
  }

  if (argc == 1)
  {
    parser->PrintMenu(std::cerr, 5, false);
    return EXIT_FAILURE;
  }
  else if (parser->GetOption("help")->GetFunction() &&
           parser->Convert<bool>(parser->GetOption("help")->GetFunction()->GetName()))
  {
    parser->PrintMenu(std::cout, 5, false);
    return EXIT_SUCCESS;
  }
  else if (parser->GetOption('h')->GetFunction() &&
           parser->Convert<bool>(parser->GetOption('h')->GetFunction()->GetName()))
  {
    parser->PrintMenu(std::cout, 5, true);
    return EXIT_SUCCESS;
  }
  // Show automatic version
  itk::ants::CommandLineParser::OptionType::Pointer versionOption = parser->GetOption("version");
  if (versionOption && versionOption->GetNumberOfFunctions())
  {
    std::string versionFunction = versionOption->GetFunction(0)->GetName();
    ConvertToLowerCase(versionFunction);
    if (versionFunction.compare("1") == 0 || versionFunction.compare("true") == 0)
    {
      // Print Version Information
      std::cout << ANTs::Version::ExtendedVersionString() << std::endl;
      return EXIT_SUCCESS;
    }
  }
  // Get dimensionality
  unsigned int dimension = 3;

  itk::ants::CommandLineParser::OptionType::Pointer dimOption = parser->GetOption("image-dimensionality");
  if (dimOption && dimOption->GetNumberOfFunctions())
  {
    dimension = parser->Convert<unsigned int>(dimOption->GetFunction(0)->GetName());
  }
  else
  {
    // Read in the first intensity image to get the image dimension.
    std::string filename;

    itk::ants::CommandLineParser::OptionType::Pointer imageOption = parser->GetOption("input-image");
    if (imageOption && imageOption->GetNumberOfFunctions() > 0)
    {
      if (imageOption->GetFunction(0)->GetNumberOfParameters() > 0)
      {
        filename = imageOption->GetFunction(0)->GetParameter(0);
      }
      else
      {
        filename = imageOption->GetFunction(0)->GetName();
      }
    }
    else
    {
      std::cerr << "No input images were specified.  Specify an input image"
                << " with the -i option" << std::endl;
      return EXIT_FAILURE;
    }
    itk::ImageIOBase::Pointer imageIO =
      itk::ImageIOFactory::CreateImageIO(filename.c_str(), itk::IOFileModeEnum::ReadMode);
    dimension = imageIO->GetNumberOfDimensions();
  }

  int returnValue = EXIT_FAILURE;

  switch (dimension)
  {
    case 2:
    {
      returnValue = N4<2>(parser);
    }
    break;
    case 3:
    {
      returnValue = N4<3>(parser);
    }
    break;
    case 4:
    {
      returnValue = N4<4>(parser);
    }
    break;
    default:
      std::cout << "Unsupported dimension" << std::endl;
      return EXIT_FAILURE;
  }
  return returnValue;
}
} // namespace ants