Nifti (RPI orientation) file exported to DICOM - wrong orientation


#1

Hello,

Just FYI.

experiment:
take this Nifti file
https://nifti.nimh.nih.gov/nifti-1/data/avg152T1_LR_nifti.nii.gz/file_view

and export to DICOM

result is horrible wrong.

It happens because GDCM uses Multi Frame Grayscale Byte SC Image IOD Modules.
Internally GDCM takes origin and adds positions one by one, calculating right handed cross product, so far GDCM does nothing wrong, but RPI orientation is not compatible with this approach. IMHO, IO could catch incompatible orientations and apply “change orientation filter” (RPI to RPS in particular case) to re-slice such images to compatible orientation.

Code to reproduce (trivial, read nifti, save dcm, nothing more):
Test

In physical space it looks like:

There are already enough messy multi-frame files, were good not produce more.

Sincerely

P.S. It is not related to discussion about sorting frames in multi-frame images several months ago. Here sorting will not help, wrong positions are written in the image.

Sorry i don have this problem in my app and will not work on the issue.

Edit:
another screenshot (exported image), from Slicer
Slicer


#2

Don’t understand me wrong, ITK handles image orientation ideally, it is not a problem with ITK’s orientation design, it is just required step before exporting to multi-frame DICOM.

If someone cares - link to demo program with required re-orientation step before exporting to multi-frame DICOM . All 48 possible orientations test images (in MHA format with R/L marks) are also inside.

Demo

#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkSpatialOrientation.h>
#include <itkOrientImageFilter.h>
#include <iostream>
#include <string>

int main(int argc, char ** argv)
{
  if (argc < 2)
  {
    std::cout << "File name is required" << std::endl;
    return 0;
  }
  typedef itk::Image<unsigned char, 3> ImageTypeUC;
  typedef itk::ImageFileReader<ImageTypeUC> ReaderType;
  typedef itk::ImageFileWriter<ImageTypeUC> WriterType;
  ReaderType::Pointer reader = ReaderType::New();
  try
  {
    reader->SetFileName(argv[1]);
    reader->Update();
  }
  catch (itk::ExceptionObject & ex)
  {
    std::cout << ex.GetDescription() << std::endl;
    return 1;
  }
  ImageTypeUC::Pointer image = reader->GetOutput();
  itk::SpatialOrientationAdapter adapter;
  unsigned int x = 0;
  const unsigned int orientation = static_cast<unsigned int>(
     adapter.FromDirectionCosines(image->GetDirection()));
  switch (orientation)
  {
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_AIL:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_AIR;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ALS:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ALI;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ARI:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ARS;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ASR:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ASL;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_IAR:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_IAL;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ILA:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ILP;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_IPL:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_IPR;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_IRP:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_IRA;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LAI:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LAS;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LIP:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LIA;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LPS:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LPI;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LSA:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LSP;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PIR:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PIL;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PLI:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PLS;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PRS:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PRI;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PSL:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PSR;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RAS:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RAI;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RIA:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RIP;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RPI:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RPS;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RSP:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RSA;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SAL:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SAR;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SLP:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SLA;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SPR:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SPL;
    break;
  case itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SRA:
    x = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SRP;
    break;
  default:
   break;
  }
  const std::string f = std::string(argv[1]) + std::string(".dcm");
  WriterType::Pointer writer = WriterType::New();
  writer->SetFileName(f.c_str());
  if (x != 0)
  {
    typedef itk::OrientImageFilter<ImageTypeUC,ImageTypeUC>
      OrientImageFilterType;
    OrientImageFilterType::Pointer filter =
      OrientImageFilterType::New();
    try
    {
      filter->SetInput(image);
      filter->UseImageDirectionOn();
      filter->SetDesiredCoordinateOrientation(
        static_cast
          <itk::SpatialOrientation::ValidCoordinateOrientationFlags>
            (x));
      filter->Update();
      writer->SetInput(filter->GetOutput());
    }
    catch (itk::ExceptionObject & ex)
    {
      std::cout << ex.GetDescription() << std::endl;
      return 1;
    }
  }
  else
  {
    writer->SetInput(image);
  }
  try
  {
    writer->Update();
  }
  catch (itk::ExceptionObject & ex)
  {
    std::cout << ex.GetDescription() << std::endl;
    return 1;
  }
  return 0;
}

