Skip to content

Commit 84e2a15

Browse files
blowekampblowekamp
committed
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.
1 parent 0291f7c commit 84e2a15

File tree

3 files changed

+365
-3
lines changed

3 files changed

+365
-3
lines changed

‎Modules/IO/GDCM/test/CMakeLists.txt‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -448,3 +448,23 @@ addcompliancetest(raw-YBR_FULL_422)
448448
addcompliancetest(RLE-RGB)
449449
addcompliancetest(HTJ2K-YBR_ICT Lily_full.mha)
450450
addcompliancetest(HTJ2K-YBR_RCT Lily_full.mha)
451+
452+
set(ITKGDCMImageIOGTests itkGDCMImageIOGTest.cxx)
453+
creategoogletestdriver(ITKGDCMImageIO "${ITKIOGDCM-Test_LIBRARIES}" "${ITKGDCMImageIOGTests}")
454+
455+
target_compile_definitions(
456+
ITKGDCMImageIOGTestDriver
457+
PRIVATE
458+
"-DITK_TEST_OUTPUT_DIR=${ITK_TEST_OUTPUT_DIR}"
459+
)
460+
461+
ExternalData_Expand_Arguments(
462+
ITKData
463+
_DICOM_SERIES_INPUT
464+
"DATA{${ITK_DATA_ROOT}/Input/DicomSeries/,REGEX:Image[0-9]+.dcm}"
465+
)
466+
target_compile_definitions(
467+
ITKGDCMImageIOGTestDriver
468+
PRIVATE
469+
"-DDICOM_SERIES_INPUT=${_DICOM_SERIES_INPUT}"
470+
)

‎Modules/IO/GDCM/test/itkGDCMImageIOGTest.cxx‎

