<commit_before>#include "antsUtilities.h"
#include "antsAllocImage.h"
#include "itkantsRegistrationHelper.h"
#include "ReadWriteData.h"
#include "TensorFunctions.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkExtractImageFilter.h"
#include "itkResampleImageFilter.h"
#include "itkVectorImage.h"
#include "itkVectorIndexSelectionCastImageFilter.h"

#include "itkAffineTransform.h"
#include "itkCompositeTransform.h"
#include "itkDisplacementFieldTransform.h"
#include "itkIdentityTransform.h"
#include "itkMatrixOffsetTransformBase.h"
#include "itkTransformFactory.h"
#include "itkTransformFileReader.h"
#include "itkTransformToDisplacementFieldFilter.h"

#include "itkBSplineInterpolateImageFunction.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkGaussianInterpolateImageFunction.h"
#include "itkInterpolateImageFunction.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
#include "itkWindowedSincInterpolateImageFunction.h"
#include "itkLabelImageGaussianInterpolateImageFunction.h"
#include "itkLabelImageGenericInterpolateImageFunction.h"

namespace ants
{
template <typename TensorImageType, typename ImageType>
void
CorrectImageTensorDirection( TensorImageType * movingTensorImage, ImageType * referenceImage )
{
  typedef typename TensorImageType::DirectionType    DirectionType;
  typedef typename DirectionType::InternalMatrixType MatrixType;

  // Assume tensors start in moving voxel space, we want to put them in fixed voxel space

  MatrixType direction = referenceImage->GetDirection().GetTranspose() * movingTensorImage->GetDirection().GetVnlMatrix();

  if( !direction.is_identity( 0.00001 ) )
    {
    itk::ImageRegionIterator<TensorImageType> It( movingTensorImage, movingTensorImage->GetBufferedRegion() );
    for( It.GoToBegin(); !It.IsAtEnd(); ++It )
      {
      typedef typename TensorImageType::PixelType                         TensorType;
      typedef typename TensorImageType::DirectionType::InternalMatrixType TensorMatrixType;

      TensorType       tensor = It.Get();
      TensorMatrixType dt;

      Vector2Matrix<TensorType, TensorMatrixType>(tensor, dt);

      dt = direction * dt * direction.transpose();

      tensor = Matrix2Vector<TensorType, TensorMatrixType>(dt);

      It.Set( tensor );
      }
    }
}

template <typename DisplacementFieldType, typename ImageType>
void
CorrectImageVectorDirection( DisplacementFieldType * movingVectorImage, ImageType * referenceImage )
{
  typedef typename DisplacementFieldType::DirectionType DirectionType;

  // Assume vectors are initially described in the voxel space of the moving image, like tensors
  typename DirectionType::InternalMatrixType direction =
    referenceImage->GetDirection().GetTranspose() * movingVectorImage->GetDirection().GetVnlMatrix();

  typedef typename DisplacementFieldType::PixelType VectorType;

  const unsigned int dimension = ImageType::ImageDimension;

  if( !direction.is_identity( 0.00001 ) )
    {
    itk::ImageRegionIterator<DisplacementFieldType> It( movingVectorImage, movingVectorImage->GetBufferedRegion() );
    for( It.GoToBegin(); !It.IsAtEnd(); ++It )
      {
      VectorType vector = It.Get();

      vnl_vector<double> internalVector( dimension );
      for( unsigned int d = 0; d < dimension; d++ )
        {
        internalVector[d] = vector[d];
        }

      internalVector.pre_multiply( direction );;
      for( unsigned int d = 0; d < dimension; d++ )
        {
        vector[d] = internalVector[d];
        }

      It.Set( vector );
      }
    }
}

template <unsigned int NDim>
unsigned int numTensorElements()
{
  return NDim + numTensorElements<NDim - 1>();
}

template <>
unsigned int numTensorElements<0>()
{
  return 0;
}

template <unsigned int NDim>
std::vector<unsigned int> tensorDiagonalArrayIndices()
{
  std::vector<unsigned int> diagElements;
  for( unsigned int d = 0; d < NDim; d++ )
    { // I think this is correct for upper-triangular ordering but only tested on 3D tensors
    diagElements.push_back(d * NDim - d * (d - 1) / 2);
    }
  return diagElements;
}

template <unsigned int NDim>
bool isDiagonalElement(std::vector<unsigned int> diagElements, unsigned int ind)
{
  for( unsigned int i = 0; i < NDim; i++ )
    {
    if ( diagElements[i] == ind )
      {
      return true;
      }
    }

  return false;
}

template <typename T, unsigned int Dimension>
int antsApplyTransforms( itk::ants::CommandLineParser::Pointer & parser, unsigned int inputImageType = 0 )
{
  typedef T                                      RealType;
  typedef T                                      PixelType;
  typedef itk::Vector<RealType, Dimension>       VectorType;

  typedef itk::SymmetricSecondRankTensor<RealType, Dimension> TensorPixelType;

  // typedef unsigned int                     LabelPixelType;
  // typedef itk::Image<PixelType, Dimension> LabelImageType;

  typedef itk::Image<PixelType, Dimension>           ImageType;
  typedef ImageType                                  ReferenceImageType;

  typedef itk::Image<PixelType, Dimension + 1>       TimeSeriesImageType;
  typedef itk::Image<VectorType, Dimension>          DisplacementFieldType;
  typedef itk::VectorImage<PixelType, Dimension>     MultiChannelImageType;
  typedef itk::Image<TensorPixelType, Dimension>     TensorImageType;

  typedef typename ants::RegistrationHelper<T, Dimension>         RegistrationHelperType;
  typedef typename RegistrationHelperType::AffineTransformType    AffineTransformType;
  typedef typename RegistrationHelperType::CompositeTransformType CompositeTransformType;

  const unsigned int NumberOfTensorElements = numTensorElements<Dimension>();

  std::vector<unsigned int> tensorDiagIndices = tensorDiagonalArrayIndices<Dimension>();

  typename TimeSeriesImageType::Pointer timeSeriesImage = nullptr;
  typename TensorImageType::Pointer tensorImage = nullptr;
  typename MultiChannelImageType::Pointer multiChannelImage = nullptr;
  typename DisplacementFieldType::Pointer vectorImage = nullptr;

  std::vector<typename ImageType::Pointer> inputImages;
  inputImages.clear();

  std::vector<typename ImageType::Pointer> outputImages;
  outputImages.clear();

  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() );
    }

  /**
   * Input object option - for now, we're limiting this to images.
   */
  typename itk::ants::CommandLineParser::OptionType::Pointer inputOption = parser->GetOption( "input" );
  typename itk::ants::CommandLineParser::OptionType::Pointer outputOption = parser->GetOption( "output" );

  if( inputImageType == 4 && inputOption && inputOption->GetNumberOfFunctions() )
    {
    if( verbose )
      {
      std::cout << "Input multichannel image: " << inputOption->GetFunction( 0 )->GetName() << std::endl;
      }
    ReadImage<MultiChannelImageType>( multiChannelImage, ( inputOption->GetFunction( 0 )->GetName() ).c_str() );
    timeSeriesImage = ConvertMultiChannelImageToTimeSeriesImage<MultiChannelImageType, TimeSeriesImageType>( multiChannelImage );
    }
  else if( inputImageType == 3 && inputOption && inputOption->GetNumberOfFunctions() )
    {
    if( verbose )
      {
      std::cout << "Input time-series image: " << inputOption->GetFunction( 0 )->GetName() << std::endl;
      }
    ReadImage<TimeSeriesImageType>( timeSeriesImage, ( inputOption->GetFunction( 0 )->GetName() ).c_str() );
    }
  else if( inputImageType == 2 && inputOption && inputOption->GetNumberOfFunctions() )
    {
    if( verbose )
      {
      std::cout << "Input tensor image: " << inputOption->GetFunction( 0 )->GetName() << std::endl;
      }
    ReadTensorImage<TensorImageType>( tensorImage, ( inputOption->GetFunction( 0 )->GetName() ).c_str(), true );
    }
  else if( inputImageType == 0 && inputOption && inputOption->GetNumberOfFunctions() )
    {
    if( verbose )
      {
      std::cout << "Input scalar image: " << inputOption->GetFunction( 0 )->GetName() << std::endl;
      }
    typename ImageType::Pointer image;
    ReadImage<ImageType>( image, ( inputOption->GetFunction( 0 )->GetName() ).c_str()  );
    inputImages.push_back( image );
    }
  else if( inputImageType == 1 && inputOption && inputOption->GetNumberOfFunctions() )
    {
    if( verbose )
      {
      std::cout << "Input vector image: " << inputOption->GetFunction( 0 )->GetName() << std::endl;
      }

    typedef itk::ImageFileReader<DisplacementFieldType> ReaderType;
    typename ReaderType::Pointer reader = ReaderType::New();
    reader->SetFileName( ( inputOption->GetFunction( 0 )->GetName() ).c_str() );

    try
      {
      vectorImage = reader->GetOutput();
      vectorImage->Update();
      vectorImage->DisconnectPipeline();
      }
    catch( ... )
      {
      if( verbose )
        {
        std::cerr << "Unable to read vector image " << reader->GetFileName() << std::endl;
        }
      return EXIT_FAILURE;
      }
    }
  else if( outputOption && outputOption->GetNumberOfFunctions() )
    {
    if( outputOption->GetFunction( 0 )->GetNumberOfParameters() > 1 &&
        parser->Convert<unsigned int>( outputOption->GetFunction( 0 )->GetParameter( 1 ) ) == 0 )
      {
      // Here we have no input image, valid usage is
      // -o [warp.nii.gz, 1] or -o Linear[affine.mat, 0]  
      // Exit failure with -o [outputImage.nii.gz,0] but no input to warp
      std::string outputOptionName = outputOption->GetFunction( 0 )->GetName();
      ConvertToLowerCase( outputOptionName );
      if( std::strcmp( outputOptionName.c_str(), "linear" ) )
        {      
        if( verbose )
          {
          std::cerr << "An input image is required." << std::endl;
          }
        return EXIT_FAILURE;
        }
      }
    }
  /**
   * Reference image option
   */
  bool needReferenceImage = true;
  if( outputOption && outputOption->GetNumberOfFunctions() )
    {
    std::string outputOptionName = outputOption->GetFunction( 0 )->GetName();
    ConvertToLowerCase( outputOptionName );
    if( !std::strcmp( outputOptionName.c_str(), "linear" ) )
      {
      needReferenceImage = false;
      }
    }

  typedef ImageType ReferenceImageType;
  typename ReferenceImageType::Pointer referenceImage;

  typename itk::ants::CommandLineParser::OptionType::Pointer referenceOption =
    parser->GetOption( "reference-image" );

  if( referenceOption && referenceOption->GetNumberOfFunctions() )
    {
    if( verbose )
      {
      std::cout << "Reference image: " << referenceOption->GetFunction( 0 )->GetName() << std::endl;
      }
    ReadImage<ReferenceImageType>( referenceImage,  ( referenceOption->GetFunction( 0 )->GetName() ).c_str() );
    }
  else if( needReferenceImage == true )
    {
    if( verbose )
      {
      std::cerr << "A reference image is required." << std::endl;
      }
    return EXIT_FAILURE;
    }

  if( inputImageType == 1 )
    {
    CorrectImageVectorDirection<DisplacementFieldType, ReferenceImageType>( vectorImage, referenceImage );
    for( unsigned int i = 0; i < Dimension; i++ )
      {
      typedef itk::VectorIndexSelectionCastImageFilter<DisplacementFieldType, ImageType> SelectorType;
      typename SelectorType::Pointer selector = SelectorType::New();
      selector->SetInput( vectorImage );
      selector->SetIndex( i );
      selector->Update();

      inputImages.push_back( selector->GetOutput() );
      }
    }
  else if( inputImageType == 2 )
    {
    CorrectImageTensorDirection<TensorImageType, ReferenceImageType>( tensorImage, referenceImage );
    for( unsigned int i = 0; i < NumberOfTensorElements; i++ )
      {
      typedef itk::VectorIndexSelectionCastImageFilter<TensorImageType, ImageType> SelectorType;
      typename SelectorType::Pointer selector = SelectorType::New();
      selector->SetInput( tensorImage );
      selector->SetIndex( i );
      selector->Update();

      inputImages.push_back( selector->GetOutput() );
      }
    }
  else if( inputImageType == 3 || inputImageType == 4 )
    {
    typename TimeSeriesImageType::RegionType extractRegion = timeSeriesImage->GetLargestPossibleRegion();
    unsigned int numberOfTimePoints = extractRegion.GetSize()[Dimension];
    int          startTimeIndex = extractRegion.GetIndex()[Dimension];

    extractRegion.SetSize( Dimension, 0 );
    for( unsigned int i = 0; i < numberOfTimePoints; i++ )
      {
      extractRegion.SetIndex( Dimension, startTimeIndex + i );

      typedef itk::ExtractImageFilter<TimeSeriesImageType, ImageType> ExtracterType;
      typename ExtracterType::Pointer extracter = ExtracterType::New();
      extracter->SetInput( timeSeriesImage );
      extracter->SetExtractionRegion( extractRegion );
      extracter->SetDirectionCollapseToSubmatrix();
      extracter->Update();

      inputImages.push_back( extracter->GetOutput() );
      }
    }

  /**
   * Transform option
   */
  // Register the matrix offset transform base class to the
  // transform factory for compatibility with the current ANTs.
  typedef itk::MatrixOffsetTransformBase<RealType, Dimension, Dimension> MatrixOffsetTransformType;
  itk::TransformFactory<MatrixOffsetTransformType>::RegisterTransform();

  typename itk::ants::CommandLineParser::OptionType::Pointer transformOption = parser->GetOption( "transform" );

  bool useStaticCastForR = false;
  typename itk::ants::CommandLineParser::OptionType::Pointer rOption =
    parser->GetOption( "static-cast-for-R" );
  if( rOption && rOption->GetNumberOfFunctions() )
    {
    useStaticCastForR = parser->Convert<bool>(  rOption->GetFunction( 0 )->GetName() );
    }

  std::vector<bool> isDerivedTransform;
  typename CompositeTransformType::Pointer compositeTransform =
    GetCompositeTransformFromParserOption<RealType, Dimension>( parser, transformOption, isDerivedTransform,
                                                                useStaticCastForR );
  if( compositeTransform.IsNull() )
    {
    return EXIT_FAILURE;
    }

  if( !compositeTransform->GetNumberOfParameters() )
    {
    if( verbose )
      {
      std::cout << "WARNING: No transforms found, using identify transform" << std::endl;
      }
    typename MatrixOffsetTransformType::Pointer idTransform = MatrixOffsetTransformType::New();
    idTransform->SetIdentity();
    compositeTransform->AddTransform( idTransform );
    }

  std::string whichInterpolator( "linear" );
  typename itk::ants::CommandLineParser::OptionType::Pointer interpolationOption = parser->GetOption( "interpolation" );
  if( interpolationOption && interpolationOption->GetNumberOfFunctions() )
    {
    whichInterpolator = interpolationOption->GetFunction( 0 )->GetName();
    ConvertToLowerCase( whichInterpolator );
    }

  const size_t VImageDimension = Dimension;
  typename ImageType::SpacingType
    cache_spacing_for_smoothing_sigmas(itk::NumericTraits<typename ImageType::SpacingType::ValueType>::ZeroValue());
  if( !std::strcmp( whichInterpolator.c_str(), "gaussian" )
      ||   !std::strcmp( whichInterpolator.c_str(), "multilabel" )
      )
    {
    cache_spacing_for_smoothing_sigmas = inputImages[0]->GetSpacing();
    }

