/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         https://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/

#include "itkArray.h"
#include "itkVectorImage.h"
#include "itkMINCImageIO.h"
#include "itkMINCImageIOFactory.h"
#include "itkIOTestHelper.h"
#include "itkMetaDataObject.h"
#include "itkObjectFactoryBase.h"

static void
RandomPix(vnl_random &                   randgen,
          itk::RGBPixel<unsigned char> & pix,
          double                         _max = itk::NumericTraits<unsigned char>::max())
{
  for (unsigned int i = 0; i < 3; ++i)
  {
    pix[i] = randgen.lrand32(_max);
  }
}

static void
RandomPix(vnl_random & randgen, itk::RGBPixel<char> & pix, double _max = itk::NumericTraits<char>::max())
{
  for (unsigned int i = 0; i < 3; ++i)
  {
    pix[i] = randgen.lrand32(_max);
  }
}


static void
RandomPix(vnl_random &                     randgen,
          itk::Vector<unsigned short, 3> & pix,
          double                           _max = itk::NumericTraits<unsigned short>::max())
{
  for (unsigned int i = 0; i < 3; ++i)
  {
    pix[i] = randgen.lrand32(_max);
  }
}

static void
RandomPix(vnl_random & randgen, itk::Vector<short, 3> & pix, double _max = itk::NumericTraits<short>::max())
{
  for (unsigned int i = 0; i < 3; ++i)
  {
    pix[i] = randgen.lrand32(_max);
  }
}

static void
RandomPix(vnl_random &                   randgen,
          itk::Vector<unsigned int, 3> & pix,
          double                         _max = itk::NumericTraits<unsigned int>::max())
{
  (void)_max;
  for (unsigned int i = 0; i < 3; ++i)
  {
    pix[i] = randgen.lrand32();
  }
}

static void
RandomPix(vnl_random & randgen, itk::Vector<int, 3> & pix, double _max = itk::NumericTraits<int>::max())
{
  (void)_max;
  for (unsigned int i = 0; i < 3; ++i)
  {
    pix[i] = randgen.lrand32();
  }
}

static void
RandomPix(vnl_random & randgen, itk::Vector<float, 3> & pix, float _max = itk::NumericTraits<float>::max())
{
  for (unsigned int i = 0; i < 3; ++i)
  {
    pix[i] = randgen.drand64(_max);
  }
}

static void
RandomPix(vnl_random & randgen, itk::Vector<double, 3> & pix, double _max = itk::NumericTraits<double>::max())
{
  for (unsigned int i = 0; i < 3; ++i)
  {
    pix[i] = randgen.drand64(_max);
  }
}

static double
abs_diff(const itk::RGBPixel<unsigned char> & pix1, const itk::RGBPixel<unsigned char> & pix2)
{
  double diff = 0.0;

  for (int i = 0; i < 3; ++i)
  {
    diff += itk::Math::abs(static_cast<double>(pix1[i] - pix2[i]));
  }
  return diff;
}

static double
abs_diff(const itk::RGBPixel<char> & pix1, const itk::RGBPixel<char> & pix2)
{
  double diff = 0.0;

  for (int i = 0; i < 3; ++i)
  {
    diff += itk::Math::abs(static_cast<double>(pix1[i] - pix2[i]));
  }
  return diff;
}


static double
abs_diff(const itk::Vector<float> & pix1, const itk::Vector<float> & pix2)
{
  double diff = 0.0;

  for (int i = 0; i < 3; ++i)
  {
    diff += itk::Math::abs(pix1[i] - pix2[i]);
  }
  return diff;
}

static double
abs_diff(const itk::Vector<double> & pix1, const itk::Vector<double> & pix2)
{
  double diff = 0.0;

  for (int i = 0; i < 3; ++i)
  {
    diff += itk::Math::abs(pix1[i] - pix2[i]);
  }
  return diff;
}


