OrientImageFilter changes the origin of the image

sarthakpati · August 30, 2019, 1:48pm

Hi,

I am trying to change the orientation of an image using the OrientImageFilter but I am running into a very weird issue. Once I convert the original image to RAI and then convert it back to it’s original orientation (in this case, it is LPI), the origin of the final “converted” LPI image and the original image is not the same:

Original Image: 86.92, 90.03, -101.8
Original_RAI_LPI: 86.94, 90.01, -101.8

Though the actual difference is very small, in my subsequent computation pipeline, it causes issues since the original image and the reconverted images are no longer in the same physical space.

I am using the following chunk of code for this conversion:

template< class TImageType = ImageTypeFloat3D >
std::pair< std::string, typename TImageType::Pointer > GetImageOrientation(const typename TImageType::Pointer inputImage, const std::string &desiredOrientation = "RAI")
{
  if (TImageType::ImageDimension != 3)
  {
    std::cerr << "This function is only defined for 3D images.\n";
    exit(EXIT_FAILURE);
  }
  auto orientFilter = itk::OrientImageFilter< TImageType, TImageType >::New();
  orientFilter->SetInput(inputImage);
  orientFilter->UseImageDirectionOn();
  orientFilter->SetDirectionTolerance(0);
  orientFilter->SetCoordinateTolerance(0);

  auto desiredOrientation_wrap = desiredOrientation;
  std::transform(desiredOrientation_wrap.begin(), desiredOrientation_wrap.end(), desiredOrientation_wrap.begin(), ::toupper);
  
  std::map< std::string, itk::SpatialOrientation::ValidCoordinateOrientationFlags > orientationMap;
  orientationMap["Axial"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RAI;
  orientationMap["Coronal"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RSA;
  orientationMap["Sagittal"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ASL;
  orientationMap["RIP"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RIP;
  orientationMap["LIP"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LIP;
  orientationMap["RSP"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RSP;
  orientationMap["LSP"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LSP;
  orientationMap["RIA"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RIA;
  orientationMap["LIA"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LIA;
  orientationMap["RSA"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RSA;
  orientationMap["LSA"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LSA;
  orientationMap["IRP"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_IRP;
  orientationMap["ILP"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ILP;
  orientationMap["SRP"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SRP;
  orientationMap["SLP"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SLP;
  orientationMap["IRA"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_IRA;
  orientationMap["ILA"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ILA;
  orientationMap["SRA"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SRA;
  orientationMap["SLA"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SLA;
  orientationMap["RPI"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RPI;
  orientationMap["LPI"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LPI;
  orientationMap["RAI"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RAI;
  orientationMap["LAI"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LAI;
  orientationMap["RPS"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RPS;
  orientationMap["LPS"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LPS;
  orientationMap["RAS"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RAS;
  orientationMap["LAS"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LAS;
  orientationMap["PRI"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PRI;
  orientationMap["PLI"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PLI;
  orientationMap["ARI"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ARI;
  orientationMap["ALI"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ALI;
  orientationMap["PRS"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PRS;
  orientationMap["PLS"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PLS;
  orientationMap["ARS"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ARS;
  orientationMap["ALS"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ALS;
  orientationMap["IPR"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_IPR;
  orientationMap["SPR"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SPR;
  orientationMap["IAR"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_IAR;
  orientationMap["SAR"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SAR;
  orientationMap["IPL"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_IPL;
  orientationMap["SPL"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SPL;
  orientationMap["IAL"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_IAL;
  orientationMap["SAL"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SAL;
  orientationMap["PIR"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PIR;
  orientationMap["PSR"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PSR;
  orientationMap["AIR"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_AIR;
  orientationMap["ASR"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ASR;
  orientationMap["PIL"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PIL;
  orientationMap["PSL"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PSL;
  orientationMap["AIL"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_AIL;
  orientationMap["ASL"] = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ASL;

  // set the desired orientation and update
  orientFilter->SetDesiredCoordinateOrientation(orientationMap[desiredOrientation_wrap]);
  orientFilter->Update();
  auto outputImage = orientFilter->GetOutput();

  std::string originalOrientation;

  for (auto it = orientationMap.begin(); it != orientationMap.end(); ++it)
  {
    if (it->second == orientFilter->GetGivenCoordinateOrientation())
    {
      originalOrientation = it->first;
    }
  }
  if (originalOrientation.empty())
  {
    originalOrientation = "Unknown";
  }

  return std::make_pair(originalOrientation, outputImage);
}

All the data (original + converted images) can be found here: https://1drv.ms/u/s!AgvRZZtXCbbOnIwGa6dEHvJ75jvSOQ?e=xf5TiS

Any help on this will be much appreciated.

Thanks,
Sarthak

EDIT: I tried resampling the image to an isotropic resolution of 1,1,1 before applying the orientation fix and that made the problem worse.

EDIT2: Run-able example can be found here: https://github.com/sarthakpati/ITKOrientFilterExample

dzenanz · August 30, 2019, 9:00pm

Do you mind turning this into a runnable example? E.g. read an image, reorient it, write it.

mihail.isakov · August 30, 2019, 9:36pm

so far i have such functions too, to apply change orientation filter and save and could not reproduce.
Results are correct, IMHO.

original
LPI  Origin: [86.9175, 90.0284, -101.794]
filter result
LPI -> RAI  Origin: [-100.04, -158.669, -101.794]
RAI -> LPI  Origin: [86.9175, 90.0284, -101.794]

saved to nii.gz and re-loaded
LPI -> RAI  Origin: [-100.04, -158.669, -101.794]
RAI -> LPI  Origin: [86.9175, 90.0284, -101.794]

sarthakpati · September 3, 2019, 12:08pm

Thanks for the response, @dzenanz. Here you go: https://github.com/sarthakpati/ITKOrientFilterExample

sarthakpati · September 3, 2019, 12:09pm

Could you please give me a link to your function(s)? I may be grossly mistaken in what I am trying to do in the code highlighted (I have updated the original post with a working example, if you are interested).

dzenanz · September 3, 2019, 3:44pm

I can reproduce it:

Original Orientation: LPI
Started verification.
Original Image Properties:
        Origin: [86.9175, 90.0284, -101.794]
Oriented Image [RAI] Properties:
        Origin: [-100.04, -158.669, -101.794]
Oriented Image [ORIGINAL] Properties:
        Origin: [86.9378, 90.0131, -101.794]
Finished successfully.

but it doesn’t occur with some other image I had lying around:

Original Orientation: RAI
Started verification.
Original Image Properties:
        Origin: [-127.5, -127.5, -127.5]
Oriented Image [RAI] Properties:
        Origin: [-127.5, -127.5, -127.5]
Oriented Image [ORIGINAL] Properties:
        Origin: [-127.5, -127.5, -127.5]
Finished successfully.

so it seems to be related to something specific in your image. I am investigating.

sarthakpati · September 3, 2019, 4:53pm

Is your image isotropic? Perhaps that might be a source of difference. In any case, thank you for looking into it! I am looking forward to hearing from you.

blowekamp · September 3, 2019, 5:00pm

I was just poking at this issue.

That is the issue! Look at this code in PermuteAxesImageFilter:

github.com

InsightSoftwareConsortium/ITK/blob/d2e511437d3d9480ae1ab7ea8beb210696b6ceee/Modules/Filtering/ImageGrid/include/itkPermuteAxesImageFilter.hxx#L149-L162


      
          for (j = 0; j < ImageDimension; j++)
          {
            // origin does not change by a Permute.  But spacing, directions,
            // size and start index do.
            outputOrigin[j] = inputOrigin[j];
          
            outputSpacing[j] = inputSpacing[m_Order[j]];
            outputSize[j] = inputSize[m_Order[j]];
            outputStartIndex[j] = inputStartIndex[m_Order[j]];
            for (i = 0; i < ImageDimension; i++)
            {
              outputDirection[i][j] = inputDirection[i][m_Order[j]];
            }
          }

For the case of anisotropic voxels the origin should change! While the location for the corner of the voxel should remains the same, the location of the center is suppose to change for the PermuteImageFiler.

sarthakpati · September 3, 2019, 5:12pm

Thanks for the response!

I understand that the origins should change for anisotropic images but when I change the image back to the original orientation, should it not be preserved? From the example I have shown (and verified by @dzenanz), there is a minor difference in the final re-oriented and original images.

blowekamp · September 3, 2019, 5:26pm

The above code is from the PermuteImageFilter which is used by the OrientImageFilter. I believe it is the cause of the bug you are seeing as it does not correctly compute origin after the permutation of the order of the axes.

sarthakpati · September 3, 2019, 5:33pm

Ah, I understand. I got thrown off by the fact that the origin changed for me.

mihail.isakov · September 3, 2019, 5:40pm

Edited: typo

FYI, my test, original image, (BTW it is signed short image, but opened as float to reproduce above test:

r@deb:~/tmp/test038/build$ ./test038 original.nii.gz 
original LPI [86.9175, 90.0284, -101.794]
filter output LPI->RAI [-100.04, -158.669, -101.794]
filter output RAI->LPI [86.9175, 90.0284, -101.794]
re-opened LPI->RAI [-100.04, -158.669, -101.794]
re-opened RAI->LPI [86.9175, 90.0284, -101.794]

Here is code

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

template<typename T> bool orient_filter(
  const typename T::Pointer & image,
  typename T::Pointer & out_image,
  int f)
{
  if (image.IsNull()) return false;
  typedef itk::OrientImageFilter<T,T> OrientImageFilterType;
  typename OrientImageFilterType::Pointer filter = OrientImageFilterType::New();
  try
  {
    filter->SetInput(image);
    filter->UseImageDirectionOn();
    filter->SetDesiredCoordinateOrientation(
      static_cast<itk::SpatialOrientation::ValidCoordinateOrientationFlags>(f));
    filter->Update();
    out_image = filter->GetOutput();
  }
  catch (itk::ExceptionObject & ex)
  {
    std::cout << ex.GetDescription() << std::endl;
    return false;
  }
  if (out_image.IsNotNull()) out_image->DisconnectPipeline();
  else return false;
  return true;
}

int main(int argc, char ** argv)
{
  typedef itk::Image<float, 3> ImageType; 
  typedef itk::ImageFileReader<ImageType> ReaderType;
  typedef itk::ImageFileWriter<ImageType> WriterType;
  ImageType::Pointer inputLPI;
  ImageType::Pointer outputRAI;
  ImageType::Pointer outputLPI;
  ImageType::Pointer reopenedRAI;
  ImageType::Pointer reopenedLPI;
  ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName(argv[1]);
  try
  {
    reader->Update();
    inputLPI = reader->GetOutput();
  }
  catch (itk::ExceptionObject & ex)
  {
    std::cout << ex.GetDescription() << std::endl;
    return 1;
  }
  std::cout << "original LPI " << inputLPI->GetOrigin()
    << std::endl;
  orient_filter<ImageType>(
    inputLPI,
    outputRAI,
    (int)itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RAI);
  std::cout << "filter output LPI->RAI " << outputRAI->GetOrigin()
    << std::endl;
  orient_filter<ImageType>(
    outputRAI,
    outputLPI,
    (int)itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LPI);
  std::cout << "filter output RAI->LPI " << outputLPI->GetOrigin()
    << std::endl;
  WriterType::Pointer writer1 = WriterType::New();
  WriterType::Pointer writer2 = WriterType::New();
  writer1->SetInput(outputRAI);
  writer1->SetFileName("outputRAI.nii.gz");
  writer2->SetInput(outputLPI);
  writer2->SetFileName("outputLPI.nii.gz");
  try
  {
    writer1->Update();
    writer2->Update();
  }
  catch (itk::ExceptionObject & ex)
  {
    std::cout << ex.GetDescription() << std::endl;
    return 1;
  }
  ReaderType::Pointer reader1 = ReaderType::New();
  ReaderType::Pointer reader2 = ReaderType::New();
  reader1->SetFileName("outputRAI.nii.gz");
  reader2->SetFileName("outputLPI.nii.gz");
  try
  {
    reader1->Update();
    reader2->Update();
    reopenedRAI = reader1->GetOutput();
    reopenedLPI = reader2->GetOutput();
  }
  catch (itk::ExceptionObject & ex)
  {
    std::cout << ex.GetDescription() << std::endl;
    return 1;
  }
  std::cout << "re-opened LPI->RAI " << reopenedRAI->GetOrigin()
    << std::endl;
  std::cout << "re-opened RAI->LPI n" << reopenedLPI->GetOrigin()
    << std::endl;
  return 0;
}

nnn

dzenanz · September 3, 2019, 5:56pm

The other image I have is isotropic (spacing 1.0).

I created an issue to track this. Do we have a course of action regarding this? Are @blowekamp or @mihail.isakov working on a fix already? If yes, please assign the issue to yourself.

mihail.isakov · September 3, 2019, 6:03pm

I still don’t see the issue, change orientation filter works well for me, s. test above, not sure, may i be oversee something, of course, example is a little too long for me dig into just now, sorry.

sarthakpati · September 3, 2019, 6:48pm

After trying to replicate what @mihail.isakov has done, I tried doing the re-conversion of the RAI image from both the itk::Image Pointer and the saved image on disk [ref] and it turns out, there are differences when the file is written from the disk (see output):

Original Orientation: LPI
Started verification.
Original Image Properties:
        Origin: [86.9175, 90.0284, -101.794]
Oriented Image [RAI] Properties:
        Origin: [-100.04, -158.669, -101.794]
Re-Oriented Image [LPI] - from Raw RAI Properties:
        Origin: [86.9175, 90.0284, -101.794]
Re-Oriented Image [LPI] - from Written RAI Properties:
        Origin: [86.9378, 90.0131, -101.794]

blowekamp · September 3, 2019, 6:56pm

Here is my untested fix:

If someone else case write a test, it would be much appreciated.

sarthakpati · September 3, 2019, 11:43pm

I pulled these fixes but the problem persists:

Original Orientation: LPI
Started verification.
Original Image Properties:
        Origin: [86.9175, 90.0284, -101.794]
Oriented Image [RAI] Properties:
        Origin: [-100.04, -158.669, -101.794]
Re-Oriented Image [LPI] - from Raw RAI Properties:
        Origin: [86.9175, 90.0284, -101.794]
Re-Oriented Image [LPI] - from Written RAI Properties:
        Origin: [86.9378, 90.0131, -101.794]

sarthakpati · September 4, 2019, 11:52am

I believe the issue is coming from the way the writer is getting called (though I have no idea why). I am trying to find the direction cosines of the same image here, one right after the orientation function and another after I have written and read it back from disk:

Oriented Image [RAI] Raw Properties:
        Origin: [-100.04, -158.669, -101.794]
        Direction Cosines:
0.999998 0.00174533 0
-0.00174533 0.999998 0
0 0 1

Oriented Image [RAI] FromFile Properties:
        Origin: [-100.04, -158.669, -101.794]
        Direction Cosines:
0.999998 0.00182698 0
-0.00182698 0.999998 0
0 0 1

mihail.isakov · September 4, 2019, 8:40pm

Thank you for looking into, interesting. I will do more tests with orientation filter, probably i’ll modify my example to apply a number of random re-orientations, e.g. 20 or so and look at origin, spacing, direction, without writing. I hope the results will be OK. I’ll post.
I use Nifti IO as ITK offers, no time to dig into, may be there is something with quaternions math, just guessing.

blowekamp · September 4, 2019, 8:43pm

After writing some tests to verify the geometry for the PermuteImageFilter, I determined that the existing computation of the origin is correct.