#include "make_interpolator_snip.tmpl"

  /**
   * Default voxel value
   */
  PixelType defaultValue = 0;
  typename itk::ants::CommandLineParser::OptionType::Pointer defaultOption =
    parser->GetOption( "default-value" );
  if( defaultOption && defaultOption->GetNumberOfFunctions() )
    {
    defaultValue = parser->Convert<PixelType>( defaultOption->GetFunction( 0 )->GetName() );
    }
  if( verbose )
    {
    if( inputImageType == 2 )
      std::cout << "Default pixel mean diffusivity: " << defaultValue << std::endl;
    else
      std::cout << "Default pixel value: " << defaultValue << std::endl;
    }
  for( unsigned int n = 0; n < inputImages.size(); n++ )
    {
    typedef itk::ResampleImageFilter<ImageType, ImageType, RealType> ResamplerType;
    typename ResamplerType::Pointer resampleFilter = ResamplerType::New();
    resampleFilter->SetInput( inputImages[n] );
    resampleFilter->SetOutputParametersFromImage( referenceImage );
    resampleFilter->SetTransform( compositeTransform );
    if ( inputImageType == 2 )
      {
      // Set background pixel values in tensor images to produce an isotropic tensor
      if ( defaultValue > 0 && isDiagonalElement<Dimension>(tensorDiagIndices, n) )
        {
        // defaultValue == MD of isotropic tensor. Resampling is done in log space
        resampleFilter->SetDefaultPixelValue( log( defaultValue ) );
        }
             else
        {
        resampleFilter->SetDefaultPixelValue( 0 );
        }
      }
    else
      {
      // for non-tensor images, set the same background value for each component
      resampleFilter->SetDefaultPixelValue( defaultValue );
      }

    interpolator->SetInputImage( inputImages[n] );
    resampleFilter->SetInterpolator( interpolator );
    if( n == 0 )
      {
      if( verbose )
        {
        std::cout << "Interpolation type: " << resampleFilter->GetInterpolator()->GetNameOfClass() << std::endl;
        }
      }
    if( inputImageType == 3 )
      {
      if( verbose )
        {
        std::cout << "  Applying transform(s) to time point " << n << " (out of " << inputImages.size() << ")." << std::endl;
        }
      }
    else if( inputImageType == 4 )
      {
      if( verbose )
        {
        std::cout << "  Applying transform(s) to channel " << n << " (out of " << inputImages.size() << ")." << std::endl;
        }
      }
    resampleFilter->Update();
    outputImages.push_back( resampleFilter->GetOutput() );
    }

  /**
   * output
   */
  if( outputOption && outputOption->GetNumberOfFunctions() )
    {
    std::string outputOptionName = outputOption->GetFunction( 0 )->GetName();
    ConvertToLowerCase( outputOptionName );
    if( !std::strcmp( outputOptionName.c_str(), "linear" ) )
      {
      if( !compositeTransform->IsLinear() )
        {
        if( verbose )
          {
          std::cerr << "The transform or set of transforms is not linear." << std::endl;
          }
        return EXIT_FAILURE;
        }
      else
        {
        typename RegistrationHelperType::Pointer helper = RegistrationHelperType::New();

        typename AffineTransformType::Pointer transform = helper->CollapseLinearTransforms( compositeTransform );

        typedef itk::TransformFileWriterTemplate<T> TransformWriterType;
        typename TransformWriterType::Pointer transformWriter = TransformWriterType::New();
        transformWriter->SetFileName( ( outputOption->GetFunction( 0 )->GetParameter( 0 ) ).c_str() );

        if( outputOption->GetFunction( 0 )->GetNumberOfParameters() > 1 &&
            parser->Convert<unsigned int>( outputOption->GetFunction( 0 )->GetParameter( 1 ) ) != 0 )
          {
          typename AffineTransformType::Pointer inverseTransform = AffineTransformType::New();
          inverseTransform->SetMatrix(
            dynamic_cast<MatrixOffsetTransformType *>( transform->GetInverseTransform().GetPointer() )->GetMatrix() );
          inverseTransform->SetOffset( -( inverseTransform->GetMatrix() * transform->GetOffset() ) );
          transformWriter->SetInput( inverseTransform );
          }
        else
          {
          transformWriter->SetInput( transform );
          }
#if ITK_VERSION_MAJOR >= 5
        transformWriter->SetUseCompression(true);
#endif
        transformWriter->Update();
        }
      }
    else if( outputOption->GetFunction( 0 )->GetNumberOfParameters() > 1 &&
             parser->Convert<unsigned int>( outputOption->GetFunction( 0 )->GetParameter( 1 ) ) != 0 )
      {
      if( verbose )
        {
        std::cout << "Output composite transform displacement field: "
                  << outputOption->GetFunction( 0 )->GetParameter( 0 ) << std::endl;
        }

      typedef typename itk::TransformToDisplacementFieldFilter<DisplacementFieldType, RealType> ConverterType;
      typename ConverterType::Pointer converter = ConverterType::New();
      converter->SetOutputOrigin( referenceImage->GetOrigin() );
      converter->SetOutputStartIndex( referenceImage->GetBufferedRegion().GetIndex() );
      converter->SetSize( referenceImage->GetBufferedRegion().GetSize() );
      converter->SetOutputSpacing( referenceImage->GetSpacing() );
      converter->SetOutputDirection( referenceImage->GetDirection() );
      converter->SetTransform( compositeTransform );
      converter->Update();

      typedef  itk::ImageFileWriter<DisplacementFieldType> DisplacementFieldWriterType;
      typename DisplacementFieldWriterType::Pointer displacementFieldWriter = DisplacementFieldWriterType::New();
      displacementFieldWriter->SetInput( converter->GetOutput() );
      displacementFieldWriter->SetFileName( ( outputOption->GetFunction( 0 )->GetParameter( 0 ) ).c_str() );
      displacementFieldWriter->Update();
      }
    else
      {
      std::string outputFileName = "";
      if( outputOption->GetFunction( 0 )->GetNumberOfParameters() > 1 &&
          parser->Convert<unsigned int>( outputOption->GetFunction( 0 )->GetParameter( 1 ) ) == 0 )
        {
        outputFileName = outputOption->GetFunction( 0 )->GetParameter( 0 );
        }
      else
        {
        outputFileName = outputOption->GetFunction( 0 )->GetName();
        }
      if( verbose )
        {
        std::cout << "Output warped image: " << outputFileName << std::endl;
        }

      if( inputImageType == 1 )
        {
        if( outputImages.size() != Dimension )
          {
          if( verbose )
            {
            std::cerr << "The number of output images does not match the number of vector components." << std::endl;
            }
          return EXIT_FAILURE;
          }

        VectorType zeroVector( 0.0 );

        typename DisplacementFieldType::Pointer outputVectorImage = DisplacementFieldType::New();
        outputVectorImage->CopyInformation( referenceImage );
        outputVectorImage->SetRegions( referenceImage->GetRequestedRegion() );
        outputVectorImage->Allocate();
        outputVectorImage->FillBuffer( zeroVector );

        itk::ImageRegionIteratorWithIndex<DisplacementFieldType> It( outputVectorImage,
                                                                     outputVectorImage->GetRequestedRegion() );
        for( It.GoToBegin(); !It.IsAtEnd(); ++It )
          {
          VectorType vector = It.Get();
          typename DisplacementFieldType::IndexType index = It.GetIndex();
          for( unsigned int n = 0; n < Dimension; n++ )
            {
            vector.SetNthComponent( n, outputImages[n]->GetPixel( index ) );
            }
          It.Set( vector );
          }
        typedef  itk::ImageFileWriter<DisplacementFieldType> WriterType;
        typename WriterType::Pointer writer = WriterType::New();
        writer->SetInput( outputVectorImage );
        writer->SetFileName( ( outputFileName ).c_str() );
        writer->Update();
        }
      else if( inputImageType == 2 )
        {
        if( outputImages.size() != NumberOfTensorElements )
          {
          if( verbose )
            {
            std::cerr << "The number of output images does not match the number of tensor elements." << std::endl;
            }
          return EXIT_FAILURE;
          }

        TensorPixelType zeroTensor( 0.0 );

        typename TensorImageType::Pointer outputTensorImage = TensorImageType::New();
        outputTensorImage->CopyInformation( referenceImage );
        outputTensorImage->SetRegions( referenceImage->GetRequestedRegion() );
        outputTensorImage->Allocate();
        outputTensorImage->FillBuffer( zeroTensor );

        itk::ImageRegionIteratorWithIndex<TensorImageType> It( outputTensorImage,
                                                               outputTensorImage->GetRequestedRegion() );
        for( It.GoToBegin(); !It.IsAtEnd(); ++It )
          {
          TensorPixelType tensor = It.Get();
          typename TensorImageType::IndexType index = It.GetIndex();
          for( unsigned int n = 0; n < NumberOfTensorElements; n++ )
            {
            tensor.SetNthComponent( n, outputImages[n]->GetPixel( index ) );
            }
          It.Set( tensor );
          }
        WriteTensorImage<TensorImageType>( outputTensorImage, ( outputFileName ).c_str(), true );
        }
      else if( inputImageType == 3 || inputImageType == 4 )
        {
        unsigned int numberOfTimePoints = timeSeriesImage->GetLargestPossibleRegion().GetSize()[Dimension];

        if( outputImages.size() != numberOfTimePoints )
          {
          if( verbose )
            {
            std::cerr << "The number of output images does not match the number of image time/channel points." << std::endl;
            }
          return EXIT_FAILURE;
          }

        typename TimeSeriesImageType::Pointer outputTimeSeriesImage = TimeSeriesImageType::New();

        typename TimeSeriesImageType::PointType origin = timeSeriesImage->GetOrigin();
        typename TimeSeriesImageType::SizeType size = timeSeriesImage->GetLargestPossibleRegion().GetSize();
        typename TimeSeriesImageType::DirectionType direction = timeSeriesImage->GetDirection();
        typename TimeSeriesImageType::IndexType index = timeSeriesImage->GetLargestPossibleRegion().GetIndex();
        typename TimeSeriesImageType::SpacingType spacing = timeSeriesImage->GetSpacing();

        for( unsigned int i = 0; i < Dimension; i++ )
          {
          origin[i] = referenceImage->GetOrigin()[i];
          size[i] = referenceImage->GetRequestedRegion().GetSize()[i];
          index[i] = referenceImage->GetRequestedRegion().GetIndex()[i];
          spacing[i] = referenceImage->GetSpacing()[i];
          for( unsigned int j = 0; j < Dimension; j++ )
            {
            direction[i][j] = referenceImage->GetDirection()[i][j];
            }
          }

        typename TimeSeriesImageType::RegionType region;
        region.SetSize( size );
        region.SetIndex( index );

        int startTimeIndex = timeSeriesImage->GetLargestPossibleRegion().GetIndex()[Dimension];

        outputTimeSeriesImage->CopyInformation( timeSeriesImage );
        outputTimeSeriesImage->SetOrigin( origin );
        outputTimeSeriesImage->SetDirection( direction );
        outputTimeSeriesImage->SetSpacing( spacing );
        outputTimeSeriesImage->SetRegions( region );
        outputTimeSeriesImage->Allocate();
        outputTimeSeriesImage->FillBuffer( 0 );

        typename ImageType::IndexType referenceIndex;

        itk::ImageRegionIteratorWithIndex<TimeSeriesImageType> It( outputTimeSeriesImage,
                                                                   outputTimeSeriesImage->GetRequestedRegion() );
        for( It.GoToBegin(); !It.IsAtEnd(); ++It )
          {
          typename TimeSeriesImageType::IndexType timeImageIndex = It.GetIndex();
          for( unsigned int i = 0; i < Dimension; i++ )
            {
            referenceIndex[i] = timeImageIndex[i];
            }
          It.Set( outputImages[timeImageIndex[Dimension] - startTimeIndex]->GetPixel( referenceIndex ) );
          }
        if( inputImageType == 3 )  
          {
          WriteImage<TimeSeriesImageType>( outputTimeSeriesImage, ( outputFileName ).c_str() );
          }
        else // inputImageType == 4
          {
          multiChannelImage = ConvertTimeSeriesImageToMultiChannelImage<TimeSeriesImageType, MultiChannelImageType>( timeSeriesImage );<commit_msg>BUG: antsApplyTransforms -e 4 would write input image rather than deformed image<commit_after>#include "antsUtilities.h"
#include "antsAllocImage.h"
#include "itkantsRegistrationHelper.h"
#include "ReadWriteData.h"
#include "TensorFunctions.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkExtractImageFilter.h"
#include "itkResampleImageFilter.h"
#include "itkVectorImage.h"
#include "itkVectorIndexSelectionCastImageFilter.h"

#include "itkAffineTransform.h"
#include "itkCompositeTransform.h"
#include "itkDisplacementFieldTransform.h"
#include "itkIdentityTransform.h"
#include "itkMatrixOffsetTransformBase.h"
#include "itkTransformFactory.h"
#include "itkTransformFileReader.h"
#include "itkTransformToDisplacementFieldFilter.h"

#include "itkBSplineInterpolateImageFunction.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkGaussianInterpolateImageFunction.h"
#include "itkInterpolateImageFunction.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
#include "itkWindowedSincInterpolateImageFunction.h"
#include "itkLabelImageGaussianInterpolateImageFunction.h"
#include "itkLabelImageGenericInterpolateImageFunction.h"

namespace ants
{
template <typename TensorImageType, typename ImageType>
void
CorrectImageTensorDirection( TensorImageType * movingTensorImage, ImageType * referenceImage )
{
  typedef typename TensorImageType::DirectionType    DirectionType;
  typedef typename DirectionType::InternalMatrixType MatrixType;

  // Assume tensors start in moving voxel space, we want to put them in fixed voxel space

  MatrixType direction = referenceImage->GetDirection().GetTranspose() * movingTensorImage->GetDirection().GetVnlMatrix();

  if( !direction.is_identity( 0.00001 ) )
    {
    itk::ImageRegionIterator<TensorImageType> It( movingTensorImage, movingTensorImage->GetBufferedRegion() );
    for( It.GoToBegin(); !It.IsAtEnd(); ++It )
      {
      typedef typename TensorImageType::PixelType                         TensorType;
      typedef typename TensorImageType::DirectionType::InternalMatrixType TensorMatrixType;

      TensorType       tensor = It.Get();
      TensorMatrixType dt;

      Vector2Matrix<TensorType, TensorMatrixType>(tensor, dt);

      dt = direction * dt * direction.transpose();

      tensor = Matrix2Vector<TensorType, TensorMatrixType>(dt);

      It.Set( tensor );
      }
    }
}

template <typename DisplacementFieldType, typename ImageType>
void
CorrectImageVectorDirection( DisplacementFieldType * movingVectorImage, ImageType * referenceImage )
{
  typedef typename DisplacementFieldType::DirectionType DirectionType;

  // Assume vectors are initially described in the voxel space of the moving image, like tensors
  typename DirectionType::InternalMatrixType direction =
    referenceImage->GetDirection().GetTranspose() * movingVectorImage->GetDirection().GetVnlMatrix();

  typedef typename DisplacementFieldType::PixelType VectorType;

  const unsigned int dimension = ImageType::ImageDimension;

  if( !direction.is_identity( 0.00001 ) )
    {
    itk::ImageRegionIterator<DisplacementFieldType> It( movingVectorImage, movingVectorImage->GetBufferedRegion() );
    for( It.GoToBegin(); !It.IsAtEnd(); ++It )
      {
      VectorType vector = It.Get();

      vnl_vector<double> internalVector( dimension );
      for( unsigned int d = 0; d < dimension; d++ )
        {
        internalVector[d] = vector[d];
        }

      internalVector.pre_multiply( direction );;
      for( unsigned int d = 0; d < dimension; d++ )
        {
        vector[d] = internalVector[d];
        }

      It.Set( vector );
      }
    }
}

template <unsigned int NDim>
unsigned int numTensorElements()
{
  return NDim + numTensorElements<NDim - 1>();
}

template <>
unsigned int numTensorElements<0>()
{
  return 0;
}

template <unsigned int NDim>
std::vector<unsigned int> tensorDiagonalArrayIndices()
{
  std::vector<unsigned int> diagElements;
  for( unsigned int d = 0; d < NDim; d++ )
    { // I think this is correct for upper-triangular ordering but only tested on 3D tensors
    diagElements.push_back(d * NDim - d * (d - 1) / 2);
    }
  return diagElements;
}

template <unsigned int NDim>
bool isDiagonalElement(std::vector<unsigned int> diagElements, unsigned int ind)
{
  for( unsigned int i = 0; i < NDim; i++ )
    {
    if ( diagElements[i] == ind )
      {
      return true;
      }
    }

  return false;
}

template <typename T, unsigned int Dimension>
int antsApplyTransforms( itk::ants::CommandLineParser::Pointer & parser, unsigned int inputImageType = 0 )
{
  typedef T                                      RealType;
  typedef T                                      PixelType;
  typedef itk::Vector<RealType, Dimension>       VectorType;

  typedef itk::SymmetricSecondRankTensor<RealType, Dimension> TensorPixelType;

  // typedef unsigned int                     LabelPixelType;
  // typedef itk::Image<PixelType, Dimension> LabelImageType;

  typedef itk::Image<PixelType, Dimension>           ImageType;
  typedef ImageType                                  ReferenceImageType;

  typedef itk::Image<PixelType, Dimension + 1>       TimeSeriesImageType;
  typedef itk::Image<VectorType, Dimension>          DisplacementFieldType;
  typedef itk::VectorImage<PixelType, Dimension>     MultiChannelImageType;
  typedef itk::Image<TensorPixelType, Dimension>     TensorImageType;

  typedef typename ants::RegistrationHelper<T, Dimension>         RegistrationHelperType;
  typedef typename RegistrationHelperType::AffineTransformType    AffineTransformType;
  typedef typename RegistrationHelperType::CompositeTransformType CompositeTransformType;

  const unsigned int NumberOfTensorElements = numTensorElements<Dimension>();

  std::vector<unsigned int> tensorDiagIndices = tensorDiagonalArrayIndices<Dimension>();

  typename TimeSeriesImageType::Pointer timeSeriesImage = nullptr;
  typename TensorImageType::Pointer tensorImage = nullptr;
  typename MultiChannelImageType::Pointer multiChannelImage = nullptr;
  typename DisplacementFieldType::Pointer vectorImage = nullptr;

  std::vector<typename ImageType::Pointer> inputImages;
  inputImages.clear();

  std::vector<typename ImageType::Pointer> outputImages;
  outputImages.clear();

  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() );
    }

  /**
   * Input object option - for now, we're limiting this to images.
   */
  typename itk::ants::CommandLineParser::OptionType::Pointer inputOption = parser->GetOption( "input" );
  typename itk::ants::CommandLineParser::OptionType::Pointer outputOption = parser->GetOption( "output" );

  if( inputImageType == 4 && inputOption && inputOption->GetNumberOfFunctions() )
    {
    if( verbose )
      {
      std::cout << "Input multichannel image: " << inputOption->GetFunction( 0 )->GetName() << std::endl;
      }
    ReadImage<MultiChannelImageType>( multiChannelImage, ( inputOption->GetFunction( 0 )->GetName() ).c_str() );
    timeSeriesImage = ConvertMultiChannelImageToTimeSeriesImage<MultiChannelImageType, TimeSeriesImageType>( multiChannelImage );
    }
  else if( inputImageType == 3 && inputOption && inputOption->GetNumberOfFunctions() )
    {
    if( verbose )
      {
      std::cout << "Input time-series image: " << inputOption->GetFunction( 0 )->GetName() << std::endl;
      }
    ReadImage<TimeSeriesImageType>( timeSeriesImage, ( inputOption->GetFunction( 0 )->GetName() ).c_str() );
    }
  else if( inputImageType == 2 && inputOption && inputOption->GetNumberOfFunctions() )
    {
    if( verbose )
      {
      std::cout << "Input tensor image: " << inputOption->GetFunction( 0 )->GetName() << std::endl;
      }
    ReadTensorImage<TensorImageType>( tensorImage, ( inputOption->GetFunction( 0 )->GetName() ).c_str(), true );
    }
  else if( inputImageType == 0 && inputOption && inputOption->GetNumberOfFunctions() )
    {
    if( verbose )
      {
      std::cout << "Input scalar image: " << inputOption->GetFunction( 0 )->GetName() << std::endl;
      }
    typename ImageType::Pointer image;
    ReadImage<ImageType>( image, ( inputOption->GetFunction( 0 )->GetName() ).c_str()  );
    inputImages.push_back( image );
    }
  else if( inputImageType == 1 && inputOption && inputOption->GetNumberOfFunctions() )
    {
    if( verbose )
      {
      std::cout << "Input vector image: " << inputOption->GetFunction( 0 )->GetName() << std::endl;
      }

    typedef itk::ImageFileReader<DisplacementFieldType> ReaderType;
    typename ReaderType::Pointer reader = ReaderType::New();
    reader->SetFileName( ( inputOption->GetFunction( 0 )->GetName() ).c_str() );

    try
      {
      vectorImage = reader->GetOutput();
      vectorImage->Update();
      vectorImage->DisconnectPipeline();
      }
    catch( ... )
      {
      if( verbose )
        {
        std::cerr << "Unable to read vector image " << reader->GetFileName() << std::endl;
        }
      return EXIT_FAILURE;
      }
    }
  else if( outputOption && outputOption->GetNumberOfFunctions() )
    {
    if( outputOption->GetFunction( 0 )->GetNumberOfParameters() > 1 &&
        parser->Convert<unsigned int>( outputOption->GetFunction( 0 )->GetParameter( 1 ) ) == 0 )
      {
      // Here we have no input image, valid usage is
      // -o [warp.nii.gz, 1] or -o Linear[affine.mat, 0]  
      // Exit failure with -o [outputImage.nii.gz,0] but no input to warp
      std::string outputOptionName = outputOption->GetFunction( 0 )->GetName();
      ConvertToLowerCase( outputOptionName );
      if( std::strcmp( outputOptionName.c_str(), "linear" ) )
        {      
        if( verbose )
          {
          std::cerr << "An input image is required." << std::endl;
          }
        return EXIT_FAILURE;
        }
      }
    }
  /**
   * Reference image option
   */
  bool needReferenceImage = true;
  if( outputOption && outputOption->GetNumberOfFunctions() )
    {
    std::string outputOptionName = outputOption->GetFunction( 0 )->GetName();
    ConvertToLowerCase( outputOptionName );
    if( !std::strcmp( outputOptionName.c_str(), "linear" ) )
      {
      needReferenceImage = false;
      }
    }

  typedef ImageType ReferenceImageType;
  typename ReferenceImageType::Pointer referenceImage;

  typename itk::ants::CommandLineParser::OptionType::Pointer referenceOption =
    parser->GetOption( "reference-image" );

  if( referenceOption && referenceOption->GetNumberOfFunctions() )
    {
    if( verbose )
      {
      std::cout << "Reference image: " << referenceOption->GetFunction( 0 )->GetName() << std::endl;
      }
    ReadImage<ReferenceImageType>( referenceImage,  ( referenceOption->GetFunction( 0 )->GetName() ).c_str() );
    }
  else if( needReferenceImage == true )
    {
    if( verbose )
      {
      std::cerr << "A reference image is required." << std::endl;
      }
    return EXIT_FAILURE;
    }

  if( inputImageType == 1 )
    {
    CorrectImageVectorDirection<DisplacementFieldType, ReferenceImageType>( vectorImage, referenceImage );
    for( unsigned int i = 0; i < Dimension; i++ )
      {
      typedef itk::VectorIndexSelectionCastImageFilter<DisplacementFieldType, ImageType> SelectorType;
      typename SelectorType::Pointer selector = SelectorType::New();
      selector->SetInput( vectorImage );
      selector->SetIndex( i );
      selector->Update();

      inputImages.push_back( selector->GetOutput() );
      }
    }
  else if( inputImageType == 2 )
    {
    CorrectImageTensorDirection<TensorImageType, ReferenceImageType>( tensorImage, referenceImage );
    for( unsigned int i = 0; i < NumberOfTensorElements; i++ )
      {
      typedef itk::VectorIndexSelectionCastImageFilter<TensorImageType, ImageType> SelectorType;
      typename SelectorType::Pointer selector = SelectorType::New();
      selector->SetInput( tensorImage );
      selector->SetIndex( i );
      selector->Update();

      inputImages.push_back( selector->GetOutput() );
      }
    }
  else if( inputImageType == 3 || inputImageType == 4 )
    {
    typename TimeSeriesImageType::RegionType extractRegion = timeSeriesImage->GetLargestPossibleRegion();
    unsigned int numberOfTimePoints = extractRegion.GetSize()[Dimension];
    int          startTimeIndex = extractRegion.GetIndex()[Dimension];

    extractRegion.SetSize( Dimension, 0 );
    for( unsigned int i = 0; i < numberOfTimePoints; i++ )
      {
      extractRegion.SetIndex( Dimension, startTimeIndex + i );

      typedef itk::ExtractImageFilter<TimeSeriesImageType, ImageType> ExtracterType;
      typename ExtracterType::Pointer extracter = ExtracterType::New();
      extracter->SetInput( timeSeriesImage );
      extracter->SetExtractionRegion( extractRegion );
      extracter->SetDirectionCollapseToSubmatrix();
      extracter->Update();

      inputImages.push_back( extracter->GetOutput() );
      }
    }

  /**
   * Transform option
   */
  // Register the matrix offset transform base class to the
  // transform factory for compatibility with the current ANTs.
  typedef itk::MatrixOffsetTransformBase<RealType, Dimension, Dimension> MatrixOffsetTransformType;
  itk::TransformFactory<MatrixOffsetTransformType>::RegisterTransform();

  typename itk::ants::CommandLineParser::OptionType::Pointer transformOption = parser->GetOption( "transform" );

  bool useStaticCastForR = false;
  typename itk::ants::CommandLineParser::OptionType::Pointer rOption =
    parser->GetOption( "static-cast-for-R" );
  if( rOption && rOption->GetNumberOfFunctions() )
    {
    useStaticCastForR = parser->Convert<bool>(  rOption->GetFunction( 0 )->GetName() );
    }

  std::vector<bool> isDerivedTransform;
  typename CompositeTransformType::Pointer compositeTransform =
    GetCompositeTransformFromParserOption<RealType, Dimension>( parser, transformOption, isDerivedTransform,
                                                                useStaticCastForR );
  if( compositeTransform.IsNull() )
    {
    return EXIT_FAILURE;
    }

  if( !compositeTransform->GetNumberOfParameters() )
    {
    if( verbose )
      {
      std::cout << "WARNING: No transforms found, using identify transform" << std::endl;
      }
    typename MatrixOffsetTransformType::Pointer idTransform = MatrixOffsetTransformType::New();
    idTransform->SetIdentity();
    compositeTransform->AddTransform( idTransform );
    }

  std::string whichInterpolator( "linear" );
  typename itk::ants::CommandLineParser::OptionType::Pointer interpolationOption = parser->GetOption( "interpolation" );
  if( interpolationOption && interpolationOption->GetNumberOfFunctions() )
    {
    whichInterpolator = interpolationOption->GetFunction( 0 )->GetName();
    ConvertToLowerCase( whichInterpolator );
    }

  const size_t VImageDimension = Dimension;
  typename ImageType::SpacingType
    cache_spacing_for_smoothing_sigmas(itk::NumericTraits<typename ImageType::SpacingType::ValueType>::ZeroValue());
  if( !std::strcmp( whichInterpolator.c_str(), "gaussian" )
      ||   !std::strcmp( whichInterpolator.c_str(), "multilabel" )
      )
    {
    cache_spacing_for_smoothing_sigmas = inputImages[0]->GetSpacing();
    }

#include "make_interpolator_snip.tmpl"

  /**
   * Default voxel value
   */
  PixelType defaultValue = 0;
  typename itk::ants::CommandLineParser::OptionType::Pointer defaultOption =
    parser->GetOption( "default-value" );
  if( defaultOption && defaultOption->GetNumberOfFunctions() )
    {
    defaultValue = parser->Convert<PixelType>( defaultOption->GetFunction( 0 )->GetName() );
    }
  if( verbose )
    {
    if( inputImageType == 2 )
      std::cout << "Default pixel mean diffusivity: " << defaultValue << std::endl;
    else
      std::cout << "Default pixel value: " << defaultValue << std::endl;
    }
  for( unsigned int n = 0; n < inputImages.size(); n++ )
    {
    typedef itk::ResampleImageFilter<ImageType, ImageType, RealType> ResamplerType;
    typename ResamplerType::Pointer resampleFilter = ResamplerType::New();
    resampleFilter->SetInput( inputImages[n] );
    resampleFilter->SetOutputParametersFromImage( referenceImage );
    resampleFilter->SetTransform( compositeTransform );
    if ( inputImageType == 2 )
      {
      // Set background pixel values in tensor images to produce an isotropic tensor
      if ( defaultValue > 0 && isDiagonalElement<Dimension>(tensorDiagIndices, n) )
        {
        // defaultValue == MD of isotropic tensor. Resampling is done in log space
        resampleFilter->SetDefaultPixelValue( log( defaultValue ) );
        }
             else
        {
        resampleFilter->SetDefaultPixelValue( 0 );
        }
      }
    else
      {
      // for non-tensor images, set the same background value for each component
      resampleFilter->SetDefaultPixelValue( defaultValue );
      }

    interpolator->SetInputImage( inputImages[n] );
    resampleFilter->SetInterpolator( interpolator );
    if( n == 0 )
      {
      if( verbose )
        {
        std::cout << "Interpolation type: " << resampleFilter->GetInterpolator()->GetNameOfClass() << std::endl;
        }
      }
    if( inputImageType == 3 )
      {
      if( verbose )
        {
        std::cout << "  Applying transform(s) to time point " << n << " (out of " << inputImages.size() << ")." << std::endl;
        }
      }
    else if( inputImageType == 4 )
      {
      if( verbose )
        {
        std::cout << "  Applying transform(s) to channel " << n << " (out of " << inputImages.size() << ")." << std::endl;
        }
      }
    resampleFilter->Update();
    outputImages.push_back( resampleFilter->GetOutput() );
    }

  /**
   * output
   */
  if( outputOption && outputOption->GetNumberOfFunctions() )
    {
    std::string outputOptionName = outputOption->GetFunction( 0 )->GetName();
    ConvertToLowerCase( outputOptionName );
    if( !std::strcmp( outputOptionName.c_str(), "linear" ) )
      {
      if( !compositeTransform->IsLinear() )
        {
        if( verbose )
          {
          std::cerr << "The transform or set of transforms is not linear." << std::endl;
          }
        return EXIT_FAILURE;
        }
      else
        {
        typename RegistrationHelperType::Pointer helper = RegistrationHelperType::New();

        typename AffineTransformType::Pointer transform = helper->CollapseLinearTransforms( compositeTransform );

        typedef itk::TransformFileWriterTemplate<T> TransformWriterType;
        typename TransformWriterType::Pointer transformWriter = TransformWriterType::New();
        transformWriter->SetFileName( ( outputOption->GetFunction( 0 )->GetParameter( 0 ) ).c_str() );

        if( outputOption->GetFunction( 0 )->GetNumberOfParameters() > 1 &&
            parser->Convert<unsigned int>( outputOption->GetFunction( 0 )->GetParameter( 1 ) ) != 0 )
          {
          typename AffineTransformType::Pointer inverseTransform = AffineTransformType::New();
          inverseTransform->SetMatrix(
            dynamic_cast<MatrixOffsetTransformType *>( transform->GetInverseTransform().GetPointer() )->GetMatrix() );
          inverseTransform->SetOffset( -( inverseTransform->GetMatrix() * transform->GetOffset() ) );
          transformWriter->SetInput( inverseTransform );
          }
        else
          {
          transformWriter->SetInput( transform );
          }
#if ITK_VERSION_MAJOR >= 5
        transformWriter->SetUseCompression(true);
#endif
        transformWriter->Update();
        }
      }
    else if( outputOption->GetFunction( 0 )->GetNumberOfParameters() > 1 &&
             parser->Convert<unsigned int>( outputOption->GetFunction( 0 )->GetParameter( 1 ) ) != 0 )
      {
      if( verbose )
        {
        std::cout << "Output composite transform displacement field: "
                  << outputOption->GetFunction( 0 )->GetParameter( 0 ) << std::endl;
        }

      typedef typename itk::TransformToDisplacementFieldFilter<DisplacementFieldType, RealType> ConverterType;
      typename ConverterType::Pointer converter = ConverterType::New();
      converter->SetOutputOrigin( referenceImage->GetOrigin() );
      converter->SetOutputStartIndex( referenceImage->GetBufferedRegion().GetIndex() );
      converter->SetSize( referenceImage->GetBufferedRegion().GetSize() );
      converter->SetOutputSpacing( referenceImage->GetSpacing() );
      converter->SetOutputDirection( referenceImage->GetDirection() );
      converter->SetTransform( compositeTransform );
      converter->Update();

      typedef  itk::ImageFileWriter<DisplacementFieldType> DisplacementFieldWriterType;
      typename DisplacementFieldWriterType::Pointer displacementFieldWriter = DisplacementFieldWriterType::New();
      displacementFieldWriter->SetInput( converter->GetOutput() );
      displacementFieldWriter->SetFileName( ( outputOption->GetFunction( 0 )->GetParameter( 0 ) ).c_str() );
      displacementFieldWriter->Update();
      }
    else
      {
      std::string outputFileName = "";
      if( outputOption->GetFunction( 0 )->GetNumberOfParameters() > 1 &&
          parser->Convert<unsigned int>( outputOption->GetFunction( 0 )->GetParameter( 1 ) ) == 0 )
        {
        outputFileName = outputOption->GetFunction( 0 )->GetParameter( 0 );
        }
      else
        {
        outputFileName = outputOption->GetFunction( 0 )->GetName();
        }
      if( verbose )
        {
        std::cout << "Output warped image: " << outputFileName << std::endl;
        }

      if( inputImageType == 1 )
        {
        if( outputImages.size() != Dimension )
          {
          if( verbose )
            {
            std::cerr << "The number of output images does not match the number of vector components." << std::endl;
            }
          return EXIT_FAILURE;
          }

        VectorType zeroVector( 0.0 );

        typename DisplacementFieldType::Pointer outputVectorImage = DisplacementFieldType::New();
        outputVectorImage->CopyInformation( referenceImage );
        outputVectorImage->SetRegions( referenceImage->GetRequestedRegion() );
        outputVectorImage->Allocate();
        outputVectorImage->FillBuffer( zeroVector );

        itk::ImageRegionIteratorWithIndex<DisplacementFieldType> It( outputVectorImage,
                                                                     outputVectorImage->GetRequestedRegion() );
        for( It.GoToBegin(); !It.IsAtEnd(); ++It )
          {
          VectorType vector = It.Get();
          typename DisplacementFieldType::IndexType index = It.GetIndex();
          for( unsigned int n = 0; n < Dimension; n++ )
            {
            vector.SetNthComponent( n, outputImages[n]->GetPixel( index ) );
            }
          It.Set( vector );
          }
        typedef  itk::ImageFileWriter<DisplacementFieldType> WriterType;
        typename WriterType::Pointer writer = WriterType::New();
        writer->SetInput( outputVectorImage );
        writer->SetFileName( ( outputFileName ).c_str() );
        writer->Update();
        }
      else if( inputImageType == 2 )
        {
        if( outputImages.size() != NumberOfTensorElements )
          {
          if( verbose )
            {
            std::cerr << "The number of output images does not match the number of tensor elements." << std::endl;
            }
          return EXIT_FAILURE;
          }

        TensorPixelType zeroTensor( 0.0 );

        typename TensorImageType::Pointer outputTensorImage = TensorImageType::New();
        outputTensorImage->CopyInformation( referenceImage );
        outputTensorImage->SetRegions( referenceImage->GetRequestedRegion() );
        outputTensorImage->Allocate();
        outputTensorImage->FillBuffer( zeroTensor );

        itk::ImageRegionIteratorWithIndex<TensorImageType> It( outputTensorImage,
                                                               outputTensorImage->GetRequestedRegion() );
        for( It.GoToBegin(); !It.IsAtEnd(); ++It )
          {
          TensorPixelType tensor = It.Get();
          typename TensorImageType::IndexType index = It.GetIndex();
          for( unsigned int n = 0; n < NumberOfTensorElements; n++ )
            {
            tensor.SetNthComponent( n, outputImages[n]->GetPixel( index ) );
            }
          It.Set( tensor );
          }
        WriteTensorImage<TensorImageType>( outputTensorImage, ( outputFileName ).c_str(), true );
        }
      else if( inputImageType == 3 || inputImageType == 4 )
        {
        unsigned int numberOfTimePoints = timeSeriesImage->GetLargestPossibleRegion().GetSize()[Dimension];

        if( outputImages.size() != numberOfTimePoints )
          {
          if( verbose )
            {
            std::cerr << "The number of output images does not match the number of image time/channel points." << std::endl;
            }
          return EXIT_FAILURE;
          }

        typename TimeSeriesImageType::Pointer outputTimeSeriesImage = TimeSeriesImageType::New();

        typename TimeSeriesImageType::PointType origin = timeSeriesImage->GetOrigin();
        typename TimeSeriesImageType::SizeType size = timeSeriesImage->GetLargestPossibleRegion().GetSize();
        typename TimeSeriesImageType::DirectionType direction = timeSeriesImage->GetDirection();
        typename TimeSeriesImageType::IndexType index = timeSeriesImage->GetLargestPossibleRegion().GetIndex();
        typename TimeSeriesImageType::SpacingType spacing = timeSeriesImage->GetSpacing();

        for( unsigned int i = 0; i < Dimension; i++ )
          {
          origin[i] = referenceImage->GetOrigin()[i];
          size[i] = referenceImage->GetRequestedRegion().GetSize()[i];
          index[i] = referenceImage->GetRequestedRegion().GetIndex()[i];
          spacing[i] = referenceImage->GetSpacing()[i];
          for( unsigned int j = 0; j < Dimension; j++ )
            {
            direction[i][j] = referenceImage->GetDirection()[i][j];
            }
          }

        typename TimeSeriesImageType::RegionType region;
        region.SetSize( size );
        region.SetIndex( index );

        int startTimeIndex = timeSeriesImage->GetLargestPossibleRegion().GetIndex()[Dimension];

        outputTimeSeriesImage->CopyInformation( timeSeriesImage );
        outputTimeSeriesImage->SetOrigin( origin );
        outputTimeSeriesImage->SetDirection( direction );
        outputTimeSeriesImage->SetSpacing( spacing );
        outputTimeSeriesImage->SetRegions( region );
        outputTimeSeriesImage->Allocate();
        outputTimeSeriesImage->FillBuffer( 0 );

        typename ImageType::IndexType referenceIndex;

        itk::ImageRegionIteratorWithIndex<TimeSeriesImageType> It( outputTimeSeriesImage,
                                                                   outputTimeSeriesImage->GetRequestedRegion() );
        for( It.GoToBegin(); !It.IsAtEnd(); ++It )
          {
          typename TimeSeriesImageType::IndexType timeImageIndex = It.GetIndex();
          for( unsigned int i = 0; i < Dimension; i++ )
            {
            referenceIndex[i] = timeImageIndex[i];
            }
          It.Set( outputImages[timeImageIndex[Dimension] - startTimeIndex]->GetPixel( referenceIndex ) );
          }
        if( inputImageType == 3 )  
          {
          WriteImage<TimeSeriesImageType>( outputTimeSeriesImage, ( outputFileName ).c_str() );
          }
        else // inputImageType == 4
          {
          multiChannelImage = ConvertTimeSeriesImageToMultiChannelImage<TimeSeriesImageType, MultiChannelImageType>( outputTimeSeriesImage );
{'fraction_non_alphanumeric': 0.11630216892709443, 'fraction_numerical': 0.004243945605258881, 'mean_word_length': 3.166904958988637, 'pattern_counts': {'":': 0, '<': 413, '<?xml version=': 0, '>': 551, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5502978', 'n_tokens_mistral': 15431, 'n_tokens_neox': 14803, 'n_words': 3193}