NiftiImageIO: Round Trip issue with Image Orientation (Patient)

[malaterre]

Description

Placeholder for possibly a longer bug report. But at least store the information for reproducing the issue.

Steps to Reproduce

So basically the steps are:

1. cp gdcm/Testing/Data/*SIEM*FileSeq* /tmp/dir/
2. gdcmanon --dumb --replace 0020,0037="0.9997586873\0.02191020543\0.001584328775\0.02196741222\-0.9971551496\-0.0721043717" /tmp/dir /tmp/dir2
3. bin/DicomSeriesReadImageWrite2 /tmp/dir2 foo.nii
4. bin/ImageReadWrite foo.nii nii.dcm

Expected behavior

Minimal approximation of the matrix orientation in between the different representations.

$ gdcmdump nii.dcm | grep Orie
        (0020,0037) ?? (DS) [.999758687284262\.021910205429893\.001584328774992\.021967412220091\-.99715514960414\-.07210437170029 ]         # 102,6 Image Orientation (Patient)

This goal is achievable since the MetaImage seems to be giving proper results:

1. bin/DicomSeriesReadImageWrite2 /tmp/dir2 foo.mha
2. bin/ImageReadWrite foo.mha mha.dcm

Leads to:

$ gdcmdump mha.dcm | grep Orie
        (0020,0037) ?? (DS) [.999758687284262\.021910205429893\.001584328774992\.021967412220091\-.99715514960414\-.07210437170029 ]         # 102,6 Image Orientation (Patient)

Actual behavior

Leads to:

$ gdcmdump nii.dcm | grep Orie
        (0020,0037) ?? (DS) [.9997590014762\.021938811333211\7.921645009e-4\.021938811254904\-.99715514907211\-.07211308646627 ]         # 98,6 Image Orientation (Patient)

Reproducibility

Versions

$ git rev-parse --short=12 HEAD
b089be7fdd16

Environment

Linux, but I am guessing this will be irrelevant.

Additional Information

$ git diff
diff --git a/Examples/IO/DicomSeriesReadImageWrite2.cxx b/Examples/IO/DicomSeriesReadImageWrite2.cxx
index be15736..d435c10 100644
--- a/Examples/IO/DicomSeriesReadImageWrite2.cxx
+++ b/Examples/IO/DicomSeriesReadImageWrite2.cxx
@@ -74,7 +74,7 @@ int main( int argc, char* argv[] )
 // Software Guide : EndLatex
 
 // Software Guide : BeginCodeSnippet
-  using PixelType = signed short;
+  using PixelType = unsigned short;
   constexpr unsigned int Dimension = 3;
 
   using ImageType = itk::Image< PixelType, Dimension >;
diff --git a/Examples/IO/ImageReadWrite.cxx b/Examples/IO/ImageReadWrite.cxx
index 0ccc48a..1eb8e71 100644
--- a/Examples/IO/ImageReadWrite.cxx
+++ b/Examples/IO/ImageReadWrite.cxx
@@ -85,8 +85,8 @@ int main( int argc, char ** argv )
   //  Software Guide : EndLatex
 
   // Software Guide : BeginCodeSnippet
-  using PixelType = short;
-  constexpr unsigned int Dimension = 2;
+  using PixelType = unsigned short;
+  constexpr unsigned int Dimension = 3;
   using ImageType = itk::Image< PixelType, Dimension >;
   // Software Guide : EndCodeSnippet
 

[malaterre]

For reference:

$ nifti_tool -disp_hdr -infiles foo.nii

header file 'foo.nii', num_fields = 43

all fields:
  name                offset  nvals  values
  ------------------- ------  -----  ------
  sizeof_hdr             0      1    348
  data_type              4     10    
  db_name               14     18    
  extents               32      1    0
  session_error         36      1    0
  regular               38      1    r
  dim_info              39      1    0
  dim                   40      8    3 128 128 4 1 1 1 1
  intent_p1             56      1    0.0
  intent_p2             60      1    0.0
  intent_p3             64      1    0.0
  intent_code           68      1    0
  datatype              70      1    512
  bitpix                72      1    16
  slice_start           74      1    0
  pixdim                76      8    1.0 1.875 1.875 5.5 0.0 0.0 0.0 0.0
  vox_offset           108      1    352.0
  scl_slope            112      1    1.0
  scl_inter            116      1    0.0
  slice_end            120      1    0
  slice_code           122      1    0
  xyzt_units           123      1    2
  cal_max              124      1    0.0
  cal_min              128      1    0.0
  slice_duration       132      1    0.0
  toffset              136      1    0.0
  glmax                140      1    0
  glmin                144      1    0
  descrip              148     80    
  aux_file             228     24    
  qform_code           252      1    1
  sform_code           254      1    0
  quatern_b            256      1    0.010977
  quatern_c            260      1    -0.999288
  quatern_d            264      1    0.036082
  qoffset_x            268      1    120.0
  qoffset_y            272      1    120.0
  qoffset_z            276      1    -30.0383
  srow_x               280      4    0.0 0.0 0.0 0.0
  srow_y               296      4    0.0 0.0 0.0 0.0
  srow_z               312      4    0.0 0.0 0.0 0.0
  intent_name          328     16    
  magic                344      4    n+1

[malaterre]

For reference I tried converting the DICOM series to NifTI using vtk-dicom, here is what I found out:

$ dicomtonifti --no-column-reordering --no-row-reordering -o bla.nii /tmp/dir2/SIEMENS_MAGNETOM-12-MONO2-FileSeq*

Which gives me:

$ nifti_tool -diff_hdr -infiles foo.nii bla.nii 
  name                offset  nvals  values
  ------------------- ------  -----  ------
  regular               38      1    r
  regular               38      1    
  pixdim                76      8    1.0 1.875 1.875 5.5 0.0 0.0 0.0 0.0
  pixdim                76      8    1.0 1.875 1.875 5.485677 1.0 1.0 1.0 1.0
  scl_slope            112      1    1.0
  scl_slope            112      1    0.0
  xyzt_units           123      1    2
  xyzt_units           123      1    18
  cal_max              124      1    0.0
  cal_max              124      1    871.0
  descrip              148     80    
  descrip              148     80    01/Mar/2001 PROT-DP/DIFF-0-500-1000-TR 1
  sform_code           254      1    0
  sform_code           254      1    2
  quatern_c            260      1    -0.999288
  quatern_c            260      1    -0.999288
  srow_x               280      4    0.0 0.0 0.0 0.0
  srow_x               280      4    -1.874547 -0.041189 -0.0 120.0
  srow_y               296      4    0.0 0.0 0.0 0.0
  srow_y               296      4    -0.041082 1.869666 -0.395637 120.0
  srow_z               312      4    0.0 0.0 0.0 0.0
  srow_z               312      4    0.002971 -0.135196 -5.471392 -30.0383

With this input NIfTI file here is what I get:

$ bin/ImageReadWrite bla.nii bla.dcm
$ gdcmdump bla.dcm | grep Ori
        (0020,0037) ?? (DS) [.999758623648593\.021907443474456\.001660831323115\.021970175225087\-.99715514931604\-.07210353384595 ]         # 102,6 Image Orientation (Patient)

[dzenanz]

Self-contained code to reproduce contributed by @issakomi:

#include "itkImage.h"
#include "itkImageFileWriter.h"
#include "itkImageFileReader.h"

int main(int argc, char ** argv)
{
  typedef itk::Image<float, 3> ImageType; 
  typedef itk::ImageFileWriter<ImageType> WriterType;
  typedef itk::ImageFileReader<ImageType> ReaderType;
  {
    ImageType::Pointer image = ImageType::New();
    ImageType::IndexType index;
    ImageType::SizeType size;
    ImageType::RegionType region;
    ImageType::SpacingType spacing;
    ImageType::PointType origin;
    ImageType::DirectionType direction;
    index.Fill(0);
    size[0] = 192;
    size[1] = 256;
    size[2] = 192;
    spacing[0] = 0.976562;
    spacing[1] = 0.976562;
    spacing[2] = 1;
    origin[0] = -100.04;
    origin[1] = -158.669;
    origin[2] = -101.794;
    direction(0,0) = 0.999998;
    direction(0,1) = 0.00174533;
    direction(0,2) = 0;
    direction(1,0) = -0.00174533;
    direction(1,1) = 0.999998;
    direction(1,2) = 0;
    direction(2,0) = 0;
    direction(2,1) = 0;
    direction(2,2) = 1;
    region.SetIndex(index);
    region.SetSize(size);
    try
    {
      image->SetRegions(region);
      image->Allocate();
      image->SetOrigin(origin);
      image->SetSpacing(spacing);
      image->SetDirection(direction);
    }
    catch (itk::ExceptionObject & ex)
    {
      std::cout << ex.GetDescription() << std::endl;
      return 1;
    }
    std::cout << "*** Origin " << image->GetOrigin() << std::endl;
    std::cout << "Direction\n" << image->GetDirection() << std::endl;
    WriterType::Pointer writer = WriterType::New();
    writer->SetInput(image);
    writer->SetFileName("output.nii.gz"); // MHA is OK
    try
    {
      writer->Write();
    }
    catch (itk::ExceptionObject & ex)
    {
      std::cout << ex.GetDescription() << std::endl;
      return 1;
    }
  }

  {
    ImageType::Pointer reopened;
    ReaderType::Pointer reader = ReaderType::New();
    reader->SetFileName("output.nii.gz"); // MHA is OK
    try
    {
      reader->Update();
      reopened = reader->GetOutput();
    }
    catch (itk::ExceptionObject & ex)
    {
      std::cout << ex.GetDescription() << std::endl;
      return 1;
    }
    std::cout << "*** Re-opened\nOrigin " << reopened->GetOrigin() << std::endl;
    std::cout << "Direction\n" << reopened->GetDirection() << std::endl;
  }
  return 0;
}

example output:

*** Origin [-100.04, -158.669, -101.794]
Direction
0.999998 0.00174533 0
-0.00174533 0.999998 0
0 0 1

*** Re-opened
Origin [-100.04, -158.669, -101.794]
Direction
0.999998 0.00182698 0
-0.00182698 0.999998 0
0 0 1

[issakomi]

Some points mitigate the issue.
Closer look at vectors provided in the very first post above

x              [   0.999759, 0.0219102, 0.00158433]
y              [  0.0219674, -0.997155, -0.0721044]

z = cross(x,y) = [3.02631e-13, 0.0721218, -0.997396 ]
dot(x,y) = 2.39187e-13
dot(x,z) = -2.1684e-19
dot(y,z) = 0

not a ideal orthogonal matrix, the error is numerically very small, but it is here, scalar product must be zero.

The data in my test above, taken from another file also has little numerical issue

dot(x,y) -9.99998e-10
dot(x,z) 0
dot(y,z) 0

If we see closer it is visible

Direction
0.999998 0.00174533 0
-0.00174533 0.999998 0
0 0 1

See at sign in the 2nd row

If the data were ideal, for example

0.999998 -0.00174533 0
-0.00174533 0.999998 0
0 0 1

or

-0.999998 0.00174533 0
-0.00174533 0.999998 0
0 0 1

Nifti IO has less problem

*** Re-opened
Origin [-100.04, -158.669, -101.794]
Direction
1 0 0
0 1 0
0 0 1

Still some loss in precision, too big, but better

[issakomi]

It is better to use updated example for reproducing

#include "itkImage.h"
#include "itkImageFileWriter.h"
#include "itkImageFileReader.h"

int main(int argc, char ** argv)
{
  typedef itk::Image<float, 3> ImageType; 
  typedef itk::ImageFileWriter<ImageType> WriterType;
  typedef itk::ImageFileReader<ImageType> ReaderType;
  {
    ImageType::Pointer image = ImageType::New();
    ImageType::IndexType index;
    ImageType::SizeType size;
    ImageType::RegionType region;
    ImageType::SpacingType spacing;
    ImageType::PointType origin;
    ImageType::DirectionType direction;
    index.Fill(0);
    size[0] = 192;
    size[1] = 256;
    size[2] = 192;
    spacing[0] = 0.976562;
    spacing[1] = 0.976562;
    spacing[2] = 1;
    origin[0] = 0;
    origin[1] = 0;
    origin[2] = 0;
    direction(0,0) = 0.999998;
    direction(0,1) = -0.00174533;
    direction(0,2) = 0;
    direction(1,0) = -0.00174533;
    direction(1,1) = 0.999998;
    direction(1,2) = 0;
    direction(2,0) = 0;
    direction(2,1) = 0;
    direction(2,2) = 1;
    region.SetIndex(index);
    region.SetSize(size);
    try
    {
      image->SetRegions(region);
      image->Allocate();
      image->SetOrigin(origin);
      image->SetSpacing(spacing);
      image->SetDirection(direction);
    }
    catch (itk::ExceptionObject & ex)
    {
      std::cout << ex.GetDescription() << std::endl;
      return 1;
    }
    std::cout << "*** Direction\n"
      << image->GetDirection() << std::endl;
    WriterType::Pointer writer = WriterType::New();
    writer->SetInput(image);
    writer->SetFileName("output.nii.gz"); // MHA is OK
    try
    {
      writer->Write();
    }
    catch (itk::ExceptionObject & ex)
    {
      std::cout << ex.GetDescription() << std::endl;
      return 1;
    }
  }

  {
    ImageType::Pointer reopened;
    ReaderType::Pointer reader = ReaderType::New();
    reader->SetFileName("output.nii.gz"); // MHA is OK
    try
    {
      reader->Update();
      reopened = reader->GetOutput();
    }
    catch (itk::ExceptionObject & ex)
    {
      std::cout << ex.GetDescription() << std::endl;
      return 1;
    }
    std::cout << "*** Re-opened\nDirection\n"
      << reopened->GetDirection() << std::endl;
  }
  return 0;
}

Output is

*** Direction
0.999998 -0.00174533 0
-0.00174533 0.999998 0
0 0 1

*** Re-opened
Direction
1 0 0
0 1 0
0 0 1

Edit: i looked closer and as already mentioned there are numerical limitations, at least as long qto_xyz matrix is used, it goes thru quaternions math, so there is no bug, but limitation. May be there are alternative ways in Nifti standard, i don't know.

[seanm]

CC @seanm @vfonov

[hjmjohnson]

Nifti stores the direction cosine as a quaternion, not an orthogonal matrix. This is a fundamental issue with the storage mechanism for nifti. There are no alternatives in the Nifti standard to address this issue :(.

[malaterre]

@hjmjohnson As a gentle reminder, this is a regression, so anything that revert the changeset would be fine for me.

AFAIK, I also doubt NIfTI only way of storing orientation is through quaternions (I believe this is an implementation detail of ITK to pick quaternion over the matrix one).

[stale]

This issue has been automatically marked as stale because it has not had recent activity. Thank you for your contributions.