BUG: Fix incorrect reading of MINC2 files, convert RAS coordinates to LPS

Modules/IO/MINC/include/itkMINCImageIO.h

@@ -123,6 +123,11 @@ class ITKIOMINC_EXPORT MINCImageIO : public ImageIOBase
   void
   Write(const void * buffer) override;
 
+  /** Set to automatically convert from PositiveCoordinateOrientation RAS to PositiveCoordinateOrientation LPS*/
+  itkSetMacro(RAStoLPS, bool);
+  itkGetConstMacro(RAStoLPS, bool);
+  itkBooleanMacro(RAStoLPS);
+
 protected:
   MINCImageIO();
   ~MINCImageIO() override;
@@ -153,6 +158,7 @@ class ITKIOMINC_EXPORT MINCImageIO : public ImageIOBase
 
   // complex type images, composed of complex numbers
   // int m_Complex;
+  bool m_RAStoLPS;
 };
 } // end namespace itk
 

Modules/IO/MINC/src/CMakeLists.txt

@@ -1,3 +1,4 @@
+
 set(ITKIOMINC_SRC itkMINCImageIO.cxx itkMINCImageIOFactory.cxx)
 
 add_library(ITKIOMINC ${ITK_LIBRARY_BUILD_TYPE} ${ITKIOMINC_SRC})

Modules/IO/MINC/src/itkMINCImageIO.cxx

@@ -29,6 +29,7 @@
 
 #include <memory> // For unique_ptr.
 