Lines changed: 334 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,334 @@
1+
/*=========================================================================
2+
*
3+
* Copyright NumFOCUS
4+
*
5+
* Licensed under the Apache License, Version 2.0 (the "License");
6+
* you may not use this file except in compliance with the License.
7+
* You may obtain a copy of the License at
8+
*
9+
* https://www.apache.org/licenses/LICENSE-2.0.txt
10+
*
11+
* Unless required by applicable law or agreed to in writing, software
12+
* distributed under the License is distributed on an "AS IS" BASIS,
13+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14+
* See the License for the specific language governing permissions and
15+
* limitations under the License.
16+
*
17+
*=========================================================================*/
18+
19+
#include "itkImage.h"
20+
#include "itkImageFileWriter.h"
21+
#include "itkGDCMImageIO.h"
22+
#include "itkGDCMSeriesFileNames.h"
23+
#include "itkMetaDataObject.h"
24+
#include "itkGTest.h"
25+
#include "itksys/SystemTools.hxx"
26+
#include "itkImageSeriesReader.h"
27+
#include <string>
28+
#include <vector>
29+
30+
31+
#define _STRING(s) #s
32+
#define TOSTRING(s) std::string(_STRING(s))
33+
34+
namespace
35+
{
36+
37+
struct ITKGDCMImageIO : public ::testing::Test
38+
{
39+
void
40+
SetUp() override
41+
{
42+
itksys::SystemTools::MakeDirectory(m_TempDir);
43+
}
44+
45+
void
46+
TearDown() override
47+
{
48+
itksys::SystemTools::RemoveADirectory(m_TempDir);
49+
}
50+
51+
const std::string m_TempDir{ TOSTRING(ITK_TEST_OUTPUT_DIR) + "/ITKGDCMImageIO" };
52+
const std::string m_DicomSeriesInput{ TOSTRING(DICOM_SERIES_INPUT) };
53+
};
54+
55+
class ITKGDCMSeriesTestData : public ITKGDCMImageIO
56+
{
57+
public:
58+
using PixelType = uint16_t;
59+
static constexpr unsigned int Dimension = 2;
60+
using ImageType = itk::Image<PixelType, Dimension>;
61+
using WriterType = itk::ImageFileWriter<ImageType>;
62+
63+
void
64+
SetUp() override
65+
{
66+
itksys::SystemTools::MakeDirectory(m_TempDir);
67+
CreateTestDicomSeries();
68+
}
69+
70+
void
71+
TearDown() override
72+
{
73+
itksys::SystemTools::RemoveADirectory(m_TempDir);
74+
}
75+
76+
protected:
77+
const std::string m_TempDir{ TOSTRING(ITK_TEST_OUTPUT_DIR) + "/ITKGDCMSeriesTestData" };
78+
std::vector<std::string> m_DicomFiles;
79+
80+
81+
private:
82+
void
83+
CreateTestDicomSeries()
84+
{
85+
// DICOM meta-data values from the report
86+
const std::vector<std::string> positions = {
87+
"-216.500\\-216.500\\70.000", // slice 1 (top)
88+
"-216.500\\-216.500\\-187.500", // slice 2 (middle)
89+
"-216.500\\-216.500\\-445.000" // slice 3 (bottom)
90+
};
91+
92+
const std::string orientation = "1.000000\\0.000000\\0.000000\\0.000000\\1.000000\\0.000000";
93+
94+
// Create a 2x2 image for each slice (equivalent to 3D [2,2,3] sliced)
95+
ImageType::SizeType size;
96+
size[0] = 2;
97+
size[1] = 2;
98+
99+
ImageType::IndexType start;
100+
start.Fill(0);
101+
102+
ImageType::RegionType region(start, size);
103+
104+
auto gdcmIO = itk::GDCMImageIO::New();
105+
gdcmIO->KeepOriginalUIDOn();
106+
auto writer = WriterType::New();
107+
writer->SetImageIO(gdcmIO);
108+
109+
// Write each slice as a DICOM file with appropriate tags
110+
for (size_t i = 0; i < positions.size(); ++i)
111+
{
112+
auto image = ImageType::New();
113+
image->SetRegions(region);
114+
image->Allocate();
115+
image->FillBuffer(100); // Just to have nonzero pixel values
116+
117+
// Get the metadata dictionary
118+
auto & dict = image->GetMetaDataDictionary();
119+
120+
// Set required DICOM tags
121+
itk::EncapsulateMetaData<std::string>(dict, "0020|0032", positions[i]); // ImagePositionPatient
122+
itk::EncapsulateMetaData<std::string>(dict, "0020|0037", orientation); // ImageOrientationPatient
123+
itk::EncapsulateMetaData<std::string>(dict, "0008|0060", "CT"); // Modality
124+
itk::EncapsulateMetaData<std::string>(dict, "0020|0013", std::to_string(i + 1)); // InstanceNumber
125+
itk::EncapsulateMetaData<std::string>(dict, "0010|0010", "Test^Patient"); // PatientName
126+
itk::EncapsulateMetaData<std::string>(
127+
dict, "0020|000e", "1.2.3.4.5.6.7.8"); // SeriesInstanceUID (same for all slices)
128+
129+
// Additional required DICOM tags for proper series
130+
itk::EncapsulateMetaData<std::string>(
131+
dict, "0008|0016", "1.2.840.10008.5.1.4.1.1.2"); // SOPClassUID (CT Image Storage)
132+
itk::EncapsulateMetaData<std::string>(
133+
dict, "0008|0018", "1.2.3.4.5.6.7.8.9." + std::to_string(i + 1)); // SOPInstanceUID
134+
itk::EncapsulateMetaData<std::string>(dict, "0020|000d", "1.2.3.4.5.6.7.8"); // StudyInstanceUID
135+
itk::EncapsulateMetaData<std::string>(dict, "0010|0020", "12345"); // PatientID
136+
itk::EncapsulateMetaData<std::string>(dict, "0008|0020", "20240101"); // StudyDate
137+
itk::EncapsulateMetaData<std::string>(dict, "0008|0030", "120000"); // StudyTime
138+
139+
const std::string filename = m_TempDir + "/slice_" + std::to_string(i + 1) + ".dcm";
140+
writer->SetFileName(filename);
141+
writer->SetInput(image);
142+
writer->Update();
143+
144+
m_DicomFiles.push_back(filename);
145+
}
146+
}
147+
};
148+
149+
} // namespace
150+
151+
TEST_F(ITKGDCMSeriesTestData, ReadSlicesReverseOrder)
152+
{
153+
// This test image series has non-unit meta-data:
154+
// spacing: [0.859375, 0.85939, 1.60016]
155+
// origin:[-112, -21.688, 126.894]
156+
// direction:
157+
// [1 0 0,
158+
// 0 0.466651 0.884442,
159+
// 0 -0.884442 0.466651]
160+
using PixelType = uint16_t;
161+
constexpr unsigned int Dimension = 3;
162+
using ImageType = itk::Image<PixelType, Dimension>;
163+
164+
using NamesGeneratorType = itk::GDCMSeriesFileNames;
165+
auto namesGenerator = NamesGeneratorType::New();
166+
namesGenerator->SetDirectory(m_DicomSeriesInput);
167+
namesGenerator->SetUseSeriesDetails(true);
168+
std::vector<std::string> fileNames = namesGenerator->GetInputFileNames();
169+
170+
using SeriesReaderType = itk::ImageSeriesReader<ImageType>;
171+
auto seriesReader = SeriesReaderType::New();
172+
seriesReader->SetFileNames(fileNames);
173+
auto gdcmIO = itk::GDCMImageIO::New();
174+
seriesReader->SetImageIO(gdcmIO);
175+
ASSERT_NO_THROW(seriesReader->UpdateLargestPossibleRegion());
176+
177+
ImageType::Pointer outputImage = seriesReader->GetOutput();
178+
outputImage->DisconnectPipeline();
179+
180+
seriesReader->ReverseOrderOn();
181+
ASSERT_NO_THROW(seriesReader->UpdateLargestPossibleRegion());
182+
ImageType::Pointer reversedOutputImage = seriesReader->GetOutput();
183+
reversedOutputImage->DisconnectPipeline();
184+
185+
std::cout << "baseline direction: " << outputImage->GetDirection() << std::endl;
186+
std::cout << "reversed direction: " << reversedOutputImage->GetDirection() << std::endl;
187+
EXPECT_EQ(outputImage->GetLargestPossibleRegion().GetSize(),
188+
reversedOutputImage->GetLargestPossibleRegion().GetSize());
189+
ITK_EXPECT_VECTOR_NEAR(outputImage->GetSpacing(), reversedOutputImage->GetSpacing(), 1e-6);
190+
EXPECT_NE(outputImage->GetOrigin(), reversedOutputImage->GetOrigin());
191+
192+
// calculate the index at the middle of the image
193+
ImageType::IndexType middleIndex;
194+
for (unsigned int d = 0; d < Dimension; ++d)
195+
{
196+
middleIndex[d] =
197+
outputImage->GetLargestPossibleRegion().GetIndex()[d] + outputImage->GetLargestPossibleRegion().GetSize()[d] / 2;
198+
}
199+
200+
const std::vector<ImageType::IndexType> testIndices = {
201+
{ { 0, 0, 0 } }, { { 1, 1, 1 } }, { { 2, 2, 2 } }, middleIndex
202+
};
203+
204+
// test that the reversed image has the same pixel values at the same physical location
205+
for (const auto & idx : testIndices)
206+
{
207+
ImageType::PointType point;
208+
outputImage->TransformIndexToPhysicalPoint(idx, point);
209+
auto reverseIdx = reversedOutputImage->TransformPhysicalPointToIndex(point);
210+
211+
std::cout << "Testing idx: " << idx << " reverseIdx: " << reverseIdx << std::endl;
212+
ASSERT_TRUE(reversedOutputImage->GetLargestPossibleRegion().IsInside(reverseIdx));
213+
EXPECT_EQ(outputImage->GetPixel(idx), reversedOutputImage->GetPixel(reverseIdx));
214+
}
215+
}
216+
217+
TEST_F(ITKGDCMSeriesTestData, CreateAndReadTestSeries)
218+
{
219+
// Verify that the DICOM files were created
220+
ASSERT_EQ(m_DicomFiles.size(), 3);
221+
222+
for (const auto & filename : m_DicomFiles)
223+
{
224+
ASSERT_TRUE(itksys::SystemTools::FileExists(filename));
225+
}
226+
227+
// Read the series using GDCMSeriesFileNames
228+
using NamesGeneratorType = itk::GDCMSeriesFileNames;
229+
auto namesGenerator = NamesGeneratorType::New();
230+
namesGenerator->SetDirectory(m_TempDir);
231+
namesGenerator->SetUseSeriesDetails(true);
232+
233+
std::vector<std::string> fileNames = namesGenerator->GetInputFileNames();
234+
ASSERT_EQ(fileNames.size(), 3);
235+
236+
// Read the series
237+
using ImageType3D = itk::Image<uint16_t, 3>;
238+
using SeriesReaderType = itk::ImageSeriesReader<ImageType3D>;
239+
auto seriesReader = SeriesReaderType::New();
240+
seriesReader->SetFileNames(fileNames);
241+
242+
auto gdcmIO = itk::GDCMImageIO::New();
243+
seriesReader->SetImageIO(gdcmIO);
244+
245+
ASSERT_NO_THROW(seriesReader->UpdateLargestPossibleRegion());
246+
247+
ImageType3D::Pointer image = seriesReader->GetOutput();
248+
249+
// Verify image properties
250+
ImageType3D::SizeType expectedSize = { { 2, 2, 3 } };
251+
EXPECT_EQ(image->GetLargestPossibleRegion().GetSize(), expectedSize);
252+
253+
// Verify pixel values (should be 100)
254+
EXPECT_EQ(image->GetPixel({ 0, 0, 0 }), 100);
255+
EXPECT_EQ(image->GetPixel({ 1, 1, 1 }), 100);
256+
}
257+
258+
TEST_F(ITKGDCMSeriesTestData, ReadSeriesTopToBottom)
259+
{
260+
// Read in top-to-bottom order (files ordered by ImagePositionPatient Z coordinate)
261+
using ImageType3D = itk::Image<uint16_t, 3>;
262+
using SeriesReaderType = itk::ImageSeriesReader<ImageType3D>;
263+
264+
// Get file list in top-to-bottom order
265+
const std::vector<std::string> & filesTopToBottom = m_DicomFiles;
266+
267+
auto seriesReader = SeriesReaderType::New();
268+
auto gdcmIO = itk::GDCMImageIO::New();
269+
seriesReader->SetImageIO(gdcmIO);
270+
seriesReader->SetFileNames(filesTopToBottom);
271+
seriesReader->ForceOrthogonalDirectionOn(); // explicitly set default
272+
273+
ASSERT_NO_THROW(seriesReader->UpdateLargestPossibleRegion());
274+
ImageType3D::Pointer image1 = seriesReader->GetOutput();
275+
276+
std::cout << "Top-to-bottom order:" << std::endl;
277+
std::cout << " Origin: " << image1->GetOrigin() << std::endl;
278+
std::cout << " Direction: " << image1->GetDirection() << std::endl;
279+
std::cout << " Spacing: " << image1->GetSpacing() << std::endl;
280+
281+
// Verify image properties
282+
ImageType3D::SizeType expectedSize = { { 2, 2, 3 } };
283+
EXPECT_EQ(image1->GetLargestPossibleRegion().GetSize(), expectedSize);
284+
285+
// The origin should be at the position of the first slice (top slice)
286+
ImageType3D::PointType expectedOrigin{ { -216.500, -216.500, 70.000 } }; // Z position of slice 1 (top)
287+
288+
ITK_EXPECT_VECTOR_NEAR(image1->GetOrigin(), expectedOrigin, 1e-3);
289+
290+
// Z spacing should be positive
291+
EXPECT_GT(image1->GetSpacing()[2], 0.0);
292+
// but the direction should have a negative Z component
293+
EXPECT_LT(image1->GetDirection()[2][2], 0.0);
294+
}
295+
296+
TEST_F(ITKGDCMSeriesTestData, ReadSeriesBottomToTop)
297+
{
298+
// Read in bottom-to-top order (files ordered by ImagePositionPatient Z coordinate)
299+
using ImageType3D = itk::Image<uint16_t, 3>;
300+
using SeriesReaderType = itk::ImageSeriesReader<ImageType3D>;
301+
302+
// Get file list in bottom-to-top order
303+
std::vector<std::string> filesBottomToTop(m_DicomFiles.rbegin(), m_DicomFiles.rend());
304+
305+
auto seriesReader = SeriesReaderType::New();
306+
auto gdcmIO = itk::GDCMImageIO::New();
307+
seriesReader->SetImageIO(gdcmIO);
308+
seriesReader->SetFileNames(filesBottomToTop);
309+
seriesReader->ForceOrthogonalDirectionOn(); // explicitly set default
310+
311+
ASSERT_NO_THROW(seriesReader->UpdateLargestPossibleRegion());
312+
ImageType3D::Pointer image2 = seriesReader->GetOutput();
313+
314+
std::cout << "Bottom-to-top order:" << std::endl;
315+
std::cout << " Origin: " << image2->GetOrigin() << std::endl;
316+
std::cout << " Direction: " << image2->GetDirection() << std::endl;
317+
std::cout << " Spacing: " << image2->GetSpacing() << std::endl;
318+
319+
// Verify image properties
320+
ImageType3D::SizeType expectedSize = { { 2, 2, 3 } };
321+
EXPECT_EQ(image2->GetLargestPossibleRegion().GetSize(), expectedSize);
322+
323+
// The origin should be at the position of the first slice (bottom slice)
324+
ImageType3D::PointType expectedOrigin{
325+
{ -216.500, -216.500, -445.000 }
326+
}; // X,Y from slice 1, Z position of slice 3 (bottom)
327+
328+
ITK_EXPECT_VECTOR_NEAR(image2->GetOrigin(), expectedOrigin, 1e-3);
329+
330+
// Z spacing should be positive (going from bottom to top)
331+
EXPECT_GT(image2->GetSpacing()[2], 0.0);
332+
// and the direction should have a positive Z component
333+
EXPECT_GT(image2->GetDirection()[2][2], 0.0);
334+
}

