ENH: Support genuine 4D (x,y,z,time) MINC images

The MINC time dimension now maps to a 4th ITK dimension; the
vector_dimension maps to components. Component count no longer selects
the written dimension kind, so time and vector survive a round-trip.

Change-Id: I155561693ee98c5677b2dfa125fc9990e7721bb2

Author: gdevenyi

Modules/IO/MINC/include/itkMINCImageIO.h

@@ -85,11 +85,12 @@ class ITKIOMINC_EXPORT MINCImageIO : public ImageIOBase
   /** \see LightObject::GetNameOfClass() */
   itkOverrideGetNameOfClassMacro(MINCImageIO);
 
-  /** Right now MINC supports up to 3D with multiple components */
+  /** MINC supports up to 3 spatial dimensions plus a time dimension (4D),
+   * each with multiple components (MINC vector_dimension). */
   bool
   SupportsDimension(unsigned long dim) override
   {
-    return dim < 4;
+    return dim < 5;
   }
 
   /*-------- This part of the interface deals with reading data. ------ */

Modules/IO/MINC/src/itkMINCImageIO.cxx

@@ -440,12 +440,14 @@ MINCImageIO::ReadImageInformation()
     itkExceptionStringMacro(" minc files without spatial dimensions are not supported!");
   }
 
-  if (m_MINCPImpl->m_DimensionIndices[0] != -1 && m_MINCPImpl->m_DimensionIndices[4] != -1)
-  {
-    itkExceptionStringMacro(" 4D minc files vector dimension are not supported currently");
-  }
+  const bool haveTimeDimension = (m_MINCPImpl->m_DimensionIndices[4] != -1);
 
-  this->SetNumberOfDimensions(spatial_dimension_count);
+  // The MINC time dimension becomes an additional (highest) ITK dimension, while
+  // the MINC vector_dimension is mapped to ITK components. The two are
+  // independent and may coexist (e.g. a 4D vector image).
+  const unsigned int itkDimensionCount = spatial_dimension_count + (haveTimeDimension ? 1 : 0);
+
+  this->SetNumberOfDimensions(itkDimensionCount);
 
   int          numberOfComponents = 1;
   unsigned int usableDimensions = 0;
@@ -457,11 +459,37 @@ MINCImageIO::ReadImageInformation()
   auto RAStofromLPS = Matrix<double, 3, 3>::GetIdentity();
   RAStofromLPS(0, 0) = -1.0;
   RAStofromLPS(1, 1) = -1.0;
-  std::vector<double> dir_cos_temp(3);
+  std::vector<double> dir_cos_temp(itkDimensionCount, 0.0);
 
   Vector<double, 3> origin{};
   Vector<double, 3> oOrigin{};
 
+  // Build the MINC apparent dimension order slowest-varying first. The time
+  // dimension is the slowest, matching the ITK buffer layout used by
+  // Read()/Write() where it is the highest ITK dimension.
+  if (haveTimeDimension)
+  {
+    m_MINCPImpl->m_MincApparentDims[usableDimensions] = m_MINCPImpl->m_MincFileDims[m_MINCPImpl->m_DimensionIndices[4]];
+    // always use positive
+    miset_dimension_apparent_voxel_order(m_MINCPImpl->m_MincApparentDims[usableDimensions], MI_POSITIVE);
+    misize_t _sz = 0;
+    miget_dimension_size(m_MINCPImpl->m_MincApparentDims[usableDimensions], &_sz);
+    double _sep = NAN;
+    miget_dimension_separation(m_MINCPImpl->m_MincApparentDims[usableDimensions], MI_ORDER_APPARENT, &_sep);
+    double _start = NAN;
+    miget_dimension_start(m_MINCPImpl->m_MincApparentDims[usableDimensions], MI_ORDER_APPARENT, &_start);
+
+    this->SetDimensions(spatial_dimension_count, static_cast<unsigned int>(_sz));
+    this->SetSpacing(spatial_dimension_count, _sep);
+    this->SetOrigin(spatial_dimension_count, _start);
+
+    std::vector<double> time_dir(itkDimensionCount, 0.0);
+    time_dir[spatial_dimension_count] = 1.0;
+    this->SetDirection(spatial_dimension_count, time_dir);
+
+    ++usableDimensions;
+  }
+
   // minc api uses inverse order of dimensions , fastest varying are last
   Vector<double, 3> sep;
   for (int i = 3; i > 0; i--)
