Itk: regression: 8d093d56 changed image spacing behaviour

[seanm]

Description

Testing recent ITK 5 against my own codebase, one of my tests now fails. I bisected it to:

commit 8d093d5600113f94b0d635184efc14161035ce81
Author: Dženan Zukić <[email protected]>
Date:   Fri Aug 3 10:19:33 2018 -0400

    ENH: compute spacing using n-1 instead 2-1 in ImageSeriesReader

    Use first and last slice for computing the spacing,
    instead of second and first - it is more numerically stable.
    Also simplify GenerateOutputInformation method somewhat.

    Change-Id: Id8076ae195e6d145d8a9051a94276ded172904d5

Steps to Reproduce

Reduced code as follows:

#include <cstdint>
#include <string>
#include <vector>

#include "itkConfigure.h"
#include "itkImageSeriesReader.h"
#include "itkGDCMImageIO.h"
#include "itkGDCMSeriesFileNames.h"

int main ()
{
    std::vector<std::string> fileNames = {
        "Unit Tests Files/DICOM/DeliberatelyMissingSlices/000001.IMA",
        "Unit Tests Files/DICOM/DeliberatelyMissingSlices/000002.IMA",
        "Unit Tests Files/DICOM/DeliberatelyMissingSlices/000004.IMA",
        "Unit Tests Files/DICOM/DeliberatelyMissingSlices/000005.IMA",
    };

    typedef itk::Image<uint16_t, 3> RR3DImageITKType;
    itk::GDCMImageIO::Pointer gdcmIO = itk::GDCMImageIO::New();
    typename itk::ImageSeriesReader<RR3DImageITKType>::Pointer reader = itk::ImageSeriesReader<RR3DImageITKType>::New();
    reader->SetFileNames(fileNames);
    reader->SetImageIO(gdcmIO);
    reader->Update();

    auto image = reader->GetOutput();
    auto spacing = image->GetSpacing();
    double x = spacing[0];
    double y = spacing[1];
    double z = spacing[2];

    bool isAsExpected = (z == 1.5);

    std::cout << "spacing is x: " << x << " y: " << y << " z: " << z << "\n";

    return isAsExpected ? 0 : 1;
}

I can't share the DICOM file here, but maybe could privately. I suspect it would occur with any though. Notice slice 3 is missing. That's the point of the test.

Expected behavior

Previous to 8d093d56 the z spacing was 1.5.

Actual behavior

Starting with 8d093d56 the z spacing is 2.0.

[seanm]

@dzenanz thoughts?

[issakomi]

IMHO, if the image is not uniform both z-spacing values (1.5 or 2.0) are wrong/doesn't matter? I guess it is ("missing slices"). The problem is - ITK IO doesn't examine DICOM series. To handle image as 3D uniform image, it should be examined first, every slice (in fact rather easy) and recognized as volume only after. If image is not a volume some variable (currently doesn't exist ), e.g. "uniform" or "volume" could be set to false to prevent/warn about wrong 3D views.

Fig. 1
Fig. 2
Fig. 3

[dzenanz]

As @issakomi said, both are wrong. What you were testing is whether the broken input was interpreted in the way you expected.

Slicer recently added some support for handling missing slices and similar things.

[seanm]

If both 1.5 and 2.0 are "wrong/doesn't matter", then going back to the answer ITK has been giving for years should be ok too, yes?

The docs here https://itk.org/Doxygen/html/classitk_1_1Image.html say "The spacing is the geometric distance between image samples." Given that definition, 1.5 is a much more correct answer than 2.0 (for this dataset).

It would help to know what numerical stability issue the commit messages refers to. Perhaps there is a way we can restore the old answer and also have improved numerical stability.

[dzenanz]

Improved numerical stability is just a minor side-effect. With oblique volumes, it is a matter of hypothenuse or cathetus being the spacing. It is demonstrated here.

If the old behavior in case of missing slices is important to you, a PR would be welcome.

[seanm]

Ah. Well, the commit message's only stated rational was numeric stability. :) Which implied to me it could be somewhat revertable...

But I guess going back to using the second and first slice would break supporting oblique volumes? Is there a way to know the volume is oblique? Perhaps then we could use the new code then, and use the old way (second and first slice) otherwise?

[dzenanz]

We could probably deduce the obliqueness from first two slices too, but that's where the numerical stability part comes into play.

[issakomi]

Is there a way to know the volume is oblique?

Image is uniform if

• all slices have same "Image Orientation (Patient)"
• distances between every pair of slices are equal
• z-direction computed from "Image Orientation (Patient)" and .real z-direction are equal

z-direction from "Image Orientation (Patient)" is right-handed cross product of x- and y-directions
real z-direction = _normalized(origin of last slice - origin of 1st slice)_, slices must be sorted by "Image position (Patient)" before

Serious DICOM viewers with 3D support do it.
Screenshot
AFAIK Slicer does too.

Sorry, my first experience with medical imaging was related to real stereotactic brain surgery, you can imagine what happens in case of mistakes.

Edit: for "paranoid" check also possible to compare expected and real z-direction for every pair of slices, but it is probably not required.

[phcerdan]

ITK has no support for different spacing between slices. Also, it does not have a clean, direct support for sheared images. But, the direction matrix and the spacing can be cleverly "tweaked" to include the shear, which is what Dzenan patch does. However this makes the Direction to lose its meaning, it is not an orthonormal director cosines matrix anymore.

One fix would be to include a shear matrix to the metadata. Or even better, and just support all the cases involving linear and affine transforms, store the final transformation matrix.
That final transformation matrix, is computed every time in ITK based on Spacing and Direction:

From the Software Guide:

If we store the transform matrix A, which will be constructed based on metadata (Spacing, Direction, Shear, etc...) I think we will gain on flexibility and will cover almost every case.

For the non-uniforms cases you are discussing, could we create slices with interpolated values to make the volume uniform? https://discourse.itk.org/t/resampling-to-isotropic-signal-processing-theory/1403
I guess that is what Slicer does.