BUG: Fix case where DICOM series order has negative spacing

[blowekamp]

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

PR Checklist

• No API changes were made (or the changes have been approved)
• No major design changes were made (or the changes have been approved)
• Added test (or behavior not changed)
• Updated API documentation (or API not changed)
• Added license to new files (if any)
• Added Python wrapping to new files (if any) as described in ITK Software Guide Section 9.5
• Added ITK examples for all new major features (if any)

Refer to the ITK Software Guide for
further development details if necessary.

[blowekamp]

relates to SimpleITK/SimpleITK#2292

[dzenanz]

Looks good after a quick look.

[blowekamp]

I update the test case to use the sample MRI series in ITK test data, and use the ReverseOrder flag. The test is compare the pixel values at the physical locations match.

[dzenanz]

Maybe we could have both of those tests?

[blowekamp]

Maybe we could have both of those tests?

The DICOM tags were not correct in the prior example test case where the DICOM files were manually created. And oddly the X,Y spacing was being swapped for some reason likely related to the TAGs. Also this test case for the file has a non identify cosine matrix, which helps verifies the correct column was multiplied and not a row by -1.

[N-Dekker]

Thanks @blowekamp I'm just trying to understand the issue. It reminds me of an ITKElastix issue reported by @thewtex:

• GroupwiseRegistration notebook: ERROR: the direction cosines matrix of the fixed image is invalid! ITKElastix#291

Which appeared related to ImageFileReader::GenerateOutputInformation():

ITK/Modules/IO/ImageBase/include/itkImageFileReader.hxx

Lines 222 to 235 in b375cef

	// Spacing is expected to be greater than 0
	// If negative, flip image direction along this axis.
	// and store this information in the metadata
	for (unsigned int i = 0; i < TOutputImage::ImageDimension; ++i)
	{
	if (spacing[i] < 0)
	{
	spacing[i] = -spacing[i];
	for (unsigned int j = 0; j < TOutputImage::ImageDimension; ++j)
	{
	direction[j][i] = -direction[j][i];
	}
	}
	}

[blowekamp]

The original issue was reported in SimpleITK here:
SimpleITK/SimpleITK#2292

The cause appeared to be that the z-slice direction was negative. This was more easily reproduce by enabling the "ReverseOrder" flag.

It may very well be likely related to this change too:

ITK/Modules/IO/ImageBase/include/itkImageFileReader.hxx

Lines 222 to 235 in b375cef

	// Spacing is expected to be greater than 0
	// If negative, flip image direction along this axis.
	// and store this information in the metadata
	for (unsigned int i = 0; i < TOutputImage::ImageDimension; ++i)
	{
	if (spacing[i] < 0)
	{
	spacing[i] = -spacing[i];
	for (unsigned int j = 0; j < TOutputImage::ImageDimension; ++j)
	{
	direction[j][i] = -direction[j][i];
	}
	}
	}

It's not clear to me if the original DICOM reported in SimpleITK had negative spacing. The "ReverseOrder" reproduction of the issue, has negative inter-slice distance. There is certainly areas that need further investigation.

The principle used in the test, that for DICOM if the "reverse order" flag is set then the physical location of the voxel should be maintained, should be correct. But DICOM is complicated.

[issakomi]

For the the change related to SpacingBetweenSlices s.
#4794 (comment)
#4794 (comment)
This may lead to issues, surprisingly negative SpacingBetweenSlices it is not uncommon (negative spacing is invalid from DICOM point of view, except for the NM storage, s. next post for the example file).

DICOM doesn't force any order of slices. For classic series it is easy to sort files (ascending along slice normal, if IPP/IOP are available). This is done in itkGDCMSeriesFileNames.cxx.

For enhanced multi-frame instances is it not really done in GDCM, but if the order is descending (very very often) GDCM returns negative spacing and the code in itkImageFileReader.hxx (posted above) is triggered and so far the image is correct, but has left-handed matrix.

P.S. I don't use this stuff in my app, so it is just FYI.

[issakomi]

For enhanced multi-frame instances is it not really done in GDCM, but if the order is descending (very very often) GDCM returns negative spacing and the code in itkImageFileReader.hxx (posted above) is triggered and so far the image is correct, but has left-handed matrix.

@N-Dekker To demonstrate how GDCM returns negative spacing and the code in itkImageFileReader.hxx is triggered:

The test file can be e.g. this one, used for regular testing.

Add e.g. std::cout << sp[0] << ' ' << sp[1] << ' ' << sp[2] << std::endl; before return in gdcmImageHelper.cxx GetSpacingValueFromSequence
1 1 -1.33

or it is possible to load the file, save as MHA and look at the line ITK_original_spacing = 1 1 -1.33 (it is important to use recent ITK, there was a bug in ITK_original_spacing).

Edit: I mean at least in this particular case this code in itkImageFileReader.hxx does a good job, even it is somehow simplified approach and left-handed matrix can cause issues. BTW, there is no guarantee that frames in enhanced multi-frame are sorted at all, but in most cases they are sorted, either ascending or descending. And there are really complex cases, like multiple stacks, time frames, etc., where the naive approach can't work. To proper handle complex enhanced multi-frame: usage of the Multi-frame Dimension Module and knowledge about every frame and the ability to sort frames, if required.

Edit: BTW, I don't know, is this negative spacing here expected by GDCM or not, but surprisingly together with the code in itkImageFileReader.hxx it seems to work.

But again, I don't use all of this and don't really care about, it is just FYI.

Edit: another multi-frame test file
1.3125 1.3125 -3

Edit: S. this post (D. Clunie) with examples of intentionally randomized frames in enhanced multi-frame for testing, if interested. Not supported by GetSpacingValueFromSequence, of course.

Edit: for completeness, this file is an example there negative SpacingBetweenSlices is valid. It is only (!) applicable to NM Reconstruction Module, "The sign of the Spacing Between Slices (0018,0088) determines the direction of stacking"
ITK_original_spacing = 4.41806 4.41806 -8.83612