static double
abs_diff(const itk::Vector<int> & pix1, const itk::Vector<int> & pix2)
{
  double diff = 0.0;

  for (int i = 0; i < 3; ++i)
  {
    diff += itk::Math::abs(static_cast<double>(pix1[i] - pix2[i]));
  }
  return diff;
}

static double
abs_diff(const itk::Vector<unsigned int> & pix1, const itk::Vector<unsigned int> & pix2)
{
  double diff = 0.0;

  for (int i = 0; i < 3; ++i)
  {
    diff += itk::Math::abs(static_cast<double>(pix1[i] - pix2[i]));
  }
  return diff;
}

static double
abs_diff(const itk::Vector<short> & pix1, const itk::Vector<short> & pix2)
{
  double diff = 0.0;

  for (short i = 0; i < 3; ++i)
  {
    diff += itk::Math::abs(static_cast<double>(pix1[i] - pix2[i]));
  }
  return diff;
}

static double
abs_diff(const itk::Vector<unsigned short> & pix1, const itk::Vector<unsigned short> & pix2)
{
  double diff = 0.0;

  for (short i = 0; i < 3; ++i)
  {
    diff += itk::Math::abs(static_cast<double>(pix1[i] - pix2[i]));
  }
  return diff;
}


static void
RandomPix(vnl_random & randgen, double & pix, double _max = itk::NumericTraits<double>::max())
{
  pix = randgen.drand64(_max);
}

static void
RandomPix(vnl_random & randgen, float & pix, float _max = itk::NumericTraits<float>::max())
{
  pix = randgen.drand64(_max);
}

static void
RandomPix(vnl_random & randgen, int & pix)
{
  pix = randgen.lrand32();
}

static void
RandomPix(vnl_random & randgen, unsigned int & pix)
{
  pix = randgen.lrand32();
}

template <typename TPixel>
static double
abs_diff(const TPixel & pix1, const TPixel & pix2)
{
  return itk::Math::abs(static_cast<double>(pix1 - pix2));
}

template <typename TPixel>
static void
RandomVectorPix(vnl_random &                        randgen,
                itk::VariableLengthVector<TPixel> & pix,
                double                              _max = itk::NumericTraits<TPixel>::max())
{
  for (size_t i = 0; i < pix.GetSize(); ++i)
  {
    pix.SetElement(i, randgen.drand64(_max));
  }
}

template <typename TPixel>
static void
RandomPix(vnl_random & randgen, TPixel & pix, double _max = itk::NumericTraits<TPixel>::max())
{
  pix = randgen.lrand32((TPixel)_max);
}

template <typename TPixel>
static bool
equal(const itk::VariableLengthVector<TPixel> & pix1, const itk::VariableLengthVector<TPixel> & pix2)
{
  for (size_t i = 0; i < pix1.GetSize(); ++i)
  {
    if (pix1[i] != pix2[i])
    {
      return false;
    }
  }
  return true;
}

template <typename TPixel>
static double
abs_vector_diff(const itk::VariableLengthVector<TPixel> & pix1, const itk::VariableLengthVector<TPixel> & pix2)
{
  double diff = 0.0;

  for (size_t i = 0; i < pix1.GetSize(); ++i)
  {
    double d = itk::Math::abs(static_cast<double>(pix1[i] - pix2[i]));
    if (d > diff)
    {
      diff = d;
    }
  }
  return diff;
}

template <typename TPixel>
static double
eql_vector_diff(const itk::Point<TPixel, 3> & v1, const itk::Point<TPixel, 3> & v2)
{
  double diff = 0.0;

  for (size_t i = 0; i < 3; ++i)
  {
    diff += (v1[i] - v2[i]) * (v1[i] - v2[i]);
  }
  return sqrt(diff);
}