+
 extern "C"
 {
   void
@@ -235,6 +236,7 @@ MINCImageIO::CloseVolume()
 
 MINCImageIO::MINCImageIO()
   : m_MINCPImpl(std::make_unique<MINCImageIOPImpl>())
+  , m_RAStoLPS(false)
 {
   for (auto & dimensionIndex : m_MINCPImpl->m_DimensionIndices)
   {
@@ -280,6 +282,7 @@ MINCImageIO::PrintSelf(std::ostream & os, Indent indent) const
 
   os << indent << "MINCPImpl: " << m_MINCPImpl.get() << std::endl;
   os << indent << "DirectionCosines: " << m_DirectionCosines << std::endl;
+  os << indent << "RAStoLPS: " << m_RAStoLPS << std::endl;
 }
 
 void
@@ -446,14 +449,20 @@ MINCImageIO::ReadImageInformation()
 
   this->SetNumberOfDimensions(spatial_dimension_count);
 
-  int numberOfComponents = 1;
-  int usable_dimensions = 0;
+  int          numberOfComponents = 1;
+  unsigned int usableDimensions = 0;
+
+  auto dir_cos = Matrix<double, 3, 3>::GetIdentity();
 
-  Matrix<double, 3, 3> dir_cos{};
-  dir_cos.SetIdentity();
+  // Conversion matrix for MINC PositiveCoordinateOrientation RAS (LeftToRight, PosteriorToAnterior, InferiorToSuperior)
+  // to ITK PositiveCoordinateOrientation LPS (RightToLeft, AnteriorToPosterior, InferiorToSuperior)
+  auto RAStofromLPS = Matrix<double, 3, 3>::GetIdentity();
+  RAStofromLPS(0, 0) = -1.0;
+  RAStofromLPS(1, 1) = -1.0;
+  std::vector<double> dir_cos_temp(3);
 
   Vector<double, 3> origin{};
-  Vector<double, 3> o_origin{};
+  Vector<double, 3> oOrigin{};
 
   // minc api uses inverse order of dimensions , fastest varying are last
   Vector<double, 3> sep;
@@ -463,19 +472,19 @@ MINCImageIO::ReadImageInformation()
     {
       // MINC2: bad design!
       // micopy_dimension(hdim[m_MINCPImpl->m_DimensionIndices[i]],&apparent_dimension_order[usable_dimensions]);
-      m_MINCPImpl->m_MincApparentDims[usable_dimensions] =
+      m_MINCPImpl->m_MincApparentDims[usableDimensions] =
         m_MINCPImpl->m_MincFileDims[m_MINCPImpl->m_DimensionIndices[i]];
       // always use positive
-      miset_dimension_apparent_voxel_order(m_MINCPImpl->m_MincApparentDims[usable_dimensions], MI_POSITIVE);
+      miset_dimension_apparent_voxel_order(m_MINCPImpl->m_MincApparentDims[usableDimensions], MI_POSITIVE);
       misize_t _sz;
-      miget_dimension_size(m_MINCPImpl->m_MincApparentDims[usable_dimensions], &_sz);
+      miget_dimension_size(m_MINCPImpl->m_MincApparentDims[usableDimensions], &_sz);
 
       double _sep;
-      miget_dimension_separation(m_MINCPImpl->m_MincApparentDims[usable_dimensions], MI_ORDER_APPARENT, &_sep);
+      miget_dimension_separation(m_MINCPImpl->m_MincApparentDims[usableDimensions], MI_ORDER_APPARENT, &_sep);
       std::vector<double> _dir(3);
-      miget_dimension_cosines(m_MINCPImpl->m_MincApparentDims[usable_dimensions], &_dir[0]);
+      miget_dimension_cosines(m_MINCPImpl->m_MincApparentDims[usableDimensions], &_dir[0]);
       double _start;
-      miget_dimension_start(m_MINCPImpl->m_MincApparentDims[usable_dimensions], MI_ORDER_APPARENT, &_start);
+      miget_dimension_start(m_MINCPImpl->m_MincApparentDims[usableDimensions], MI_ORDER_APPARENT, &_start);
 
       for (int j = 0; j < 3; ++j)
       {
@@ -486,52 +495,61 @@ MINCImageIO::ReadImageInformation()
       sep[i - 1] = _sep;
 
       this->SetDimensions(i - 1, static_cast<unsigned int>(_sz));
-      this->SetDirection(i - 1, _dir);
       this->SetSpacing(i - 1, _sep);
 
-      ++usable_dimensions;
+      ++usableDimensions;
+    }
+  }
+
+
+  // Transform MINC PositiveCoordinateOrientation RAS coordinates to
+  // internal ITK PositiveCoordinateOrientation LPS Coordinates
+  if (this->m_RAStoLPS)
+    dir_cos = RAStofromLPS * dir_cos;
+
+  // Transform origin coordinates
+  oOrigin = dir_cos * origin;
+
+  for (int i = 0; i < spatial_dimension_count; ++i)
+  {
+    this->SetOrigin(i, oOrigin[i]);
+    for (unsigned int j = 0; j < 3; j++)
+    {
+      dir_cos_temp[j] = dir_cos[j][i];
     }
+    this->SetDirection(i, dir_cos_temp);
   }
 
   if (m_MINCPImpl->m_DimensionIndices[0] != -1) // have vector dimension
   {
     // micopy_dimension(m_MINCPImpl->m_MincFileDims[m_MINCPImpl->m_DimensionIndices[0]],&apparent_dimension_order[usable_dimensions]);
-    m_MINCPImpl->m_MincApparentDims[usable_dimensions] =
-      m_MINCPImpl->m_MincFileDims[m_MINCPImpl->m_DimensionIndices[0]];
+    m_MINCPImpl->m_MincApparentDims[usableDimensions] = m_MINCPImpl->m_MincFileDims[m_MINCPImpl->m_DimensionIndices[0]];
     // always use positive, vector dimension does not supposed to have notion of positive step size, so leaving as is
     // miset_dimension_apparent_voxel_order(m_MINCPImpl->m_MincApparentDims[usable_dimensions],MI_POSITIVE);
     misize_t _sz;
-    miget_dimension_size(m_MINCPImpl->m_MincApparentDims[usable_dimensions], &_sz);
+    miget_dimension_size(m_MINCPImpl->m_MincApparentDims[usableDimensions], &_sz);
     numberOfComponents = _sz;
-    ++usable_dimensions;
+    ++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[usable_dimensions] =
-      m_MINCPImpl->m_MincFileDims[m_MINCPImpl->m_DimensionIndices[4]];
+    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[usable_dimensions], MI_POSITIVE);
+    miset_dimension_apparent_voxel_order(m_MINCPImpl->m_MincApparentDims[usableDimensions], MI_POSITIVE);
     misize_t _sz;
-    miget_dimension_size(m_MINCPImpl->m_MincApparentDims[usable_dimensions], &_sz);
+    miget_dimension_size(m_MINCPImpl->m_MincApparentDims[usableDimensions], &_sz);
     numberOfComponents = _sz;
-    ++usable_dimensions;
+    ++usableDimensions;
   }
 
   // Set apparent dimension order to the MINC2 api
-  if (miset_apparent_dimension_order(m_MINCPImpl->m_Volume, usable_dimensions, m_MINCPImpl->m_MincApparentDims) < 0)
+  if (miset_apparent_dimension_order(m_MINCPImpl->m_Volume, usableDimensions, m_MINCPImpl->m_MincApparentDims) < 0)
   {
     itkExceptionMacro(" Can't set apparent dimension order!");
   }
 
-  o_origin = dir_cos * origin;
-
-  for (int i = 0; i < spatial_dimension_count; ++i)
-  {
-    this->SetOrigin(i, o_origin[i]);
-  }
-
   miclass_t volume_data_class;
 
   if (miget_data_class(m_MINCPImpl->m_Volume, &volume_data_class) < 0)
@@ -660,6 +678,9 @@ MINCImageIO::ReadImageInformation()
   const std::string classname(GetNameOfClass());
   //  EncapsulateMetaData<std::string>(thisDic,ITK_InputFilterName,
   // classname);
+  // preserve information if the volume was PositiveCoordinateOrientation RAS to PositiveCoordinateOrientation LPS
+  // converted
+  EncapsulateMetaData<bool>(thisDic, "RAStoLPS", m_RAStoLPS);
 
   // store history
   size_t minc_history_length = 0;
@@ -944,22 +965,34 @@ MINCImageIO::WriteImageInformation()
   }
 
   // allocating dimensions
