BUG: Fix case where DICOM series order has negative spacing

If the spacing between slices is negative, then direction cosine
matrix was not correct with reguards to the sign.

This can be easily created with the "ReverseOrder" flag with the
series reader and DICOM data. With the correction, a DICOM series
mantains its correct physical location when this flag is enabled.

Author: blowekamp

Modules/IO/GDCM/test/CMakeLists.txt

@@ -448,3 +448,28 @@ addcompliancetest(raw-YBR_FULL_422)
 addcompliancetest(RLE-RGB)
 addcompliancetest(HTJ2K-YBR_ICT Lily_full.mha)
 addcompliancetest(HTJ2K-YBR_RCT Lily_full.mha)
+
+set(ITKGDCMImageIOGTests itkGDCMImageIOGTest.cxx)
+creategoogletestdriver(ITKGDCMImageIO "${ITKIOGDCM-Test_LIBRARIES}" "${ITKGDCMImageIOGTests}")
+
+ExternalData_Expand_Arguments(
+  ITKData
+  _DICOM_SERIES_INPUT
+  "DATA{${ITK_DATA_ROOT}/Input/DicomSeries/,REGEX:Image[0-9]+.dcm}"
+)
+target_compile_definitions(
+  ITKGDCMImageIOGTestDriver
+  PRIVATE
+    "DICOM_SERIES_INPUT=${_DICOM_SERIES_INPUT}"
+  PRIVATE
+    "ITK_TEST_OUTPUT_DIR=${ITK_TEST_OUTPUT_DIR}"
+)
+
+set_property(
+  TARGET
+    ITKGDCMImageIOGTestDriver
+  APPEND
+  PROPERTY
+    DEPENDS
+      ITKData
+)

Modules/IO/GDCM/test/itkGDCMImageIOGTest.cxx