@@ -511,7 +539,7 @@ MINCImageIO::ReadImageInformation()
   for (int i = 0; i < spatial_dimension_count; ++i)
   {
     this->SetOrigin(i, oOrigin[i]);
-    for (unsigned int j = 0; j < 3; j++)
+    for (unsigned int j = 0; j < static_cast<unsigned int>(spatial_dimension_count); ++j)
     {
       dir_cos_temp[j] = dir_cos[j][i];
     }
@@ -530,18 +558,6 @@ MINCImageIO::ReadImageInformation()
     ++usableDimensions;
   }
 
-  if (m_MINCPImpl->m_DimensionIndices[4] != -1) // have time dimension
-  {
-    // micopy_dimension(hdim[m_MINCPImpl->m_DimensionIndices[4]],&apparent_dimension_order[usable_dimensions]);
-    m_MINCPImpl->m_MincApparentDims[usableDimensions] = m_MINCPImpl->m_MincFileDims[m_MINCPImpl->m_DimensionIndices[4]];
-    // always use positive
-    miset_dimension_apparent_voxel_order(m_MINCPImpl->m_MincApparentDims[usableDimensions], MI_POSITIVE);
-    misize_t _sz = 0;
-    miget_dimension_size(m_MINCPImpl->m_MincApparentDims[usableDimensions], &_sz);
-    numberOfComponents = _sz;
-    ++usableDimensions;
-  }
-
   // Set apparent dimension order to the MINC2 api
   if (miset_apparent_dimension_order(m_MINCPImpl->m_Volume, usableDimensions, m_MINCPImpl->m_MincApparentDims) < 0)
   {
@@ -893,33 +909,9 @@ MINCImageIO::WriteImageInformation()
   MetaDataDictionary & thisDic = GetMetaDataDictionary();
 
   unsigned int minc_dimensions = 0;
-  if (nComp > 3) // last dimension will be either vector or time
-  {
-    micreate_dimension(MItime,
-                       MI_DIMCLASS_TIME,
-                       MI_DIMATTR_REGULARLY_SAMPLED,
-                       nComp,
-                       &m_MINCPImpl->m_MincApparentDims[m_MINCPImpl->m_NDims - minc_dimensions - 1]);
-
-    double tstart = 0.0;
-    if (!ExposeMetaData<double>(thisDic, "tstart", tstart))
-    {
-      tstart = 0.0;
-    }
-
-    miset_dimension_start(m_MINCPImpl->m_MincApparentDims[m_MINCPImpl->m_NDims - minc_dimensions - 1], tstart);
-
-    double tstep = 1.0;
-    if (!ExposeMetaData<double>(thisDic, "tstep", tstep))
-    {
-      tstep = 1.0;
-    }
-
-    miset_dimension_separation(m_MINCPImpl->m_MincApparentDims[m_MINCPImpl->m_NDims - minc_dimensions - 1], tstep);
-
-    ++minc_dimensions;
-  }
-  else if (nComp > 1)
+  // ITK components are always written as a MINC vector_dimension (fastest
+  // varying). A time series is an ITK dimension, not components (see below).
+  if (nComp > 1)
   {
     micreate_dimension(MIvector_dimension,
                        MI_DIMCLASS_RECORD,
@@ -956,10 +948,20 @@ MINCImageIO::WriteImageInformation()
     ++minc_dimensions;
   }
 
-  if (nDims > 3)
+  if (nDims > 4)
   {
     MINCIOFreeTmpDimHandle(minc_dimensions, m_MINCPImpl->m_MincApparentDims);
-    itkExceptionStringMacro("Unfortunately, only up to 3D volume are supported now.");
+    itkExceptionStringMacro("Unfortunately, only up to 4D volumes (3 spatial + time) are supported.");
+  }
+
+  if (nDims > 3) // the 4th ITK dimension is written as a MINC time dimension (slowest varying)
+  {
+    micreate_dimension(MItime,
+                       MI_DIMCLASS_TIME,
+                       MI_DIMATTR_REGULARLY_SAMPLED,
+                       this->GetDimensions(3),
+                       &m_MINCPImpl->m_MincApparentDims[m_MINCPImpl->m_NDims - minc_dimensions - 1]);
+    ++minc_dimensions;
   }
 
   // allocating dimensions
@@ -994,21 +996,18 @@ MINCImageIO::WriteImageInformation()
   for (unsigned int i = 0; i < nDims; ++i)
   {
     const unsigned int j = i + (nComp > 1 ? 1 : 0);
-    double             dir_cos[3];
-    for (unsigned int k = 0; k < 3; ++k)
+    const unsigned int apparentIndex = minc_dimensions - j - 1;
+    miset_dimension_separation(m_MINCPImpl->m_MincApparentDims[apparentIndex], this->GetSpacing(i));
+    miset_dimension_start(m_MINCPImpl->m_MincApparentDims[apparentIndex], origin[i]);
+    if (i < 3) // spatial dimension: also set the direction cosines (time has none)
     {
-      if (k < nDims)
+      double dir_cos[3];
+      for (unsigned int k = 0; k < 3; ++k)
       {
-        dir_cos[k] = directionCosineMatrix[i][k];
-      }
-      else
-      {
-        dir_cos[k] = 0.0;
+        dir_cos[k] = (k < nDims) ? directionCosineMatrix[i][k] : 0.0;
       }
+      miset_dimension_cosines(m_MINCPImpl->m_MincApparentDims[apparentIndex], dir_cos);
     }
-    miset_dimension_separation(m_MINCPImpl->m_MincApparentDims[minc_dimensions - j - 1], this->GetSpacing(i));
-    miset_dimension_start(m_MINCPImpl->m_MincApparentDims[minc_dimensions - j - 1], origin[i]);
-    miset_dimension_cosines(m_MINCPImpl->m_MincApparentDims[minc_dimensions - j - 1], dir_cos);
   }
 
   // TODO: fix this to appropriate
@@ -1107,14 +1106,14 @@ MINCImageIO::WriteImageInformation()
             break;
           case 't':
           case 'T':
-            if (nComp <= 1)
+            if (nDims <= 3) // no time dimension present
             {
               itkDebugMacro("Dimension order is incorrect " << dimension_order);
               dimorder_good = false;
             }
             else
             {
-              j = m_MINCPImpl->m_NDims - 1;
+              j = 0; // time is the slowest-varying apparent dimension
             }
             break;
           case 'x':

Modules/IO/MINC/test/CMakeLists.txt

@@ -129,6 +129,7 @@ itk_add_test(
     ${ITK_TEST_OUTPUT_DIR}/dti_sample.mnc
   ITK_REMOVE_TEMPORARY_TEST_FILES
     ${ITK_TEST_OUTPUT_DIR}/dti_sample.mnc
+    ${ITK_TEST_OUTPUT_DIR}/dti_sample.mnc_synthetic.mnc
 )
 
 itk_add_test(

Modules/IO/MINC/test/itkMINCImageIOTest2.cxx

@@ -45,7 +45,7 @@ itkMINCImageIOTest2(int argc, char * argv[])
   ITK_EXERCISE_BASIC_OBJECT_METHODS(mincIO1, MINCImageIO, ImageIOBase);
 
 
-  constexpr unsigned int     supportedDimCount{ 4 }; // includes the degenerate 0-dimensional case
+  constexpr unsigned int     supportedDimCount{ 5 }; // 3 spatial + time, plus the degenerate 0-dimensional case
   std::vector<unsigned long> supportedDims(supportedDimCount);
   std::iota(std::begin(supportedDims), std::end(supportedDims), 0);
   for (const auto & value : supportedDims)

Modules/IO/MINC/test/itkMINCImageIOTest_4D.cxx

@@ -16,12 +16,14 @@
  *
  *=========================================================================*/
 
-#include "itkImageAlgorithm.h"
 #include <iostream>
 
 #include "itkMINCImageIOFactory.h"
 #include "itkImageFileReader.h"
 #include "itkImageFileWriter.h"
+#include "itkImageRegionIterator.h"
+#include "itkImageRegionConstIterator.h"
+#include "itkMath.h"
 #include "itkTestingMacros.h"
 
 
@@ -38,7 +40,10 @@ itkMINCImageIOTest_4D(int argc, char * argv[])
 
   itk::MINCImageIOFactory::RegisterOneFactory();
 
-  using ImageType = itk::VectorImage<float, 3>;
+  // A MINC volume with x/y/z spatial dimensions plus a time dimension is read as
+  // a genuine 4D itk::Image, with the time step/start mapped to the 4th
+  // dimension spacing/origin.
+  using ImageType = itk::Image<float, 4>;
 
   using ReaderType = itk::ImageFileReader<ImageType>;
   using WriterType = itk::ImageFileWriter<ImageType>;
@@ -52,11 +57,72 @@ itkMINCImageIOTest_4D(int argc, char * argv[])
 
   ITK_TRY_EXPECT_NO_EXCEPTION(writer->Update());
 
-
   const ImageType::ConstPointer image = reader->GetOutput();
-
   image->Print(std::cout);
 
+  // Round-trip a synthetic 4D volume with a non-trivial time step/start and
+  // verify the time geometry survives, together with the voxel values.
+  const std::string synthetic = std::string(argv[2]) + "_synthetic.mnc";
+
+  auto            synthImage = ImageType::New();
+  constexpr float timeStep = 2.5F;
+  constexpr float timeStart = 10.0F;
+
+  ImageType::RegionType region;
+  region.SetSize({ { 3, 4, 5, 2 } });
+  synthImage->SetRegions(region);
+
+  ImageType::SpacingType spacing;
+  spacing[0] = 1.0;
+  spacing[1] = 1.0;
+  spacing[2] = 1.0;
+  spacing[3] = timeStep;
+  synthImage->SetSpacing(spacing);
+
+  ImageType::PointType origin;
+  origin[0] = 0.0;
+  origin[1] = 0.0;
+  origin[2] = 0.0;
+  origin[3] = timeStart;
+  synthImage->SetOrigin(origin);
+
+  synthImage->Allocate();
+
+  {
+    float                               v = 0.0F;
+    itk::ImageRegionIterator<ImageType> it(synthImage, synthImage->GetLargestPossibleRegion());
+    for (; !it.IsAtEnd(); ++it, ++v)
+    {
+      it.Set(v);
+    }
+  }
+
+  auto synthWriter = WriterType::New();
+  synthWriter->SetFileName(synthetic);
+  synthWriter->SetInput(synthImage);
+  ITK_TRY_EXPECT_NO_EXCEPTION(synthWriter->Update());
+
+  auto synthReader = ReaderType::New();
+  synthReader->SetFileName(synthetic);
+  ITK_TRY_EXPECT_NO_EXCEPTION(synthReader->Update());
+
+  const ImageType::ConstPointer synthBack = synthReader->GetOutput();
+
+  ITK_TEST_EXPECT_EQUAL(synthBack->GetLargestPossibleRegion().GetSize()[3], 2);
+  ITK_TEST_EXPECT_TRUE(itk::Math::FloatAlmostEqual(static_cast<float>(synthBack->GetSpacing()[3]), timeStep));
+  ITK_TEST_EXPECT_TRUE(itk::Math::FloatAlmostEqual(static_cast<float>(synthBack->GetOrigin()[3]), timeStart));
+
+  itk::ImageRegionConstIterator<ImageType> inIt(synthImage, synthImage->GetLargestPossibleRegion());
+  itk::ImageRegionConstIterator<ImageType> outIt(synthBack, synthBack->GetLargestPossibleRegion());
+  for (; !inIt.IsAtEnd(); ++inIt, ++outIt)
+  {
+    if (itk::Math::NotAlmostEquals(inIt.Get(), outIt.Get()))
+    {
+      std::cerr << "Voxel mismatch after 4D round-trip at " << inIt.GetIndex() << ": " << inIt.Get()
+                << " != " << outIt.Get() << std::endl;
+      return EXIT_FAILURE;
+    }
+  }
 
   std::cout << "Test finished." << std::endl;
   return EXIT_SUCCESS;