-  vnl_matrix<double> dircosmatrix(nDims, nDims);
-  dircosmatrix.set_identity();
+  vnl_matrix<double> directionCosineMatrix(nDims, nDims);
+  directionCosineMatrix.set_identity();
   vnl_vector<double> origin(nDims);
 
+  // MINC stores direction cosines in PositiveCoordinateOrientation RAS
+  // need to convert to PositiveCoordinateOrientation LPS for ITK
+  vnl_matrix<double> RAS_tofrom_LPS(nDims, nDims);
+  RAS_tofrom_LPS.set_identity();
+  RAS_tofrom_LPS(0, 0) = -1.0;
+  RAS_tofrom_LPS(1, 1) = -1.0;
+
   for (unsigned int i = 0; i < nDims; ++i)
   {
     for (unsigned int j = 0; j < nDims; ++j)
     {
-      dircosmatrix[i][j] = this->GetDirection(i)[j];
+      directionCosineMatrix[i][j] = this->GetDirection(i)[j];
     }
     origin[i] = this->GetOrigin(i);
   }
 
-  const vnl_matrix<double> inverseDirectionCosines{ vnl_matrix_inverse<double>(dircosmatrix).as_matrix() };
+  const vnl_matrix<double> inverseDirectionCosines{ vnl_matrix_inverse<double>(directionCosineMatrix).as_matrix() };
   origin *= inverseDirectionCosines; // transform to minc convention
 