@@ -0,0 +1,333 @@
+/*=========================================================================
+ *
+ *  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.
+ *
+ *=========================================================================*/
+
+#include "itkImage.h"
+#include "itkImageFileWriter.h"
+#include "itkGDCMImageIO.h"
+#include "itkGDCMSeriesFileNames.h"
+#include "itkMetaDataObject.h"
+#include "itkGTest.h"
+#include "itksys/SystemTools.hxx"
+#include "itkImageSeriesReader.h"
+#include <string>
+#include <vector>
+
+
+#define _STRING(s) #s
+#define TOSTRING(s) std::string(_STRING(s))
+
+namespace
+{
+
+struct ITKGDCMImageIO : public ::testing::Test
+{
+  void
+  SetUp() override
+  {
+    itksys::SystemTools::MakeDirectory(m_TempDir);
+  }
+
+  void
+  TearDown() override
+  {
+    itksys::SystemTools::RemoveADirectory(m_TempDir);
+  }
+
+  const std::string m_TempDir{ TOSTRING(ITK_TEST_OUTPUT_DIR) + "/ITKGDCMImageIO" };
+  const std::string m_DicomSeriesInput{ TOSTRING(DICOM_SERIES_INPUT) };
+};
+
+class ITKGDCMSeriesTestData : public ITKGDCMImageIO
+{
+public:
+  using PixelType = uint16_t;
+  static constexpr unsigned int Dimension = 2;
+  using ImageType = itk::Image<PixelType, Dimension>;
+  using WriterType = itk::ImageFileWriter<ImageType>;
+
+  void
+  SetUp() override
+  {
+    itksys::SystemTools::MakeDirectory(m_TempDir);
+    CreateTestDicomSeries();
+  }
+
+  void
+  TearDown() override
+  {
+    itksys::SystemTools::RemoveADirectory(m_TempDir);
+  }
+
+protected:
+  const std::string        m_TempDir{ TOSTRING(ITK_TEST_OUTPUT_DIR) + "/ITKGDCMSeriesTestData" };
+  std::vector<std::string> m_DicomFiles;
+
+
+private:
+  void
+  CreateTestDicomSeries()
+  {
+    // DICOM meta-data values from the report
+    const std::vector<std::string> positions = {
+      "-216.500\\-216.500\\70.000",   // slice 1 (top)
+      "-216.500\\-216.500\\-187.500", // slice 2 (middle)
+      "-216.500\\-216.500\\-445.000"  // slice 3 (bottom)
+    };
+
+    const std::string orientation = "1.000000\\0.000000\\0.000000\\0.000000\\1.000000\\0.000000";
+
+    // Create a 2x2 image for each slice (equivalent to 3D [2,2,3] sliced)
+    ImageType::SizeType size;
+    size[0] = 2;
+    size[1] = 2;
+
+    ImageType::IndexType start;
+    start.Fill(0);
+
+    ImageType::RegionType region(start, size);
+
+    auto gdcmIO = itk::GDCMImageIO::New();
+    gdcmIO->KeepOriginalUIDOn();
+    auto writer = WriterType::New();
+    writer->SetImageIO(gdcmIO);
+
+    // Write each slice as a DICOM file with appropriate tags
+    for (size_t i = 0; i < positions.size(); ++i)
+    {
+      auto image = ImageType::New();
+      image->SetRegions(region);
+      image->Allocate();
+      image->FillBuffer(100); // Just to have nonzero pixel values
+
+      // Get the metadata dictionary
+      auto & dict = image->GetMetaDataDictionary();
+
+      // Set required DICOM tags
+      itk::EncapsulateMetaData<std::string>(dict, "0020|0032", positions[i]);          // ImagePositionPatient
+      itk::EncapsulateMetaData<std::string>(dict, "0020|0037", orientation);           // ImageOrientationPatient
+      itk::EncapsulateMetaData<std::string>(dict, "0008|0060", "CT");                  // Modality
+      itk::EncapsulateMetaData<std::string>(dict, "0020|0013", std::to_string(i + 1)); // InstanceNumber
+      itk::EncapsulateMetaData<std::string>(dict, "0010|0010", "Test^Patient");        // PatientName
+      itk::EncapsulateMetaData<std::string>(
+        dict, "0020|000e", "1.2.3.4.5.6.7.8"); // SeriesInstanceUID (same for all slices)
+
+      // Additional required DICOM tags for proper series
+      itk::EncapsulateMetaData<std::string>(
+        dict, "0008|0016", "1.2.840.10008.5.1.4.1.1.2"); // SOPClassUID (CT Image Storage)
+      itk::EncapsulateMetaData<std::string>(
+        dict, "0008|0018", "1.2.3.4.5.6.7.8.9." + std::to_string(i + 1));          // SOPInstanceUID
+      itk::EncapsulateMetaData<std::string>(dict, "0020|000d", "1.2.3.4.5.6.7.8"); // StudyInstanceUID
+      itk::EncapsulateMetaData<std::string>(dict, "0010|0020", "12345");           // PatientID
+      itk::EncapsulateMetaData<std::string>(dict, "0008|0020", "20240101");        // StudyDate
+      itk::EncapsulateMetaData<std::string>(dict, "0008|0030", "120000");          // StudyTime
+
+      const std::string filename = m_TempDir + "/slice_" + std::to_string(i + 1) + ".dcm";
+      writer->SetFileName(filename);
+      writer->SetInput(image);
+      writer->Update();
+
+      m_DicomFiles.push_back(filename);
+    }
+  }
+};
+
+} // namespace
+
+TEST_F(ITKGDCMSeriesTestData, ReadSlicesReverseOrder)
+{
+  // This test image series has non-unit meta-data:
+  // spacing: [0.859375, 0.85939, 1.60016]
+  // origin:[-112, -21.688, 126.894]
+  // direction:
+  //  [1 0 0,
+  //   0 0.466651 0.884442,
+  //   0 -0.884442 0.466651]
+  constexpr unsigned int VolumeDimension = 3;
+  using VolumeImageType = itk::Image<PixelType, VolumeDimension>;
+
+  using NamesGeneratorType = itk::GDCMSeriesFileNames;
+  auto namesGenerator = NamesGeneratorType::New();
+  namesGenerator->SetDirectory(m_DicomSeriesInput);
+  namesGenerator->SetUseSeriesDetails(true);
+  std::vector<std::string> fileNames = namesGenerator->GetInputFileNames();
+
+  using SeriesReaderType = itk::ImageSeriesReader<VolumeImageType>;
+  auto seriesReader = SeriesReaderType::New();
+  seriesReader->SetFileNames(fileNames);
+  auto gdcmIO = itk::GDCMImageIO::New();
+  seriesReader->SetImageIO(gdcmIO);
+  ASSERT_NO_THROW(seriesReader->UpdateLargestPossibleRegion());
+
+  VolumeImageType::Pointer outputImage = seriesReader->GetOutput();
+  outputImage->DisconnectPipeline();
+
+  seriesReader->ReverseOrderOn();
+  ASSERT_NO_THROW(seriesReader->UpdateLargestPossibleRegion());
+  VolumeImageType::Pointer reversedOutputImage = seriesReader->GetOutput();
+  reversedOutputImage->DisconnectPipeline();
+
+  std::cout << "baseline direction: " << outputImage->GetDirection() << std::endl;
+  std::cout << "reversed direction: " << reversedOutputImage->GetDirection() << std::endl;
+  EXPECT_EQ(outputImage->GetLargestPossibleRegion().GetSize(),
+            reversedOutputImage->GetLargestPossibleRegion().GetSize());
+  ITK_EXPECT_VECTOR_NEAR(outputImage->GetSpacing(), reversedOutputImage->GetSpacing(), 1e-6);
+  EXPECT_NE(outputImage->GetOrigin(), reversedOutputImage->GetOrigin());
+
+  // calculate the index at the middle of the image
+  VolumeImageType::IndexType middleIndex;
+  for (unsigned int d = 0; d < VolumeDimension; ++d)
+  {
+    middleIndex[d] =
+      outputImage->GetLargestPossibleRegion().GetIndex()[d] + outputImage->GetLargestPossibleRegion().GetSize()[d] / 2;
+  }
+
+  const std::vector<VolumeImageType::IndexType> testIndices = {
+    { { 0, 0, 0 } }, { { 1, 1, 1 } }, { { 2, 2, 2 } }, middleIndex
+  };
+
+  // test that the reversed image has the same pixel values at the same physical location
+  for (const auto & idx : testIndices)
+  {
+    VolumeImageType::PointType point;
+    outputImage->TransformIndexToPhysicalPoint(idx, point);
+    auto reverseIdx = reversedOutputImage->TransformPhysicalPointToIndex(point);
+
+    std::cout << "Testing idx: " << idx << " reverseIdx: " << reverseIdx << std::endl;
+    ASSERT_TRUE(reversedOutputImage->GetLargestPossibleRegion().IsInside(reverseIdx));
+    EXPECT_EQ(outputImage->GetPixel(idx), reversedOutputImage->GetPixel(reverseIdx));
+  }
+}
+
+TEST_F(ITKGDCMSeriesTestData, CreateAndReadTestSeries)
+{
+  // Verify that the DICOM files were created
+  ASSERT_EQ(m_DicomFiles.size(), 3);
+
+  for (const auto & filename : m_DicomFiles)
+  {
+    ASSERT_TRUE(itksys::SystemTools::FileExists(filename));
+  }
+
+  // Read the series using GDCMSeriesFileNames
+  using NamesGeneratorType = itk::GDCMSeriesFileNames;
+  auto namesGenerator = NamesGeneratorType::New();
+  namesGenerator->SetDirectory(m_TempDir);
+  namesGenerator->SetUseSeriesDetails(true);
+
+  std::vector<std::string> fileNames = namesGenerator->GetInputFileNames();
+  ASSERT_EQ(fileNames.size(), 3);
+
+  // Read the series
+  using ImageType3D = itk::Image<uint16_t, 3>;
+  using SeriesReaderType = itk::ImageSeriesReader<ImageType3D>;
+  auto seriesReader = SeriesReaderType::New();
+  seriesReader->SetFileNames(fileNames);
+
+  auto gdcmIO = itk::GDCMImageIO::New();
+  seriesReader->SetImageIO(gdcmIO);
+
+  ASSERT_NO_THROW(seriesReader->UpdateLargestPossibleRegion());
+
+  ImageType3D::Pointer image = seriesReader->GetOutput();
+
+  // Verify image properties
+  ImageType3D::SizeType expectedSize = { { 2, 2, 3 } };
+  EXPECT_EQ(image->GetLargestPossibleRegion().GetSize(), expectedSize);
+
+  // Verify pixel values (should be 100)
+  EXPECT_EQ(image->GetPixel({ 0, 0, 0 }), 100);
+  EXPECT_EQ(image->GetPixel({ 1, 1, 1 }), 100);
+}
+
+TEST_F(ITKGDCMSeriesTestData, ReadSeriesTopToBottom)
+{
+  // Read in top-to-bottom order (files ordered by ImagePositionPatient Z coordinate)
+  using ImageType3D = itk::Image<uint16_t, 3>;
+  using SeriesReaderType = itk::ImageSeriesReader<ImageType3D>;
+
+  // Get file list in top-to-bottom order
+  const std::vector<std::string> & filesTopToBottom = m_DicomFiles;
+
+  auto seriesReader = SeriesReaderType::New();
+  auto gdcmIO = itk::GDCMImageIO::New();
+  seriesReader->SetImageIO(gdcmIO);
+  seriesReader->SetFileNames(filesTopToBottom);
+  seriesReader->ForceOrthogonalDirectionOn(); // explicitly set default
+
+  ASSERT_NO_THROW(seriesReader->UpdateLargestPossibleRegion());
+  ImageType3D::Pointer image1 = seriesReader->GetOutput();
+
+  std::cout << "Top-to-bottom order:" << std::endl;
+  std::cout << "  Origin: " << image1->GetOrigin() << std::endl;
+  std::cout << "  Direction: " << image1->GetDirection() << std::endl;
+  std::cout << "  Spacing: " << image1->GetSpacing() << std::endl;
+
+  // Verify image properties
+  ImageType3D::SizeType expectedSize = { { 2, 2, 3 } };
+  EXPECT_EQ(image1->GetLargestPossibleRegion().GetSize(), expectedSize);
+
+  // The origin should be at the position of the first slice (top slice)
+  ImageType3D::PointType expectedOrigin{ { -216.500, -216.500, 70.000 } }; // Z position of slice 1 (top)
+
+  ITK_EXPECT_VECTOR_NEAR(image1->GetOrigin(), expectedOrigin, 1e-3);
+
+  // Z spacing should be positive
+  EXPECT_GT(image1->GetSpacing()[2], 0.0);
+  // but the direction should have a negative Z component
+  EXPECT_LT(image1->GetDirection()[2][2], 0.0);
+}
+
+TEST_F(ITKGDCMSeriesTestData, ReadSeriesBottomToTop)
+{
+  // Read in bottom-to-top order (files ordered by ImagePositionPatient Z coordinate)
+  using ImageType3D = itk::Image<uint16_t, 3>;
+  using SeriesReaderType = itk::ImageSeriesReader<ImageType3D>;
+
+  // Get file list in bottom-to-top order
+  std::vector<std::string> filesBottomToTop(m_DicomFiles.rbegin(), m_DicomFiles.rend());
+
+  auto seriesReader = SeriesReaderType::New();
+  auto gdcmIO = itk::GDCMImageIO::New();
+  seriesReader->SetImageIO(gdcmIO);
+  seriesReader->SetFileNames(filesBottomToTop);
+  seriesReader->ForceOrthogonalDirectionOn(); // explicitly set default
+
+  ASSERT_NO_THROW(seriesReader->UpdateLargestPossibleRegion());
+  ImageType3D::Pointer image2 = seriesReader->GetOutput();
+
+  std::cout << "Bottom-to-top order:" << std::endl;
+  std::cout << "  Origin: " << image2->GetOrigin() << std::endl;
+  std::cout << "  Direction: " << image2->GetDirection() << std::endl;
+  std::cout << "  Spacing: " << image2->GetSpacing() << std::endl;
+
+  // Verify image properties
+  ImageType3D::SizeType expectedSize = { { 2, 2, 3 } };
+  EXPECT_EQ(image2->GetLargestPossibleRegion().GetSize(), expectedSize);
+
+  // The origin should be at the position of the first slice (bottom slice)
+  ImageType3D::PointType expectedOrigin{
+    { -216.500, -216.500, -445.000 }
+  }; // X,Y from slice 1, Z position of slice 3 (bottom)
+
+  ITK_EXPECT_VECTOR_NEAR(image2->GetOrigin(), expectedOrigin, 1e-3);
+
+  // Z spacing should be positive (going from bottom to top)
+  EXPECT_GT(image2->GetSpacing()[2], 0.0);
+  // and the direction should have a positive Z component
+  EXPECT_GT(image2->GetDirection()[2][2], 0.0);
+}