(Matt McCormick) #3

Hello @mihail.isakov,

Thanks for looking into this and the nice example data and demo code :+1: :candy: .

Do you think we should add the proposed pipeline to itk::GDCMImageIO here?:


#4

Honestly, i didn’t know where to do the step - in itk::GDCMImageIO
void GDCMImageIO::Write(const void *buffer)
works already with void * buffer
and direction is taken from metadata
ExposeMetaData< DoubleMatrixType >( dict, key, directionMatrix );
or already set, probably here is “too late”, i have to learn how IO factory works exactly :nerd_face: :thinking:

Thank you for reply


#5

Looked at image IO.
IMHO, filter could be applied here
itkImageFileWriter.hxx
if IO is GDCM image IO (DCMTK too? never tried) and if image dimension is 3.
Changing of metadata is not required.

For justification (another one):
For DICOM file 6 values from direction matrix will be taken (x and y direction), because 3rd is defined in DICOM as right handed cross product. From GDCM IO:

  // Do the direction now:
  // if the meta dictionary contains the tag "0020 0037", use it
  const bool hasIOP = ExposeMetaData<std::string>(dict, "0020|0037",tempString);
  if (hasIOP)
  {
    double directions[6];
    sscanf(tempString.c_str(), "%lf\\%lf\\%lf\\%lf\\%lf\\%lf", &(directions[0]), &(directions[1]), &(directions[2]),&(directions[3]),&(directions[4]),&(directions[5]));
    image.SetDirectionCosines(0, directions[0]);
    image.SetDirectionCosines(1, directions[1]);
    image.SetDirectionCosines(2, directions[2]);
    image.SetDirectionCosines(3, directions[3]);
    image.SetDirectionCosines(4, directions[4]);
    image.SetDirectionCosines(5, directions[5]);
  }
  else
  {
    image.SetDirectionCosines(0, m_Direction[0][0]);
    image.SetDirectionCosines(1, m_Direction[0][1]);
    if ( m_Direction.size() == 3 )
      {
      image.SetDirectionCosines(2, m_Direction[0][2]);
      }
    else
      {
      image.SetDirectionCosines(2, 0);
      }
    image.SetDirectionCosines(3, m_Direction[1][0]);
    image.SetDirectionCosines(4, m_Direction[1][1]);
    if ( m_Direction.size() == 3 )
      {
      image.SetDirectionCosines(5, m_Direction[1][2]);
      }
    else
      {
      image.SetDirectionCosines(5, 0);
      }
  }

E.g. RAI and RAS:

1 0 0
0 1 0
0 0 1

and

1 0 0
0 1 0
0 0 -1

For both orientations, the same Image Orientation Patient “1, 0, 0, 0, 1, 0” will be assumed and for multi-frame image GDCM will calculate wrong slice origins for RAS image.

Sincerely


(Matt McCormick) #6

It is tricky since itkImageFileWriter.hxx cannot depend on itkGDCMImageIO.hitkGDCMImageIO.h may not be available if the ITKIOGDCM module is not enabled, and the ITKIOImageBase module cannot depend on on the ITKIOGDCM module (the reverse dependency exists).

In one approach, if DCMTK behaves the same, we could apply the fix for output files with the .dcm.

Or, we could apply the fix in itkGDCMImageIO.cxx, even though it is more difficult since the image has been decomposed.


#7

Oh, yes. Important to know, thank you.
I’ve just tried to add

std::cout << m_ImageIO->GetNameOfClass() << std::endl;

at this point, line 174, before line with
auto * nonConstInput = const_cast< InputImageType * >( input );

seems m_ImageIO knows name of class, without including new headers.
Probably it is enough, i’ll try and report.

Sincerely


#8

@matt.mccormick you are right, of course. ImageFileWriter is bad place, the problem is - some things happen in ::Write(), some in ::GenerateData(), always accessing image like
const InputImageType *input = this->GetInput();
and buffered region, streaming, pipeline, many possible issues. It is possible, but will require too many changes, not worth… Probably possible ways are either do the job in GDCM IO class, as you have suggested, also tricky, or just give a warning (or error). I shall play again little later… Thank you