// TODO: properly implement storage type in MINC2 IO
template <typename TPixel, int VDimension>
int
MINCReadWriteTest(const char * fileName, const char * minc_storage_type, double tolerance = 0.0)
{
  int success(EXIT_SUCCESS);

  using ImageType = typename itk::Image<TPixel, VDimension>;

  typename ImageType::SizeType    size;
  typename ImageType::IndexType   index;
  typename ImageType::SpacingType spacing;
  typename ImageType::PointType   origin;

  std::cout << "Testing:" << fileName << std::endl;
  for (unsigned int i = 0; i < VDimension; ++i)
  {
    size[i] = 5;
    index[i] = 0;
    spacing[i] = 1.0 + static_cast<double>(i);
    origin[i] = static_cast<double>(i) * 5.0;
  }

  auto im = ImageType::New();

  typename ImageType::RegionType region;
  region.SetSize(size);
  region.SetIndex(index);
  im->SetRegions(region);
  im->SetSpacing(spacing);
  im->SetOrigin(origin);
  im->Allocate();

  itk::Matrix<itk::SpacePrecisionType, VDimension, VDimension> mat;

  mat.SetIdentity();

  if (VDimension == 3)
  { // there are problems with 4D direction cosines!
    // 30deg rotation
    mat[1][1] = mat[0][0] = 0.866025403784439;
    mat[0][1] = -0.5;
    mat[1][0] = 0.5;
    im->SetDirection(mat);
  }
  //
  // add some unique metadata
  itk::MetaDataDictionary & metaDict(im->GetMetaDataDictionary());

  itk::EncapsulateMetaData<double>(metaDict, "acquisition:TestDouble", 1.23);
  itk::EncapsulateMetaData<float>(metaDict, "acquisition:TestFloat", 1.2f);
  itk::EncapsulateMetaData<int>(metaDict, "acquisition:TestInt", 4);

  itk::Array<double> metaDataDoubleArray(5);
  metaDataDoubleArray[0] = 3.1;
  metaDataDoubleArray[1] = 1.2;
  metaDataDoubleArray[2] = 4.3;
  metaDataDoubleArray[3] = 5.4;
  metaDataDoubleArray[4] = 2.5;
  itk::EncapsulateMetaData<itk::Array<double>>(metaDict, "acquisition:TestDoubleArray", metaDataDoubleArray);

  itk::Array<float> metaDataFloatArray(5);
  metaDataFloatArray[0] = 3.1f;
  metaDataFloatArray[1] = 1.2f;
  metaDataFloatArray[2] = 4.3f;
  metaDataFloatArray[3] = 5.4f;
  metaDataFloatArray[4] = 2.5f;
  itk::EncapsulateMetaData<itk::Array<float>>(metaDict, "acquisition:TestFloatArray", metaDataFloatArray);

  itk::Array<int> metaDataIntArray(5);
  metaDataIntArray[0] = 3;
  metaDataIntArray[1] = 1;
  metaDataIntArray[2] = 4;
  metaDataIntArray[3] = 5;
  metaDataIntArray[4] = 2;
  itk::EncapsulateMetaData<itk::Array<int>>(metaDict, "acquisition:TestIntArray", metaDataIntArray);

  std::string metaDataStdString("Test std::string");
  itk::EncapsulateMetaData<std::string>(metaDict, "acquisition:StdString", metaDataStdString);

  //
  // fill image buffer
  vnl_random                          randgen(12345678);
  itk::ImageRegionIterator<ImageType> it(im, im->GetLargestPossibleRegion());

  for (it.GoToBegin(); !it.IsAtEnd(); ++it)
  {
    TPixel pix;
    if (tolerance > 0.0)
    {
      RandomPix(randgen, pix, 100);
    }
    else
    {
      RandomPix(randgen, pix);
    }
    it.Set(pix);
  }

  // set minc file storage type
  // TODO: implement this feature
  // minc::set_minc_storage_type(im,minc_storage_type,minc_storage_type!=typeid(unsigned char).name());
  itk::EncapsulateMetaData<std::string>(metaDict, "storage_data_type", minc_storage_type);

  typename ImageType::Pointer im2;
  try
  {
    itk::IOTestHelper::WriteImage<ImageType, itk::MINCImageIO>(im, std::string(fileName));
    im2 = itk::IOTestHelper::ReadImage<ImageType>(std::string(fileName));
  }
  catch (const itk::ExceptionObject & err)
  {
    std::cout << "itkMINCImageIOTest" << std::endl
              << "Exception Object caught: " << std::endl
              << err << " in " << fileName << std::endl;
    return EXIT_FAILURE;
  }

  if (eql_vector_diff(im->GetOrigin(), im2->GetOrigin()) > 1e-6)
  {
    std::cout << "Origin read " << im2->GetOrigin() << " doesn't match origin written" << im->GetOrigin() << " in "
              << fileName << std::endl;
    return EXIT_FAILURE;
  }
  if (im->GetSpacing() != im2->GetSpacing())
  {
    std::cout << "Spacing read " << im2->GetSpacing() << " doesn't match spacing written" << im->GetSpacing() << " in "
              << fileName << std::endl;
    return EXIT_FAILURE;
  }
  if (im->GetDirection() != im2->GetDirection())
  {
    std::cout << "Direction read " << im2->GetDirection() << " doesn't match direction written" << im->GetDirection()
              << " in " << fileName << std::endl;
    return EXIT_FAILURE;
  }

  //
  // Check MetaData
  itk::MetaDataDictionary & metaDict2(im2->GetMetaDataDictionary());

  double metaDataDouble = 0.0;
  if (!itk::ExposeMetaData<double>(metaDict2, "acquisition:TestDouble", metaDataDouble) || metaDataDouble != 1.23)
  {
    std::cerr << "Failure reading metaData "
              << "TestDouble " << std::endl;
    success = EXIT_FAILURE;
  }

  float metaDataFloat = 0.0f;
  if (!itk::ExposeMetaData<float>(metaDict2, "acquisition:TestFloat", metaDataFloat) || metaDataFloat != 1.2f)
  {
    std::cerr << "Failure reading metaData "
              << "acquisition:TestFloat " << std::endl;
    std::cerr << "Expected 1.2f, read=" << metaDataFloat << " difference=" << (metaDataFloat - 1.2f) << std::endl;
    success = EXIT_FAILURE;
  }

  int metaDataInt{};
  if (!itk::ExposeMetaData<int>(metaDict2, "acquisition:TestInt", metaDataInt) || metaDataInt != 4)
  {
    std::cerr << "Failure reading metaData "
              << "TestInt " << std::endl;
    success = EXIT_FAILURE;
  }

  itk::Array<double> metaDataDoubleArray2;
  if (!itk::ExposeMetaData<itk::Array<double>>(metaDict2, "acquisition:TestDoubleArray", metaDataDoubleArray2) ||
      metaDataDoubleArray2 != metaDataDoubleArray)
  {
    std::cerr << "Failure reading metaData "
              << "acquisition:TestDoubleArray " << std::endl;
    std::cerr << "metaDataDoubleArray=";
    for (size_t i = 0; i < metaDataDoubleArray.size(); ++i)
    {
      std::cerr << metaDataDoubleArray[i] << ' ';
    }
    std::cerr << std::endl;

    std::cerr << "metaDataDoubleArray2=";
    for (size_t i = 0; i < metaDataDoubleArray2.size(); ++i)
    {
      std::cerr << metaDataDoubleArray2[i] << ' ';
    }
    std::cerr << std::endl;

    success = EXIT_FAILURE;
  }

  itk::Array<float> metaDataFloatArray2;
  if (!itk::ExposeMetaData<itk::Array<float>>(metaDict2, "acquisition:TestFloatArray", metaDataFloatArray2) ||
      metaDataFloatArray2 != metaDataFloatArray)
  {
    std::cerr << "Failure reading metaData "
              << "acquisition:TestFloatArray " << std::endl;
    std::cerr << "metaDataFloatArray=";
    for (size_t i = 0; i < metaDataFloatArray.size(); ++i)
    {
      std::cerr << metaDataFloatArray[i] << ' ';
    }
    std::cerr << std::endl;

    std::cerr << "metaDataFloatArray2=";
    for (size_t i = 0; i < metaDataFloatArray2.size(); ++i)
    {
      std::cerr << metaDataFloatArray2[i] << ' ';
    }
    std::cerr << std::endl;

    success = EXIT_FAILURE;
  }

  itk::Array<int> metaDataIntArray2;
  if (!itk::ExposeMetaData<itk::Array<int>>(metaDict2, "acquisition:TestIntArray", metaDataIntArray2) ||
      metaDataIntArray2 != metaDataIntArray)
  {
    std::cerr << "Failure reading metaData "
              << "acquisition:TestIntArray " << std::endl;
    success = EXIT_FAILURE;
  }

  std::string metaDataStdString2("");
  if (!itk::ExposeMetaData<std::string>(metaDict2, "acquisition:StdString", metaDataStdString2) ||
      metaDataStdString2 != metaDataStdString)
  {
    std::cerr << "Failure reading metaData "
              << "StdString " << metaDataStdString2 << ' ' << metaDataStdString << std::endl;
    success = EXIT_FAILURE;
  }

  itk::ImageRegionIterator<ImageType> it2(im2, im2->GetLargestPossibleRegion());
  if (tolerance == 0.0)
  {
    for (it.GoToBegin(), it2.GoToBegin(); !it.IsAtEnd() && !it2.IsAtEnd(); ++it, ++it2)
    {
      if (it.Value() != it2.Value())
      {
        std::cout << "Original Pixel (" << it.Value() << ") doesn't match read-in Pixel (" << it2.Value() << " ) "
                  << std::endl
                  << " in " << fileName << std::endl;
        success = EXIT_FAILURE;
        break;
      }
    }
  }
  else
  { // account for rounding errors
    for (it.GoToBegin(), it2.GoToBegin(); !it.IsAtEnd() && !it2.IsAtEnd(); ++it, ++it2)
    {
      if (abs_diff(it.Value(), it2.Value()) > tolerance)
      {
        std::cout << "Original Pixel (" << it.Value() << ") doesn't match read-in Pixel (" << it2.Value() << " ) "
                  << " in " << fileName << std::endl;
        success = EXIT_FAILURE;
        break;
      }
    }
  }
  // Remove(fileName);
  return success;
}