+
+  // Convert ITK direction cosines from PositiveCoordinateOrientation LPS to PositiveCoordinateOrientation RAS
+  if (this->m_RAStoLPS)
+    directionCosineMatrix *= RAS_tofrom_LPS;
+
   for (unsigned int i = 0; i < nDims; ++i)
   {
     const unsigned int j = i + (nComp > 1 ? 1 : 0);
@@ -968,7 +1001,7 @@ MINCImageIO::WriteImageInformation()
     {
       if (k < nDims)
       {
-        dir_cos[k] = dircosmatrix[i][k];
+        dir_cos[k] = directionCosineMatrix[i][k];
       }
       else
       {
@@ -988,38 +1021,32 @@ MINCImageIO::WriteImageInformation()
   {
     case IOComponentEnum::UCHAR:
       m_MINCPImpl->m_Volume_type = MI_TYPE_UBYTE;
-      // m_MINCPImpl->m_Volume_class=MI_CLASS_INT;
       break;
     case IOComponentEnum::CHAR:
       m_MINCPImpl->m_Volume_type = MI_TYPE_BYTE;
-      // m_MINCPImpl->m_Volume_class=MI_CLASS_INT;
       break;
     case IOComponentEnum::USHORT:
       m_MINCPImpl->m_Volume_type = MI_TYPE_USHORT;
-      // m_MINCPImpl->m_Volume_class=MI_CLASS_INT;
       break;
     case IOComponentEnum::SHORT:
       m_MINCPImpl->m_Volume_type = MI_TYPE_SHORT;
-      // m_MINCPImpl->m_Volume_class=MI_CLASS_INT;
       break;
     case IOComponentEnum::UINT:
       m_MINCPImpl->m_Volume_type = MI_TYPE_UINT;
-      // m_MINCPImpl->m_Volume_class=MI_CLASS_INT;
       break;
     case IOComponentEnum::INT:
       m_MINCPImpl->m_Volume_type = MI_TYPE_INT;
-      // m_MINCPImpl->m_Volume_class=MI_CLASS_INT;
       break;
       //     case IOComponentEnum::ULONG://TODO: make sure we are cross-platform here!
       //       volume_data_type=MI_TYPE_ULONG;
       //       break;
       //     case IOComponentEnum::LONG://TODO: make sure we are cross-platform here!
       //       volume_data_type=MI_TYPE_LONG;
       //       break;
-    case IOComponentEnum::FLOAT: // TODO: make sure we are cross-platform here!
+    case IOComponentEnum::FLOAT:
       m_MINCPImpl->m_Volume_type = MI_TYPE_FLOAT;
       break;
-    case IOComponentEnum::DOUBLE: // TODO: make sure we are cross-platform here!
+    case IOComponentEnum::DOUBLE:
       m_MINCPImpl->m_Volume_type = MI_TYPE_DOUBLE;
       break;
     default:
@@ -1059,7 +1086,6 @@ MINCImageIO::WriteImageInformation()
   if (ExposeMetaData<std::string>(thisDic, "dimension_order", dimension_order))
   {
     // the format should be ((+|-)(X|Y|Z|V|T))*
-    // std::cout<<"Restoring original dimension order:"<<dimension_order<<std::endl;
     if (dimension_order.length() == (minc_dimensions * 2))
     {
       dimorder_good = true;
@@ -1296,6 +1322,13 @@ MINCImageIO::WriteImageInformation()
       // TODO: figure out what to do with it
     }
   }
+
+  // preserve information of MINC PositiveCoordinateOrientation RAS to ITK PositiveCoordinateOrientation LPS conversion
+  {
+    int tmp = (int)this->m_RAStoLPS;
+    miset_attr_values(m_MINCPImpl->m_Volume, MI_TYPE_INT, "itk", "RAStoLPS", 1, &tmp);
+  }
+
   mifree_volume_props(hprops);
 }
 

Modules/IO/MINC/src/itkMINCImageIOConfigurePrivate.h.in

@@ -0,0 +1,27 @@
+/*=========================================================================
+ *
+ *  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.
+ *
+ *=========================================================================*/
+#ifndef itkMINCImageIOConfigurePrivate_h
+#define itkMINCImageIOConfigurePrivate_h
+
+// This file is intended to hold preprocessor macros only used internally by
+// IO/MINC module and should not be installed.
+
+// Convert from PositiveCoordinateOrientation RAS to PositiveCoordinateOrientation LPS by default
+#cmakedefine01 ITK_MINC_IO_RAS_TO_LPS
+
+#endif // itkMINCImageIOConfigurePrivate_h