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.

Author: blowekamp

Modules/IO/GDCM/test/CMakeLists.txt

@@ -439,3 +439,10 @@ 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}")
+
+target_compile_definitions(ITKGDCMImageIOGTestDriver PRIVATE "-DITK_TEST_OUTPUT_DIR=${ITK_TEST_OUTPUT_DIR}")

Modules/IO/GDCM/test/itkGDCMImageIOGTest.cxx

@@ -0,0 +1,157 @@
+/*=========================================================================
+ *
+ *  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 "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);
+  }
+
+  std::string m_TempDir = TOSTRING(ITK_TEST_OUTPUT_DIR) + "/ITKGDCMImageIO";
+};
+
+} // namespace
+
+TEST_F(ITKGDCMImageIO, ReadSlicesNegativeSpacing)
+{
+  // Tests the case where the slices are written with negative space between slices,
+  // and the reader should be able to read them back correctly with negative in the direction cosine matrix.
+
+
+  using PixelType = uint16_t;
+  constexpr unsigned int Dimension = 3;
+  using ImageType = itk::Image<PixelType, Dimension>;
+  using SliceType = itk::Image<PixelType, 2>;
+
+  // Create a 3D image of size [2,2,3] and fill with 100
+  ImageType::SizeType   size = { { 2, 2, 3 } };
+  ImageType::RegionType region;
+  region.SetSize(size);
+  auto image = ImageType::New();
+  image->SetRegions(region);
+  image->Allocate();
+  image->FillBuffer(100);
+
+  // DICOM meta-data values
+  const std::vector<std::string> positions = { "-216.500\\-216.500\\70.000",
+                                               "-216.500\\-216.500\\-187.500",
+                                               "-216.500\\-216.500\\-445.000" };
+  const std::string              orientation = "1.000000\\0.000000\\0.000000\\0.000000\\1.000000\\0.000000";
+  const std::string              pixelSpacing = "1.1\\1.2";
+
+  std::string seriesDir = m_TempDir + "/" + ::testing::UnitTest::GetInstance()->current_test_info()->name();
+  itksys::SystemTools::MakeDirectory(seriesDir);
+
+  using WriterType = itk::ImageFileWriter<SliceType>;
+  using ImageIOType = itk::GDCMImageIO;
+
+  for (unsigned int i = 0; i < size[2]; ++i)
+  {
+    // Extract 2D slice
+    SliceType::RegionType sliceRegion;
+    SliceType::SizeType   sliceSize = { { size[0], size[1] } };
+    SliceType::IndexType  sliceIndex = { { 0, 0 } };
+    sliceRegion.SetSize(sliceSize);
+    sliceRegion.SetIndex(sliceIndex);
+
+    SliceType::Pointer slice = SliceType::New();
+    slice->SetRegions(sliceRegion);
+    slice->Allocate();
+
+    slice->FillBuffer(5 + i);
+
+    // Set image spacing
+    ImageType::SpacingType spacing{ { 1.1, 1.2, 257.5 } };
+    image->SetSpacing(spacing);
+
+
+    // Set required DICOM tags
+    itk::MetaDataDictionary & dict = slice->GetMetaDataDictionary();
+    itk::EncapsulateMetaData<std::string>(dict, "0020|0032", positions[i]);
+    itk::EncapsulateMetaData<std::string>(dict, "0020|0037", orientation);
+    itk::EncapsulateMetaData<std::string>(dict, "0028|0030", pixelSpacing);
+
+    // Write the slice
+    auto               writer = WriterType::New();
+    auto               gdcmIO = ImageIOType::New();
+    std::ostringstream filename;
+    filename << seriesDir << "/slice_" << (i + 1) << ".dcm";
+    writer->SetFileName(filename.str());
+    writer->SetInput(slice);
+    writer->SetImageIO(gdcmIO);
+    ASSERT_NO_THROW(writer->Update());
+  }
+
+  std::vector<std::string> fileNames;
+  for (unsigned int i = 0; i < size[2]; ++i)
+  {
+    std::ostringstream filename;
+    filename << seriesDir << "/slice_" << (i + 1) << ".dcm";
+    fileNames.push_back(filename.str());
+  }
+
+  // Reader the slices back
+  using SeriesReaderType = itk::ImageSeriesReader<ImageType>;
+  auto seriesReader = SeriesReaderType::New();
+  seriesReader->SetFileNames(fileNames);
+  auto gdcmIO = ImageIOType::New();
+  seriesReader->SetImageIO(gdcmIO);
+  ASSERT_NO_THROW(seriesReader->Update());
+  auto outputImage = seriesReader->GetOutput();
+
+  // Check the spacing, origin, and direction
+  ImageType::SpacingType   spacing = outputImage->GetSpacing();
+  ImageType::PointType     origin = outputImage->GetOrigin();
+  ImageType::DirectionType direction = outputImage->GetDirection();
+  ImageType::DirectionType expectedDirection;
+  expectedDirection.SetIdentity();
+  expectedDirection[2][2] = -1.0; // Set the z-direction to -1.0
+  ImageType::PointType   expectedOrigin{ { -216.5, -216.5, 70.0 } };
+  ImageType::SpacingType expectedSpacing{ { 1.1, 1.2, 257.5 } };
+
+  std::cout << "Spacing: " << spacing << std::endl;
+
+  EXPECT_EQ(spacing, expectedSpacing);
+  EXPECT_EQ(origin, expectedOrigin);
+  EXPECT_EQ(direction, expectedDirection);
+}

Modules/IO/ImageBase/include/itkImageSeriesReader.hxx

@@ -184,7 +184,14 @@ ImageSeriesReader<TOutputImage>::GenerateOutputInformation()
     {
       spacing[this->m_NumberOfDimensionsInImage] = dirNnorm / (numberOfFiles - 1);
       this->m_SpacingDefined = true;
-      if (!m_ForceOrthogonalDirection)
+      if (m_ForceOrthogonalDirection)
+      {
+        for (unsigned int j = 0; j < TOutputImage::ImageDimension; ++j)
+        {
+          direction[j][this->m_NumberOfDimensionsInImage] *= std::signbit(dirN[j]) ? -1.0 : 1.0;
+        }
+      }
+      else
       {
         for (unsigned int j = 0; j < TOutputImage::ImageDimension; ++j)
         {