template <typename TPixel, int VDimension>
int
MINCReadWriteTestVector(const char * fileName,
                        size_t       vector_length,
                        const char * minc_storage_type,
                        double       tolerance = 0.0)
{
  int success(EXIT_SUCCESS);

  using ImageType = typename itk::VectorImage<TPixel, VDimension>;
  using InternalPixelType = typename itk::VectorImage<TPixel, VDimension>::PixelType;

  typename ImageType::SizeType    size;
  typename ImageType::IndexType   index;
  typename ImageType::SpacingType spacing;
  typename ImageType::PointType   origin;

  std::cout << "Testing:" << fileName << std::endl;

  for (unsigned int i = 0; i < VDimension; ++i)
  {
    size[i] = 5;
    index[i] = 0;
    spacing[i] = 1.0 + static_cast<double>(i);
    origin[i] = static_cast<double>(i) * 5.0;
  }
  auto im = ImageType::New();

  // itk::IOTestHelper::AllocateImageFromRegionAndSpacing<ImageType>(imageRegion,spacing);
  typename ImageType::RegionType region;
  region.SetSize(size);
  region.SetIndex(index);
  im->SetRegions(region);
  im->SetSpacing(spacing);
  im->SetOrigin(origin);
  im->SetVectorLength(vector_length);
  im->Allocate();

  itk::Matrix<itk::SpacePrecisionType, VDimension, VDimension> mat;

  mat.SetIdentity();

  if (VDimension == 3)
  { // there are problems with 4D direction cosines!
    // 30deg rotation
    mat[1][1] = mat[0][0] = 0.866025403784439;
    mat[0][1] = -0.5;
    mat[1][0] = 0.5;
    im->SetDirection(mat);
  }
  //
  // add some unique metadata
  itk::MetaDataDictionary & metaDict(im->GetMetaDataDictionary());

  itk::EncapsulateMetaData<double>(metaDict, "acquisition:TestDouble", 1.23);
  itk::EncapsulateMetaData<int>(metaDict, "acquisition:TestInt", 4);

  itk::Array<double> metaDataDoubleArray(5);
  metaDataDoubleArray[0] = 3.1;
  metaDataDoubleArray[1] = 1.2;
  metaDataDoubleArray[2] = 4.3;
  metaDataDoubleArray[3] = 5.4;
  metaDataDoubleArray[4] = 2.5;
  itk::EncapsulateMetaData<itk::Array<double>>(metaDict, "acquisition:TestDoubleArray", metaDataDoubleArray);

  itk::Array<int> metaDataIntArray(5);
  metaDataIntArray[0] = 3;
  metaDataIntArray[1] = 1;
  metaDataIntArray[2] = 4;
  metaDataIntArray[3] = 5;
  metaDataIntArray[4] = 2;
  itk::EncapsulateMetaData<itk::Array<int>>(metaDict, "acquisition:TestIntArray", metaDataIntArray);

  std::string metaDataStdString("Test std::string");
  itk::EncapsulateMetaData<std::string>(metaDict, "acquisition:StdString", metaDataStdString);

  //
  // fill image buffer
  vnl_random                          randgen(12345678);
  itk::ImageRegionIterator<ImageType> it(im, im->GetLargestPossibleRegion());

  for (it.GoToBegin(); !it.IsAtEnd(); ++it)
  {
    InternalPixelType pix(vector_length);
    if (tolerance > 0.0)
    {
      RandomVectorPix<TPixel>(randgen, pix, 100.0);
    }
    else
    {
      RandomVectorPix<TPixel>(randgen, pix);
    }
    it.Set(pix);
  }

  // set minc file storage type
  // minc::set_minc_storage_type(im,minc_storage_type,minc_storage_type!=typeid(unsigned char).name());
  itk::EncapsulateMetaData<std::string>(metaDict, "storage_data_type", minc_storage_type);

  typename ImageType::Pointer im2;

  try
  {
    itk::IOTestHelper::WriteImage<ImageType, itk::MINCImageIO>(im, std::string(fileName));
    im2 = itk::IOTestHelper::ReadImage<ImageType>(std::string(fileName));
  }
  catch (const itk::ExceptionObject & err)
  {
    std::cout << "itkMINCImageIOTest" << std::endl
              << "Exception Object caught: " << std::endl
              << err << " in " << fileName << std::endl;
    return EXIT_FAILURE;
  }

  if (eql_vector_diff(im->GetOrigin(), im2->GetOrigin()) > 1e-6) // account for
                                                                 // rounding errors
  {
    std::cout << "Origin read " << im2->GetOrigin() << " doesn't match origin written" << im->GetOrigin() << " in "
              << fileName << std::endl;
    return EXIT_FAILURE;
  }
  if (im->GetSpacing() != im2->GetSpacing())
  {
    std::cout << "Spacing read " << im2->GetSpacing() << " doesn't match spacing written" << im->GetSpacing() << " in "
              << fileName << std::endl;
    return EXIT_FAILURE;
  }
  if (im->GetDirection() != im2->GetDirection())
  {
    std::cout << "Direction read " << im2->GetDirection() << " doesn't match direction written" << im->GetDirection()
              << " in " << fileName << std::endl;
    return EXIT_FAILURE;
  }

  // Check MetaData
  itk::MetaDataDictionary & metaDict2(im2->GetMetaDataDictionary());

  double metaDataDouble = 0.0;
  if (!itk::ExposeMetaData<double>(metaDict2, "acquisition:TestDouble", metaDataDouble) || metaDataDouble != 1.23)
  {
    std::cerr << "Failure reading metaData "
              << "TestDouble " << std::endl;
    success = EXIT_FAILURE;
  }

  int metaDataInt{};
  if (!itk::ExposeMetaData<int>(metaDict2, "acquisition:TestInt", metaDataInt) || metaDataInt != 4)
  {
    std::cerr << "Failure reading metaData "
              << "TestInt " << std::endl;
    success = EXIT_FAILURE;
  }

  itk::Array<double> metaDataDoubleArray2;
  if (!itk::ExposeMetaData<itk::Array<double>>(metaDict2, "acquisition:TestDoubleArray", metaDataDoubleArray2) ||
      metaDataDoubleArray2 != metaDataDoubleArray)
  {
    std::cerr << "Failure reading metaData "
              << "TestDoubleArray " << std::endl;
    success = EXIT_FAILURE;
  }

  itk::Array<int> metaDataIntArray2;
  if (!itk::ExposeMetaData<itk::Array<int>>(metaDict2, "acquisition:TestIntArray", metaDataIntArray2) ||
      metaDataIntArray2 != metaDataIntArray)
  {
    std::cerr << "Failure reading metaData "
              << "TestIntArray " << std::endl;
    success = EXIT_FAILURE;
  }

  std::string metaDataStdString2("");
  if (!itk::ExposeMetaData<std::string>(metaDict2, "acquisition:StdString", metaDataStdString2) ||
      metaDataStdString2 != metaDataStdString)
  {
    std::cerr << "Failure reading metaData "
              << "StdString " << metaDataStdString2 << ' ' << metaDataStdString << std::endl;
    success = EXIT_FAILURE;
  }

  itk::ImageRegionIterator<ImageType> it2(im2, im2->GetLargestPossibleRegion());
  InternalPixelType                   pix1, pix2;
  if (tolerance == 0.0)
  {
    for (it.GoToBegin(), it2.GoToBegin(); !it.IsAtEnd() && !it2.IsAtEnd(); ++it, ++it2)
    {
      pix1 = it.Get();
      pix2 = it2.Get();
      if (!equal<TPixel>(pix1, pix2))
      {
        std::cout << "Original Pixel (" << pix1 << ") doesn't match read-in Pixel (" << pix2 << " ) "
                  << " in " << fileName << std::endl;
        success = EXIT_FAILURE;
        break;
      }
    }
  }
  else
  { // account for rounding errors
    for (it.GoToBegin(), it2.GoToBegin(); !it.IsAtEnd() && !it2.IsAtEnd(); ++it, ++it2)
    {
      pix1 = it.Get();
      pix2 = it2.Get();
      if (abs_vector_diff<TPixel>(pix1, pix2) > tolerance)
      {
        std::cout << "Original Pixel (" << pix1 << ") doesn't match read-in Pixel (" << pix2 << " ) "
                  << " in " << fileName << std::endl;
        success = EXIT_FAILURE;
        break;
      }
    }
  }
  // Remove(fileName);
  return success;
}

