From d14a0b2815a13eec4a9b0e82c9d3f02730088c07 Mon Sep 17 00:00:00 2001 From: Bradley Lowekamp Date: Wed, 29 Apr 2026 14:36:23 -0500 Subject: [PATCH 1/2] 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. --- Modules/IO/GDCM/test/CMakeLists.txt | 25 ++ Modules/IO/GDCM/test/itkGDCMImageIOGTest.cxx | 333 ++++++++++++++++++ .../include/itkImageSeriesReader.hxx | 30 +- 3 files changed, 385 insertions(+), 3 deletions(-) create mode 100644 Modules/IO/GDCM/test/itkGDCMImageIOGTest.cxx diff --git a/Modules/IO/GDCM/test/CMakeLists.txt b/Modules/IO/GDCM/test/CMakeLists.txt index 630fbded88e..5c2e107c232 100644 --- a/Modules/IO/GDCM/test/CMakeLists.txt +++ b/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 +) diff --git a/Modules/IO/GDCM/test/itkGDCMImageIOGTest.cxx b/Modules/IO/GDCM/test/itkGDCMImageIOGTest.cxx new file mode 100644 index 00000000000..089ca26bc78 --- /dev/null +++ b/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 +#include + + +#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; + using WriterType = itk::ImageFileWriter; + + 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 m_DicomFiles; + + +private: + void + CreateTestDicomSeries() + { + // DICOM meta-data values from the report + const std::vector 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(dict, "0020|0032", positions[i]); // ImagePositionPatient + itk::EncapsulateMetaData(dict, "0020|0037", orientation); // ImageOrientationPatient + itk::EncapsulateMetaData(dict, "0008|0060", "CT"); // Modality + itk::EncapsulateMetaData(dict, "0020|0013", std::to_string(i + 1)); // InstanceNumber + itk::EncapsulateMetaData(dict, "0010|0010", "Test^Patient"); // PatientName + itk::EncapsulateMetaData( + 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( + dict, "0008|0016", "1.2.840.10008.5.1.4.1.1.2"); // SOPClassUID (CT Image Storage) + itk::EncapsulateMetaData( + dict, "0008|0018", "1.2.3.4.5.6.7.8.9." + std::to_string(i + 1)); // SOPInstanceUID + itk::EncapsulateMetaData(dict, "0020|000d", "1.2.3.4.5.6.7.8"); // StudyInstanceUID + itk::EncapsulateMetaData(dict, "0010|0020", "12345"); // PatientID + itk::EncapsulateMetaData(dict, "0008|0020", "20240101"); // StudyDate + itk::EncapsulateMetaData(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; + + using NamesGeneratorType = itk::GDCMSeriesFileNames; + auto namesGenerator = NamesGeneratorType::New(); + namesGenerator->SetDirectory(m_DicomSeriesInput); + namesGenerator->SetUseSeriesDetails(true); + std::vector fileNames = namesGenerator->GetInputFileNames(); + + using SeriesReaderType = itk::ImageSeriesReader; + 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 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 fileNames = namesGenerator->GetInputFileNames(); + ASSERT_EQ(fileNames.size(), 3); + + // Read the series + using ImageType3D = itk::Image; + using SeriesReaderType = itk::ImageSeriesReader; + 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; + using SeriesReaderType = itk::ImageSeriesReader; + + // Get file list in top-to-bottom order + const std::vector & 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; + using SeriesReaderType = itk::ImageSeriesReader; + + // Get file list in bottom-to-top order + std::vector 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); +} diff --git a/Modules/IO/ImageBase/include/itkImageSeriesReader.hxx b/Modules/IO/ImageBase/include/itkImageSeriesReader.hxx index cba45ad6a25..35fa7222871 100644 --- a/Modules/IO/ImageBase/include/itkImageSeriesReader.hxx +++ b/Modules/IO/ImageBase/include/itkImageSeriesReader.hxx @@ -172,8 +172,8 @@ ImageSeriesReader::GenerateOutputInformation() // Override the position if there is an ITK_ImageOrigin ExposeMetaData>(lastReader->GetImageIO()->GetMetaDataDictionary(), key, positionN); - // Compute and set the inter slice spacing - // and last (usually third) axis of direction + // Compute and set the inter-slice spacing + // and last (usually third) axis of a direction Vector dirN; for (unsigned int j = 0; j < TOutputImage::ImageDimension; ++j) { @@ -187,9 +187,33 @@ ImageSeriesReader::GenerateOutputInformation() } else { + // always positive spacing, corrected in the direction spacing[this->m_NumberOfDimensionsInImage] = dirNnorm / (numberOfFiles - 1); this->m_SpacingDefined = true; - if (!m_ForceOrthogonalDirection) + if (m_ForceOrthogonalDirection) + { + // DICOM files only specify two direction vectors, and the third is derived as orthogonal to those. Hence only + // when loading DICOM files, the third direction is forced to be orthogonal from the DICOM ImageIO. However, + // with the forced positive spacing the sign of the direction must be checked. + + + // Compute dot product between dirN and the current direction column vector + SpacingScalarType dotProduct = 0.0; + for (unsigned int j = 0; j < TOutputImage::ImageDimension; ++j) + { + dotProduct += dirN[j] * direction[j][this->m_NumberOfDimensionsInImage]; + } + + // If dot product is negative, flip the sign of the entire direction column vector + if (dotProduct < 0.0) + { + for (unsigned int j = 0; j < TOutputImage::ImageDimension; ++j) + { + direction[j][this->m_NumberOfDimensionsInImage] *= -1.0; + } + } + } + else { for (unsigned int j = 0; j < TOutputImage::ImageDimension; ++j) { From fcecee9e4c6af3281a7808d7a168a52ef583e1bf Mon Sep 17 00:00:00 2001 From: Hans Johnson Date: Wed, 29 Apr 2026 14:36:23 -0500 Subject: [PATCH 2/2] =?UTF-8?q?ENH:=20Tests=20=E2=80=94=20Extensive=20regr?= =?UTF-8?q?ession=20suite=20for=20PR=205357=20direction-cosine=20fix?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Modules/IO/GDCM/test/CMakeLists.txt | 6 +- .../GDCM/test/itkGDCMSeriesDirectionGTest.cxx | 826 ++++++++++++++++++ 2 files changed, 831 insertions(+), 1 deletion(-) create mode 100644 Modules/IO/GDCM/test/itkGDCMSeriesDirectionGTest.cxx diff --git a/Modules/IO/GDCM/test/CMakeLists.txt b/Modules/IO/GDCM/test/CMakeLists.txt index 5c2e107c232..7ba5b6f57d9 100644 --- a/Modules/IO/GDCM/test/CMakeLists.txt +++ b/Modules/IO/GDCM/test/CMakeLists.txt @@ -449,7 +449,11 @@ addcompliancetest(RLE-RGB) addcompliancetest(HTJ2K-YBR_ICT Lily_full.mha) addcompliancetest(HTJ2K-YBR_RCT Lily_full.mha) -set(ITKGDCMImageIOGTests itkGDCMImageIOGTest.cxx) +set( + ITKGDCMImageIOGTests + itkGDCMImageIOGTest.cxx + itkGDCMSeriesDirectionGTest.cxx +) creategoogletestdriver(ITKGDCMImageIO "${ITKIOGDCM-Test_LIBRARIES}" "${ITKGDCMImageIOGTests}") ExternalData_Expand_Arguments( diff --git a/Modules/IO/GDCM/test/itkGDCMSeriesDirectionGTest.cxx b/Modules/IO/GDCM/test/itkGDCMSeriesDirectionGTest.cxx new file mode 100644 index 00000000000..bb3f9133eee --- /dev/null +++ b/Modules/IO/GDCM/test/itkGDCMSeriesDirectionGTest.cxx @@ -0,0 +1,826 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +// Extended regression suite for ImageSeriesReader direction-cosine handling +// when the user-supplied file ordering implies a slice-stacking direction +// opposite to the GDCM-derived 3rd column of the direction matrix. +// +// Bug references: +// - SimpleITK#2292 — bottom-to-top input produced direction[2][2] = -1 +// when origin was simultaneously set to the bottom slice +// (origin/direction inconsistency). +// - ITKElastix#291 — direction-cosine 3rd column not [0,0,1] when expected. +// - ITK PR #5357 — proposed fix in itkImageSeriesReader.hxx. +// +// The fix flips the sign of the 3rd direction column when +// dot(positionN - position1, direction[*][N-1]) < 0. +// These tests pin the contract from multiple angles: +// * axis-aligned vs oblique row/col orientation +// * axial / sagittal / coronal patient frames +// * forward, reversed, and ReverseOrderOn() input +// * orthonormality, right-handedness (det == +1) +// * ForceOrthogonalDirection on/off symmetry +// * physical-point invariance for every voxel under reordering + +#include "itkImage.h" +#include "itkImageFileWriter.h" +#include "itkGDCMImageIO.h" +#include "itkGDCMSeriesFileNames.h" +#include "itkMetaDataObject.h" +#include "itkGTest.h" +#include "itkMath.h" +#include "itksys/SystemTools.hxx" +#include "itkImageSeriesReader.h" +#include "itkImageRegionConstIteratorWithIndex.h" +#include "itkImageRegionIteratorWithIndex.h" +#include +#include +#include +#include +#include +#include + +#define _STRING(s) #s +#define TOSTRING(s) std::string(_STRING(s)) + +namespace +{ + +constexpr double kGeomTol = 1e-3; // mm — DICOM IPP precision +constexpr double kUnitTol = 1e-6; // unit-vector tolerance + +using PixelT = uint16_t; +using Slice2D = itk::Image; +using Volume3D = itk::Image; + +// Format a 3-component vector as DICOM DS string ("a\b\c"). +std::string +formatDS(double a, double b, double c) +{ + std::ostringstream os; + os.precision(6); + os << std::fixed << a << '\\' << b << '\\' << c; + return os.str(); +} + +std::string +formatOrientation(const std::array & row, const std::array & col) +{ + std::ostringstream os; + os.precision(6); + os << std::fixed << row[0] << '\\' << row[1] << '\\' << row[2] << '\\' << col[0] << '\\' << col[1] << '\\' << col[2]; + return os.str(); +} + +// Cross product for std::array. +std::array +cross(const std::array & a, const std::array & b) +{ + return { { a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0] } }; +} + +double +norm3(const std::array & v) +{ + return std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); +} + +std::array +normalize3(const std::array & v) +{ + const double n = norm3(v); + return { { v[0] / n, v[1] / n, v[2] / n } }; +} + +// Synthesize a DICOM series by writing one 2D slice per file with controlled +// ImagePositionPatient / ImageOrientationPatient tags. Returns full file +// paths in the order written. +struct SeriesSpec +{ + std::string dir; + std::array rowOrient; // 0020,0037 first triplet + std::array colOrient; // 0020,0037 second triplet + std::array pos0; // first slice IPP + std::array stepIPP; // delta between successive slices + unsigned int numSlices = 5; + itk::Size<2> pixelSize = { { 4, 3 } }; + std::string seriesUID = "1.2.826.0.1.5357"; + std::string modality = "CT"; +}; + +std::vector +writeSyntheticSeries(const SeriesSpec & spec) +{ + itksys::SystemTools::MakeDirectory(spec.dir); + const std::string orientation = formatOrientation(spec.rowOrient, spec.colOrient); + + Slice2D::SizeType sz = spec.pixelSize; + Slice2D::IndexType start; + start.Fill(0); + Slice2D::RegionType region(start, sz); + + std::vector files; + files.reserve(spec.numSlices); + + for (unsigned int i = 0; i < spec.numSlices; ++i) + { + auto image = Slice2D::New(); + image->SetRegions(region); + image->Allocate(); + + // Encode the slice index into the pixel buffer so we can verify pixel + // identity after physical-point lookup: pixel = 1000*slice + 10*y + x. + itk::ImageRegionIteratorWithIndex it(image, region); + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + const auto idx = it.GetIndex(); + it.Set( + static_cast(1000U * (i + 1U) + 10U * static_cast(idx[1]) + static_cast(idx[0]))); + } + + auto & dict = image->GetMetaDataDictionary(); + + const std::string ipp = formatDS( + spec.pos0[0] + i * spec.stepIPP[0], spec.pos0[1] + i * spec.stepIPP[1], spec.pos0[2] + i * spec.stepIPP[2]); + + itk::EncapsulateMetaData(dict, "0020|0032", ipp); + itk::EncapsulateMetaData(dict, "0020|0037", orientation); + itk::EncapsulateMetaData(dict, "0008|0060", spec.modality); + itk::EncapsulateMetaData(dict, "0020|0013", std::to_string(i + 1)); + itk::EncapsulateMetaData(dict, "0010|0010", "PR5357^Test"); + itk::EncapsulateMetaData(dict, "0020|000e", spec.seriesUID); + itk::EncapsulateMetaData(dict, "0020|000d", spec.seriesUID + ".1"); + itk::EncapsulateMetaData(dict, "0010|0020", "PR5357"); + itk::EncapsulateMetaData(dict, "0008|0020", "20260101"); + itk::EncapsulateMetaData(dict, "0008|0030", "120000"); + itk::EncapsulateMetaData(dict, "0008|0016", "1.2.840.10008.5.1.4.1.1.2"); + itk::EncapsulateMetaData(dict, "0008|0018", spec.seriesUID + "." + std::to_string(i + 1)); + + auto gdcmIO = itk::GDCMImageIO::New(); + gdcmIO->KeepOriginalUIDOn(); + + auto writer = itk::ImageFileWriter::New(); + writer->SetImageIO(gdcmIO); + + std::ostringstream fname; + fname << spec.dir << "/slice_" << std::setw(4) << std::setfill('0') << i << ".dcm"; + writer->SetFileName(fname.str()); + writer->SetInput(image); + writer->Update(); + files.push_back(fname.str()); + } + return files; +} + +// Read a series with the given file ordering and ForceOrthogonalDirection setting. +Volume3D::Pointer +readSeries(const std::vector & files, bool forceOrtho = true, bool reverseOrder = false) +{ + auto reader = itk::ImageSeriesReader::New(); + auto io = itk::GDCMImageIO::New(); + reader->SetImageIO(io); + reader->SetFileNames(files); + if (forceOrtho) + { + reader->ForceOrthogonalDirectionOn(); + } + else + { + reader->ForceOrthogonalDirectionOff(); + } + if (reverseOrder) + { + reader->ReverseOrderOn(); + } + reader->UpdateLargestPossibleRegion(); + Volume3D::Pointer out = reader->GetOutput(); + out->DisconnectPipeline(); + return out; +} + +// Determinant of a 3×3 itk::Matrix-like type with operator[][]. +template +double +det3(const M & m) +{ + return m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1]) - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) + + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]); +} + +// Assert columns are unit length and mutually orthogonal. Determinant must be +// ±1 (orthonormal); ITK allows left-handed direction matrices, so the sign is +// not constrained. +template +void +expectOrthonormal(const M & d) +{ + for (unsigned int c = 0; c < 3; ++c) + { + double n2 = 0.0; + for (unsigned int r = 0; r < 3; ++r) + { + n2 += d[r][c] * d[r][c]; + } + EXPECT_NEAR(n2, 1.0, kUnitTol) << "column " << c << " not unit length"; + } + for (unsigned int c1 = 0; c1 < 3; ++c1) + { + for (unsigned int c2 = c1 + 1; c2 < 3; ++c2) + { + double dot = 0.0; + for (unsigned int r = 0; r < 3; ++r) + { + dot += d[r][c1] * d[r][c2]; + } + EXPECT_NEAR(dot, 0.0, kUnitTol) << "columns " << c1 << "," << c2 << " not orthogonal"; + } + } + EXPECT_NEAR(std::abs(det3(d)), 1.0, kUnitTol) << "direction det not ±1"; +} + +// Verify that for every voxel in 'a', its physical point maps to a valid voxel +// in 'b' carrying the identical pixel value. +void +expectPhysicalPointInvariant(const Volume3D::Pointer & a, const Volume3D::Pointer & b) +{ + ASSERT_EQ(a->GetLargestPossibleRegion().GetSize(), b->GetLargestPossibleRegion().GetSize()); + itk::ImageRegionConstIteratorWithIndex it(a, a->GetLargestPossibleRegion()); + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + Volume3D::PointType p; + a->TransformIndexToPhysicalPoint(it.GetIndex(), p); + auto bIdx = b->TransformPhysicalPointToIndex(p); + ASSERT_TRUE(b->GetLargestPossibleRegion().IsInside(bIdx)) + << "physical point " << p << " from index " << it.GetIndex() << " maps outside b"; + EXPECT_EQ(it.Get(), b->GetPixel(bIdx)) + << "pixel mismatch at physical point " << p << " a-idx " << it.GetIndex() << " b-idx " << bIdx; + } +} + +struct GDCMSeriesDirection : public ::testing::Test +{ + void + SetUp() override + { + itksys::SystemTools::MakeDirectory(m_TempDir); + } + void + TearDown() override + { + itksys::SystemTools::RemoveADirectory(m_TempDir); + } + std::string + makeCaseDir(const std::string & tag) + { + const std::string d = m_TempDir + "/" + tag; + itksys::SystemTools::MakeDirectory(d); + return d; + } + const std::string m_TempDir{ TOSTRING(ITK_TEST_OUTPUT_DIR) + "/GDCMSeriesDirection" }; +}; + +} // namespace + + +// --------------------------------------------------------------------------- +// 1. The exact SimpleITK#2292 reproduction: monotonically decreasing Z, +// axis-aligned axial orientation, given to the reader bottom-to-top. +// --------------------------------------------------------------------------- +TEST_F(GDCMSeriesDirection, SimpleITK2292_BottomToTop_OriginAndDirectionAgree) +{ + SeriesSpec spec; + spec.dir = makeCaseDir("simpleitk2292_btt"); + spec.rowOrient = { { 1.0, 0.0, 0.0 } }; + spec.colOrient = { { 0.0, 1.0, 0.0 } }; + spec.pos0 = { { -216.5, -216.5, 70.0 } }; // top slice + spec.stepIPP = { { 0.0, 0.0, -1.25 } }; // moving down in Z + spec.numSlices = 9; + + auto files = writeSyntheticSeries(spec); + // Hand the reader the BOTTOM slice first (reverse the natural write order). + std::vector btt(files.rbegin(), files.rend()); + + auto img = readSeries(btt); + expectOrthonormal(img->GetDirection()); + + // Bottom slice was written at z = 70 + 8*(-1.25) = 60. + EXPECT_NEAR(img->GetOrigin()[2], 60.0, kGeomTol); + // Direction Z must point UP (positive) to be consistent with origin = bottom + // slice and forced-positive spacing. + EXPECT_GT(img->GetDirection()[2][2], 0.0) << "Origin is the bottom slice; direction Z must be +1 to reach the top."; + + // Walking from origin by spacing*(N-1) along direction[*][2] must land on the top slice. + Volume3D::IndexType last; + last[0] = 0; + last[1] = 0; + last[2] = static_cast(img->GetLargestPossibleRegion().GetSize()[2]) - 1; + Volume3D::PointType p; + img->TransformIndexToPhysicalPoint(last, p); + EXPECT_NEAR(p[2], 70.0, kGeomTol); +} + +// --------------------------------------------------------------------------- +// 2. Same series, top-to-bottom input — origin is the top slice (Z=70), +// direction Z must be negative. +// --------------------------------------------------------------------------- +TEST_F(GDCMSeriesDirection, SimpleITK2292_TopToBottom_OriginAndDirectionAgree) +{ + SeriesSpec spec; + spec.dir = makeCaseDir("simpleitk2292_ttb"); + spec.rowOrient = { { 1.0, 0.0, 0.0 } }; + spec.colOrient = { { 0.0, 1.0, 0.0 } }; + spec.pos0 = { { -216.5, -216.5, 70.0 } }; + spec.stepIPP = { { 0.0, 0.0, -1.25 } }; + spec.numSlices = 9; + + auto files = writeSyntheticSeries(spec); + auto img = readSeries(files); + expectOrthonormal(img->GetDirection()); + + EXPECT_NEAR(img->GetOrigin()[2], 70.0, kGeomTol); + EXPECT_LT(img->GetDirection()[2][2], 0.0); +} + +// --------------------------------------------------------------------------- +// 3. Forward and reverse readings must place every voxel at the same physical +// point — this is the contract the patch promises. +// --------------------------------------------------------------------------- +TEST_F(GDCMSeriesDirection, PhysicalPointInvariance_AxialReversal) +{ + SeriesSpec spec; + spec.dir = makeCaseDir("axial_invariance"); + spec.rowOrient = { { 1.0, 0.0, 0.0 } }; + spec.colOrient = { { 0.0, 1.0, 0.0 } }; + spec.pos0 = { { 0.0, 0.0, 0.0 } }; + spec.stepIPP = { { 0.0, 0.0, 2.5 } }; + spec.numSlices = 6; + + auto files = writeSyntheticSeries(spec); + std::vector reversed(files.rbegin(), files.rend()); + + auto fwd = readSeries(files); + auto rev = readSeries(reversed); + + expectOrthonormal(fwd->GetDirection()); + expectOrthonormal(rev->GetDirection()); + ITK_EXPECT_VECTOR_NEAR(fwd->GetSpacing(), rev->GetSpacing(), kGeomTol); + expectPhysicalPointInvariant(fwd, rev); +} + +// --------------------------------------------------------------------------- +// 4. ReverseOrderOn() flag — same input, two interpretations, must agree on +// physical-point identity. This is the trivially-reproducible case +// @blowekamp identified. +// --------------------------------------------------------------------------- +TEST_F(GDCMSeriesDirection, PhysicalPointInvariance_ReverseOrderFlag) +{ + SeriesSpec spec; + spec.dir = makeCaseDir("reverse_order_flag"); + spec.rowOrient = { { 1.0, 0.0, 0.0 } }; + spec.colOrient = { { 0.0, 1.0, 0.0 } }; + spec.pos0 = { { 10.0, -5.0, 100.0 } }; + spec.stepIPP = { { 0.0, 0.0, 3.0 } }; + spec.numSlices = 7; + + auto files = writeSyntheticSeries(spec); + auto fwd = readSeries(files, true, false); + auto rev = readSeries(files, true, true); + + expectOrthonormal(fwd->GetDirection()); + expectOrthonormal(rev->GetDirection()); + ITK_EXPECT_VECTOR_NEAR(fwd->GetSpacing(), rev->GetSpacing(), kGeomTol); + EXPECT_NE(fwd->GetOrigin(), rev->GetOrigin()); + expectPhysicalPointInvariant(fwd, rev); +} + +// --------------------------------------------------------------------------- +// 5. Sagittal acquisition: row=AP, col=SI, slice-stack = LR. +// Reverse the file order; direction's 1st column (which by GDCM +// construction equals row-orientation) is unchanged, but the slice-stack +// column (3rd) must flip sign. +// --------------------------------------------------------------------------- +TEST_F(GDCMSeriesDirection, Sagittal_AxisAligned_Reversal) +{ + SeriesSpec spec; + spec.dir = makeCaseDir("sagittal"); + spec.rowOrient = { { 0.0, 1.0, 0.0 } }; // patient anterior + spec.colOrient = { { 0.0, 0.0, -1.0 } }; // patient inferior (DICOM convention) + spec.pos0 = { { -100.0, -120.0, 80.0 } }; + spec.stepIPP = { { 1.0, 0.0, 0.0 } }; // moving in +X + spec.numSlices = 5; + + auto files = writeSyntheticSeries(spec); + std::vector rev(files.rbegin(), files.rend()); + + auto fwd = readSeries(files); + auto rv = readSeries(rev); + expectOrthonormal(fwd->GetDirection()); + expectOrthonormal(rv->GetDirection()); + + // dirN sign-flip must reflect in the 3rd column. + EXPECT_NEAR(fwd->GetDirection()[0][2], -rv->GetDirection()[0][2], kUnitTol); + EXPECT_NEAR(fwd->GetDirection()[1][2], -rv->GetDirection()[1][2], kUnitTol); + EXPECT_NEAR(fwd->GetDirection()[2][2], -rv->GetDirection()[2][2], kUnitTol); + expectPhysicalPointInvariant(fwd, rv); +} + +// --------------------------------------------------------------------------- +// 6. Coronal acquisition: row=LR, col=SI, slice-stack = AP. +// --------------------------------------------------------------------------- +TEST_F(GDCMSeriesDirection, Coronal_AxisAligned_Reversal) +{ + SeriesSpec spec; + spec.dir = makeCaseDir("coronal"); + spec.rowOrient = { { 1.0, 0.0, 0.0 } }; + spec.colOrient = { { 0.0, 0.0, -1.0 } }; + spec.pos0 = { { -50.0, 30.0, 5.0 } }; + spec.stepIPP = { { 0.0, -2.0, 0.0 } }; // posterior-to-anterior reversed + spec.numSlices = 6; + + auto files = writeSyntheticSeries(spec); + std::vector rev(files.rbegin(), files.rend()); + + auto fwd = readSeries(files); + auto rv = readSeries(rev); + expectOrthonormal(fwd->GetDirection()); + expectOrthonormal(rv->GetDirection()); + expectPhysicalPointInvariant(fwd, rv); +} + +// --------------------------------------------------------------------------- +// 7. Oblique acquisition (no zeros, no ones in the orientation triplets). +// This is the explicit case @thewtex requested. Row and column +// orientations are derived from a small pitch-yaw rotation of the +// canonical axial frame. +// --------------------------------------------------------------------------- +TEST_F(GDCMSeriesDirection, Oblique_NonAxisAligned_Reversal) +{ + // Build an oblique row/col by rotating the axial (e_x, e_y) about an axis + // that mixes all three world components. Use angles chosen so no element + // of the resulting row/col is 0 or ±1. + const double a = 0.27; // ~15.5° + const double b = 0.41; // ~23.5° + const double ca = std::cos(a), sa = std::sin(a); + const double cb = std::cos(b), sb = std::sin(b); + // Compose Rz(a) * Rx(b) applied to (1,0,0) and (0,1,0). + std::array row = { { ca, sa * cb, sa * sb } }; + std::array col = { { -sa, ca * cb, ca * sb } }; + row = normalize3(row); + col = normalize3(col); + // Slice stacking direction = row × col, then offset every slice along it. + const std::array n = normalize3(cross(row, col)); + + SeriesSpec spec; + spec.dir = makeCaseDir("oblique"); + spec.rowOrient = row; + spec.colOrient = col; + spec.pos0 = { { -45.5, 12.25, 88.125 } }; + spec.stepIPP = { { 1.7 * n[0], 1.7 * n[1], 1.7 * n[2] } }; + spec.numSlices = 5; + + auto files = writeSyntheticSeries(spec); + std::vector rev(files.rbegin(), files.rend()); + + auto fwd = readSeries(files); + auto rv = readSeries(rev); + expectOrthonormal(fwd->GetDirection()); + expectOrthonormal(rv->GetDirection()); + + // The first two columns derive from the (unchanged) DICOM ImageOrientationPatient + // tags and must NOT change between forward and reversed readings. + for (unsigned int r = 0; r < 3; ++r) + { + EXPECT_NEAR(fwd->GetDirection()[r][0], rv->GetDirection()[r][0], kUnitTol) + << "row-orient (col 0) drifted at row " << r; + EXPECT_NEAR(fwd->GetDirection()[r][1], rv->GetDirection()[r][1], kUnitTol) + << "col-orient (col 1) drifted at row " << r; + } + // The 3rd column must flip sign under reversal. + for (unsigned int r = 0; r < 3; ++r) + { + EXPECT_NEAR(fwd->GetDirection()[r][2], -rv->GetDirection()[r][2], kUnitTol) + << "slice-stack (col 2) sign-flip failed at row " << r; + } + expectPhysicalPointInvariant(fwd, rv); +} + +// --------------------------------------------------------------------------- +// 8. Reviewer's specific request: third-column equals cross-product of the +// first two columns (after the fix is applied). +// --------------------------------------------------------------------------- +TEST_F(GDCMSeriesDirection, ThirdColumnEqualsCrossProduct_Oblique) +{ + const double a = 0.31, b = 0.19; + const double ca = std::cos(a), sa = std::sin(a); + const double cb = std::cos(b), sb = std::sin(b); + std::array row = normalize3({ { ca, sa * cb, sa * sb } }); + std::array col = normalize3({ { -sa, ca * cb, ca * sb } }); + std::array n = normalize3(cross(row, col)); + + SeriesSpec spec; + spec.dir = makeCaseDir("crossprod"); + spec.rowOrient = row; + spec.colOrient = col; + spec.pos0 = { { 0.0, 0.0, 0.0 } }; + spec.stepIPP = { { 2.0 * n[0], 2.0 * n[1], 2.0 * n[2] } }; + spec.numSlices = 4; + + auto files = writeSyntheticSeries(spec); + for (bool reverse : { false, true }) + { + std::vector ord = files; + if (reverse) + { + std::reverse(ord.begin(), ord.end()); + } + auto img = readSeries(ord); + const auto & d = img->GetDirection(); + std::array c0 = { { d[0][0], d[1][0], d[2][0] } }; + std::array c1 = { { d[0][1], d[1][1], d[2][1] } }; + std::array c2 = { { d[0][2], d[1][2], d[2][2] } }; + auto c0xc1 = cross(c0, c1); + // c0×c1 must be parallel to c2 (either +c2 or −c2 depending on handedness). + const double sign = (c0xc1[0] * c2[0] + c0xc1[1] * c2[1] + c0xc1[2] * c2[2]) >= 0 ? 1.0 : -1.0; + EXPECT_NEAR(c0xc1[0], sign * c2[0], kUnitTol) << "reverse=" << reverse; + EXPECT_NEAR(c0xc1[1], sign * c2[1], kUnitTol) << "reverse=" << reverse; + EXPECT_NEAR(c0xc1[2], sign * c2[2], kUnitTol) << "reverse=" << reverse; + } +} + +// --------------------------------------------------------------------------- +// 9. ForceOrthogonalDirectionOff: the legacy code path must remain unchanged +// (this test pins behavior — it should NOT be affected by the patch). +// --------------------------------------------------------------------------- +TEST_F(GDCMSeriesDirection, ForceOrthogonalDirectionOff_UsesDirectIPPDifference) +{ + SeriesSpec spec; + spec.dir = makeCaseDir("force_ortho_off"); + spec.rowOrient = { { 1.0, 0.0, 0.0 } }; + spec.colOrient = { { 0.0, 1.0, 0.0 } }; + spec.pos0 = { { 0.0, 0.0, 0.0 } }; + spec.stepIPP = { { 0.0, 0.0, -2.0 } }; + spec.numSlices = 4; + + auto files = writeSyntheticSeries(spec); + auto img = readSeries(files, /*forceOrtho=*/false); + // With ForceOrthogonalDirection OFF, the code sets direction[*][N-1] = dirN/|dirN|, + // so Z component should literally equal sign(stepIPP). + EXPECT_NEAR(img->GetDirection()[2][2], -1.0, kUnitTol); + expectOrthonormal(img->GetDirection()); +} + +// --------------------------------------------------------------------------- +// 10. Symmetry: with ForceOrthogonalDirection ON vs OFF, the resulting +// direction matrices must agree on simple axis-aligned cases (the GDCM +// cross-product third column equals the IPP-derived third column up to +// sign, and the patch reconciles the sign). +// --------------------------------------------------------------------------- +TEST_F(GDCMSeriesDirection, ForceOrthogonalOnOff_AgreeOnAxisAligned) +{ + SeriesSpec spec; + spec.dir = makeCaseDir("ortho_onoff_agree"); + spec.rowOrient = { { 1.0, 0.0, 0.0 } }; + spec.colOrient = { { 0.0, 1.0, 0.0 } }; + spec.pos0 = { { 1.0, 2.0, 3.0 } }; + spec.stepIPP = { { 0.0, 0.0, -1.5 } }; + spec.numSlices = 5; + + auto files = writeSyntheticSeries(spec); + auto on = readSeries(files, true); + auto off = readSeries(files, false); + for (unsigned int r = 0; r < 3; ++r) + { + for (unsigned int c = 0; c < 3; ++c) + { + EXPECT_NEAR(on->GetDirection()[r][c], off->GetDirection()[r][c], kUnitTol) + << "drift at [" << r << "][" << c << "]"; + } + } + ITK_EXPECT_VECTOR_NEAR(on->GetSpacing(), off->GetSpacing(), kUnitTol); + ITK_EXPECT_VECTOR_NEAR(on->GetOrigin(), off->GetOrigin(), kGeomTol); +} + +// --------------------------------------------------------------------------- +// 11. Two-slice series — the minimum non-degenerate case. Spacing must be +// correctly derived; reversal must still produce physical-point invariance. +// --------------------------------------------------------------------------- +TEST_F(GDCMSeriesDirection, TwoSlices_MinimumValidCase) +{ + SeriesSpec spec; + spec.dir = makeCaseDir("two_slices"); + spec.rowOrient = { { 1.0, 0.0, 0.0 } }; + spec.colOrient = { { 0.0, 1.0, 0.0 } }; + spec.pos0 = { { 0.0, 0.0, 100.0 } }; + spec.stepIPP = { { 0.0, 0.0, -7.5 } }; + spec.numSlices = 2; + + auto files = writeSyntheticSeries(spec); + std::vector rev(files.rbegin(), files.rend()); + auto fwd = readSeries(files); + auto rv = readSeries(rev); + EXPECT_NEAR(fwd->GetSpacing()[2], 7.5, kGeomTol); + EXPECT_NEAR(rv->GetSpacing()[2], 7.5, kGeomTol); + expectPhysicalPointInvariant(fwd, rv); +} + +// --------------------------------------------------------------------------- +// 13. Many slices, irregular but monotonic step (sub-millimetre, large series). +// --------------------------------------------------------------------------- +TEST_F(GDCMSeriesDirection, ManySlices_FineSpacing) +{ + SeriesSpec spec; + spec.dir = makeCaseDir("many_slices"); + spec.rowOrient = { { 1.0, 0.0, 0.0 } }; + spec.colOrient = { { 0.0, 1.0, 0.0 } }; + spec.pos0 = { { 0.0, 0.0, 0.0 } }; + spec.stepIPP = { { 0.0, 0.0, 0.625 } }; + spec.numSlices = 64; + + auto files = writeSyntheticSeries(spec); + std::vector rev(files.rbegin(), files.rend()); + auto fwd = readSeries(files); + auto rv = readSeries(rev); + EXPECT_NEAR(fwd->GetSpacing()[2], 0.625, kGeomTol); + EXPECT_NEAR(rv->GetSpacing()[2], 0.625, kGeomTol); + expectOrthonormal(fwd->GetDirection()); + expectOrthonormal(rv->GetDirection()); + expectPhysicalPointInvariant(fwd, rv); +} + +// --------------------------------------------------------------------------- +// 14. Real ITK test data: ITKData/Input/DicomSeries (used by the original PR +// as ReadSlicesReverseOrder). Run with both default and reversed input +// orderings; physical points must be invariant. +// --------------------------------------------------------------------------- +TEST_F(GDCMSeriesDirection, RealDicomSeries_PhysicalPointInvariant) +{ + const std::string dataDir = TOSTRING(DICOM_SERIES_INPUT); + if (!itksys::SystemTools::FileIsDirectory(dataDir)) + { + GTEST_SKIP() << "DICOM_SERIES_INPUT not available: " << dataDir; + } + auto names = itk::GDCMSeriesFileNames::New(); + names->SetDirectory(dataDir); + names->SetUseSeriesDetails(true); + std::vector files = names->GetInputFileNames(); + ASSERT_GT(files.size(), 1u); + + auto fwd = readSeries(files, true, false); + auto rev = readSeries(files, true, true); + + expectOrthonormal(fwd->GetDirection()); + expectOrthonormal(rev->GetDirection()); + ITK_EXPECT_VECTOR_NEAR(fwd->GetSpacing(), rev->GetSpacing(), kGeomTol); + expectPhysicalPointInvariant(fwd, rev); +} + +// --------------------------------------------------------------------------- +// 15. ITKElastix#291 contract: for an axis-aligned axial acquisition with +// monotonically increasing Z, the third direction column is exactly +// [0, 0, 1] regardless of how the file list was sorted. +// --------------------------------------------------------------------------- +TEST_F(GDCMSeriesDirection, ITKElastix291_TimeSeriesShapedDirection) +{ + SeriesSpec spec; + spec.dir = makeCaseDir("elastix291"); + spec.rowOrient = { { 1.0, 0.0, 0.0 } }; + spec.colOrient = { { 0.0, 1.0, 0.0 } }; + spec.pos0 = { { 0.0, 0.0, 0.0 } }; + spec.stepIPP = { { 0.0, 0.0, 1.0 } }; + spec.numSlices = 8; + + auto files = writeSyntheticSeries(spec); + for (bool reverse : { false, true }) + { + std::vector ord = files; + if (reverse) + { + std::reverse(ord.begin(), ord.end()); + } + auto img = readSeries(ord); + EXPECT_NEAR(img->GetDirection()[0][2], 0.0, kUnitTol); + EXPECT_NEAR(img->GetDirection()[1][2], 0.0, kUnitTol); + // Origin sits at the first listed slice; direction Z sign must reach the last listed slice. + if (!reverse) + { + EXPECT_NEAR(img->GetDirection()[2][2], 1.0, kUnitTol); + EXPECT_NEAR(img->GetOrigin()[2], 0.0, kGeomTol); + } + else + { + EXPECT_NEAR(img->GetDirection()[2][2], -1.0, kUnitTol); + EXPECT_NEAR(img->GetOrigin()[2], 7.0, kGeomTol); + } + expectOrthonormal(img->GetDirection()); + } +} + +// --------------------------------------------------------------------------- +// 16. Origin equals the first listed slice's IPP — the patch must not move +// the origin, only adjust the direction. +// --------------------------------------------------------------------------- +TEST_F(GDCMSeriesDirection, OriginEqualsFirstListedSliceIPP) +{ + SeriesSpec spec; + spec.dir = makeCaseDir("origin_first"); + spec.rowOrient = { { 1.0, 0.0, 0.0 } }; + spec.colOrient = { { 0.0, 1.0, 0.0 } }; + spec.pos0 = { { 11.5, -22.5, 33.5 } }; + spec.stepIPP = { { 0.0, 0.0, -3.0 } }; + spec.numSlices = 4; + + auto files = writeSyntheticSeries(spec); + std::vector rev(files.rbegin(), files.rend()); + + auto fwd = readSeries(files); + EXPECT_NEAR(fwd->GetOrigin()[0], 11.5, kGeomTol); + EXPECT_NEAR(fwd->GetOrigin()[1], -22.5, kGeomTol); + EXPECT_NEAR(fwd->GetOrigin()[2], 33.5, kGeomTol); + + auto rv = readSeries(rev); + // First listed slice in 'rev' is the LAST written slice → z = 33.5 + 3*(-3) = 24.5 + EXPECT_NEAR(rv->GetOrigin()[2], 24.5, kGeomTol); +} + +// --------------------------------------------------------------------------- +// 17. Spacing is always positive, regardless of input order or orientation. +// --------------------------------------------------------------------------- +TEST_F(GDCMSeriesDirection, SpacingAlwaysPositive) +{ + for (double dz : { -1.25, +1.25 }) + { + SeriesSpec spec; + spec.dir = makeCaseDir(dz < 0 ? "spacing_pos_neg" : "spacing_pos_pos"); + spec.rowOrient = { { 1.0, 0.0, 0.0 } }; + spec.colOrient = { { 0.0, 1.0, 0.0 } }; + spec.pos0 = { { 0.0, 0.0, 0.0 } }; + spec.stepIPP = { { 0.0, 0.0, dz } }; + spec.numSlices = 4; + + auto files = writeSyntheticSeries(spec); + auto img = readSeries(files); + EXPECT_GT(img->GetSpacing()[0], 0.0); + EXPECT_GT(img->GetSpacing()[1], 0.0); + EXPECT_GT(img->GetSpacing()[2], 0.0); + EXPECT_NEAR(img->GetSpacing()[2], 1.25, kGeomTol); + } +} + +// --------------------------------------------------------------------------- +// 18. Stress: random orientations + random IPP step direction. The contract +// "physical points are invariant under input reversal" must hold for all +// of them. +// --------------------------------------------------------------------------- +TEST_F(GDCMSeriesDirection, RandomizedStress_PhysicalPointInvariance) +{ + std::mt19937 rng(20260429u); + std::uniform_real_distribution ang(-1.5, 1.5); + std::uniform_real_distribution off(-50.0, 50.0); + std::uniform_real_distribution sp(0.5, 4.0); + + for (unsigned int trial = 0; trial < 8; ++trial) + { + const double a = ang(rng), b = ang(rng); + const double ca = std::cos(a), sa = std::sin(a); + const double cb = std::cos(b), sb = std::sin(b); + std::array row = normalize3({ { ca, sa * cb, sa * sb } }); + std::array col = normalize3({ { -sa, ca * cb, ca * sb } }); + // Re-orthogonalize col against row (Gram–Schmidt) to keep DICOM happy. + double dot = row[0] * col[0] + row[1] * col[1] + row[2] * col[2]; + col = normalize3({ { col[0] - dot * row[0], col[1] - dot * row[1], col[2] - dot * row[2] } }); + auto n = normalize3(cross(row, col)); + const double step = sp(rng); + + SeriesSpec spec; + spec.dir = makeCaseDir("rand_" + std::to_string(trial)); + spec.rowOrient = row; + spec.colOrient = col; + spec.pos0 = { { off(rng), off(rng), off(rng) } }; + spec.stepIPP = { { step * n[0], step * n[1], step * n[2] } }; + spec.numSlices = 5; + + auto files = writeSyntheticSeries(spec); + std::vector rev(files.rbegin(), files.rend()); + + auto fwd = readSeries(files); + auto rv = readSeries(rev); + expectOrthonormal(fwd->GetDirection()); + expectOrthonormal(rv->GetDirection()); + expectPhysicalPointInvariant(fwd, rv); + } +}