‎Modules/IO/ImageBase/include/itkImageSeriesReader.hxx‎

Lines changed: 11 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -172,8 +172,8 @@ ImageSeriesReader<TOutputImage>::GenerateOutputInformation()
172172
// Override the position if there is an ITK_ImageOrigin
173173
ExposeMetaData<Array<SpacingScalarType>>(lastReader->GetImageIO()->GetMetaDataDictionary(), key, positionN);
174174

175-
// Compute and set the inter slice spacing
176-
// and last (usually third) axis of direction
175+
// Compute and set the inter-slice spacing
176+
// and last (usually third) axis of a direction
177177
Vector<SpacingScalarType, TOutputImage::ImageDimension> dirN;
178178
for (unsigned int j = 0; j < TOutputImage::ImageDimension; ++j)
179179
{
@@ -187,9 +187,17 @@ ImageSeriesReader<TOutputImage>::GenerateOutputInformation()
187187
}
188188
else
189189
{
190+
// always positive spacing, corrected in the direction
190191
spacing[this->m_NumberOfDimensionsInImage] = dirNnorm / (numberOfFiles - 1);
191192
this->m_SpacingDefined = true;
192-
if (!m_ForceOrthogonalDirection)
193+
if (m_ForceOrthogonalDirection)
194+
{
195+
for (unsigned int j = 0; j < TOutputImage::ImageDimension; ++j)
196+
{
197+
direction[j][this->m_NumberOfDimensionsInImage] *= std::signbit(dirN[j]) ? -1.0 : 1.0;
198+
}
199+
}
200+
else
193201
{
194202
for (unsigned int j = 0; j < TOutputImage::ImageDimension; ++j)
195203
{

0 commit comments

Comments
 (0)