int
itkMINCImageIOTest(int argc, char * argv[])
{
  std::string prefix("");

  if (argc > 1)
  {
    prefix = *++argv;
    --argc;
    itksys::SystemTools::ChangeDirectory(prefix.c_str());
  }

  itk::ObjectFactoryBase::RegisterFactory(itk::MINCImageIOFactory::New(),
                                          itk::ObjectFactoryBase::InsertionPositionEnum::INSERT_AT_FRONT);

  int result(0);
  // stright forward test
  result += MINCReadWriteTest<unsigned char, 3>("3DUCharImage.mnc", typeid(unsigned char).name());
  result += MINCReadWriteTest<char, 3>("3DCharImage.mnc", typeid(unsigned char).name());
  result += MINCReadWriteTest<unsigned short, 3>("3DUShortImage.mnc", typeid(short).name());
  result += MINCReadWriteTest<short, 3>("3DShortImage.mnc", typeid(unsigned short).name());
  result += MINCReadWriteTest<unsigned int, 3>("3DUIntImage.mnc", typeid(int).name());
  result += MINCReadWriteTest<int, 3>("3DIntImage.mnc", typeid(unsigned int).name());
  result += MINCReadWriteTest<float, 3>("3DFloatImage.mnc", typeid(float).name());
  result += MINCReadWriteTest<double, 3>("3DDoubleImage.mnc", typeid(double).name());

  result += MINCReadWriteTest<itk::RGBPixel<unsigned char>, 3>("3DRGBImage.mnc", typeid(unsigned char).name());
  result += MINCReadWriteTest<itk::RGBPixel<char>, 3>("3DRGBImage2.mnc", typeid(unsigned char).name());
  result += MINCReadWriteTest<itk::Vector<float, 3>, 3>("3DVectorImage_float.mnc", typeid(float).name());
  result += MINCReadWriteTest<itk::Vector<double, 3>, 3>("3DVectorImage_double.mnc", typeid(double).name());

  result += MINCReadWriteTest<itk::Vector<short, 3>, 3>("3DVectorImage_short.mnc", typeid(short).name());
  result +=
    MINCReadWriteTest<itk::Vector<unsigned short, 3>, 3>("3DVectorImage_ushort.mnc", typeid(unsigned short).name());

  result += MINCReadWriteTest<itk::Vector<int, 3>, 3>("3DVectorImage_int.mnc", typeid(int).name());
  result += MINCReadWriteTest<itk::Vector<unsigned int, 3>, 3>("3DVectorImage_uint.mnc", typeid(unsigned int).name());

  // expecting rounding errors
  result += MINCReadWriteTest<float, 3>("3DFloatImage_byte.mnc", typeid(unsigned char).name(), 0.2);
  result += MINCReadWriteTest<float, 3>("3DFloatImage_short.mnc", typeid(short).name(), 0.01);

  result += MINCReadWriteTest<double, 3>("3DDoubleImage_double_byte.mnc", typeid(unsigned char).name(), 0.2);
  result += MINCReadWriteTest<double, 3>("3DDoubleImage_double_short.mnc", typeid(short).name(), 0.01);

  result += MINCReadWriteTest<itk::Vector<float, 3>, 3>(
    "3DVectorImage_float_byte.mnc", typeid(unsigned char).name(), 0.5 * sqrt(3.0));
  result += MINCReadWriteTest<itk::Vector<float, 3>, 3>(
    "3DVectorImage_float_short.mnc", typeid(short).name(), 0.05 * sqrt(3.0));

  // testing variable vector case
  // stright forward test
  result += MINCReadWriteTestVector<unsigned char, 3>("4DUCharImage.mnc", 10, typeid(unsigned char).name(), 0.0001);
  result += MINCReadWriteTestVector<float, 3>("4DFloatImage.mnc", 10, typeid(float).name(), 0.0001);
  result += MINCReadWriteTestVector<double, 3>("4DDoubleImage.mnc", 10, typeid(double).name(), 0.0001);

  // expecting rounding errors
  result += MINCReadWriteTestVector<float, 3>("4DFloatImage_byte.mnc", 10, typeid(unsigned char).name(), 0.2);
  result += MINCReadWriteTestVector<float, 3>("4DFloatImage_short.mnc", 10, typeid(short).name(), 0.01);

  result += MINCReadWriteTestVector<double, 3>("4DDoubleImage_byte.mnc", 10, typeid(unsigned char).name(), 0.2);
  result += MINCReadWriteTestVector<double, 3>("4DDoubleImage_short.mnc", 10, typeid(short).name(), 0.01);

  //   result += MINCReadWriteTest<unsigned char,4>("4DUCharImage.mnc",typeid(unsigned char).name(),0.5);
  //   result += MINCReadWriteTest<float,4>("4DFloatImage.mnc",typeid(float).name(),0.0001);
  //
  //   result += MINCReadWriteTest<itk::RGBPixel<unsigned char>,4 >("4DRGBImage.mnc",typeid(unsigned char).name(),0.5);
  //   result += MINCReadWriteTest<itk::Vector<float,3>,3 >("4DVectorImage.mnc",typeid(float).name(),0.0001);

  return result != 0;
}