search

spinalcordtoolbox / spinalcordtoolbox

Comprehensive and open-source library of analysis tools for MRI of the spinal cord.
https://spinalcordtoolbox.com
GNU Lesser General Public License v3.0
200 stars 101 forks source link

Output images have their unused `dim`/`pixdim` values overwritten with `1.0` #4025

Open valosekj opened 1 year ago

valosekj commented 1 year ago

Description

When running:

sct_deepseg_sc -i sub-001_T1w.nii.gz -c t1

The output SC seg has a different pixdim that the input image:

$ sct_image -i sub-001_T1w.nii.gz -header | grep pixdim
pixdim      [-1.0, 0.351562, 0.351562, 2.5, 0.889114, 0.0, 0.0, 0.0]
$ sct_image -i sub-001_T1w_seg.nii.gz -header | grep pixdim
pixdim      [-1.0, 0.351562, 0.351563, 2.499999, 1.0, 1.0, 1.0, 1.0]

Notice the difference 2.5, 0.889114, 0.0, 0.0, 0.0 vs 2.499999, 1.0, 1.0, 1.0, 1.0.

~This difference is consequently raising errors (here and here) when using nnU-Net package.~ Update: See https://github.com/spinalcordtoolbox/spinalcordtoolbox/issues/4025#issuecomment-1445193122.

Example data available here.

joshuacwnewton commented 1 year ago

I can reproduce this issue using testing data, too: 

master|joshua@tadpole:~/repos/spinalcordtoolbox$ cd data/sct_example_data/t1
master|joshua@tadpole:~/repos/spinalcordtoolbox/data/sct_example_data/t1$ sct_deepseg_sc -i t1.nii.gz -c t1

--
Spinal Cord Toolbox (git-master-9f4d05ad2ce9bd48fd5a396310597fb5db645253)

sct_deepseg_sc -i t1.nii.gz -c t1
--

[...] 

Done! To view results, type:
fsleyes /home/joshua/repos/spinalcordtoolbox/data/sct_example_data/t1/t1.nii.gz -cm greyscale /home/joshua/repos/spinalcordtoolbox/data/sct_example_data/t1/t1_seg.nii.gz -cm red -a 70.0 &

master|joshua@tadpole:~/repos/spinalcordtoolbox/data/sct_example_data/t1$ sct_image -i t1.nii.gz -header | grep dim
dim             [3, 192, 260, 320, 1, 0, 0, 0]
pixdim          [1.0, 0.999996, 1.0, 1.0, 2.0, 0.0, 0.0, 0.0]
phase_dim       1
freq_dim        2
slice_dim       3
master|joshua@tadpole:~/repos/spinalcordtoolbox/data/sct_example_data/t1$ sct_image -i t1_seg.nii.gz -header | grep dim
dim             [3, 192, 260, 320, 1, 1, 1, 1]
pixdim          [1.0, 0.999996, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
phase_dim       0
freq_dim        0
slice_dim       0

I'll dig into the cause next.

joshuacwnewton commented 1 year ago

This behavior comes from nibabel:

See also:

The most relevant quote from that issue is:

The standard provides no interpretation for dimensions beyond dim[0] or guidance on what they should be, so there is no correct/incorrect setting. They should not be considered to have any impact on the interpretation of the data at all.

joshuacwnewton commented 1 year ago

Your input image is 3D, meaning that only the values [0:3] should be relevant; [4:7] should be "unused". (i.e. set to 0, 1, whatever...)

So, my question is: Why does the nnU-Net software care about the "unused" dim/pixdim fields? This seems overly sensitive to me.

joshuacwnewton commented 1 year ago

To fix this on SCT's end, we would have to store the original values of dim[4:7] and pixdim[4:7] ourselves in our Image objects. (But, I'm not sure if nibabel even provides access to the other values?)

Unfortunately, even if we are able to access the old dim/pixdim values, we are blocked by another issue: https://github.com/spinalcordtoolbox/spinalcordtoolbox/issues/2462#issuecomment-1396266576. Right now, SCT's Image objects mishandle dim[4:7] and pixdim[4:7], so we would have to rewrite that part of Image, which is a pretty significant effort. :sweat:

valosekj commented 1 year ago

Thanks for your swift digging into this issue, @joshuacwnewton!

So, my question is: Why does the nnU-Net software care about the "unused" dim/pixdim fields? This seems overly sensitive to me. Especially given nipy/nibabel#1034 (comment):

The standard provides no interpretation for dimensions beyond dim[0] or guidance on what they should be, so there is no correct/incorrect setting. They should not be considered to have any impact on the interpretation of the data at all.

This is what I am also wondering! I will explore the nnU-Net more deeply to find out the reason!

valosekj commented 1 year ago

After deeper digging, it seems that the nnUNet errors are not caused by pixdim how I initially thought:

Notice the difference 2.5, 0.889114, 0.0, 0.0, 0.0 vs 2.499999, 1.0, 1.0, 1.0, 1.0. This difference is consequently raising errors (here and here) when using nnU-Net package.

It also seems that the problem is not related to nibabel but to SimpleITK.

Investigation The issue comes from [`SimpleITK.Image.GetDirection()`](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1Image.html#a7bd35f108a2df5c966285405af409548) method, which is [used by nnU-Net](https://github.com/MIC-DKFZ/nnUNet/blob/d1431c223407c3a66222409cb67c19ae87b2ce6a/nnunet/preprocessing/sanity_checks.py#L45) to verify geometry between image and segmentation. See the difference in ITK `direction` between `sub-001_T1w.nii.gz` and `sub-001_T1w_seg.nii.gz` (data available [here](https://www.dropbox.com/sh/piwt2050fjbljqg/AADhZH9BnQFWcVlQFoEhBZb7a?dl=0)) using [read_itk.py](https://gist.github.com/valosekj/a03195d9060b0e164faff95102129feb#file-read_itk-py): ```console $ python ./read_itk.py sub-001_T1w.nii.gz sub-001_T1w_seg.nii.gz origin: (-90.7466812133789, 103.52926635742188, -57.17771911621094) (-90.7466812133789, 103.52926635742188, -57.17771911621094) spacing: (0.3515625, 0.3515625, 2.5) (0.3515625, 0.3515625298023224, 2.4999990463256836) direction: (0.9999999977894414, -6.478192795954982e-09, 6.649148620669932e-05, 2.0651083342197973e-05, -0.9505162149977084, -0.3106749763676334, -6.320324401990454e-05, -0.31067495075469176, 0.950516204300512) (0.9999999999453735, 1.0322302829269651e-05, 1.6441191221002595e-06, 1.032230276040639e-05, -0.9505162024305571, -0.31067497705397884, 1.6441191998447065e-06, -0.31067498903265794, 0.9505162064003994) size: (512, 512, 42) (512, 512, 42) ``` A possible fix is to set the direction from `sub-001_T1w.nii.gz` to `sub-001_T1w_seg.nii.gz` using [`SimpleITK.Image.SetDirection()`](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1Image.html#a04d6c8b3777b454b6da1ee1d3e7ed56f) by [change_itk_direction.py](https://gist.github.com/valosekj/a03195d9060b0e164faff95102129feb#file-change_itk_direction-py): ```console python ./change_itk_direction.py sub-001_T1w.nii.gz sub-001_T1w_seg.nii.gz ``` Then the direction is (almost) the same, and `nnU-Net` is satisfied: ```console $ python ./read_itk.py sub-001_T1w.nii.gz sub-001_T1w_seg_updated_dir.nii.gz origin: (-90.7466812133789, 103.52926635742188, -57.17771911621094) (-90.7466812133789, 103.52926635742188, -57.17771911621094) spacing: (0.3515625, 0.3515625, 2.5) (0.3515625, 0.3515625298023224, 2.4999990463256836) direction: (0.9999999977894414, -6.478192795954982e-09, 6.649148620669932e-05, 2.0651083342197973e-05, -0.9505162149977084, -0.3106749763676334, -6.320324401990454e-05, -0.31067495075469176, 0.950516204300512) (0.9999999977894414, -6.4774766359226164e-09, 6.649148620669932e-05, 2.0651080755190826e-05, -0.9505162169214625, -0.3106749763676334, -6.320324401990454e-05, -0.31067494486892744, 0.950516204300512) size: (512, 512, 42) (512, 512, 42) ``` Further reading: [Simple ITK Fundamental Concepts](https://simpleitk.readthedocs.io/en/master/fundamentalConcepts.html)

Thanks to my colleague @KaterinaKrejci231054 for help with debugging!

joshuacwnewton commented 1 year ago

So, to summarize this issue, there are actually 2 separate differences between the anat/seg images:

  1. A change to the dim/pixdim values. (nibabel)
  2. A discrepancy in the direction matrix. (SimpleITK)

However, only 2. is what causes the error in nnU-Net. (In other words, 1. was a red herring.)

I see! Great investigation! :tada:

Follow-up questions - What does the `direction` matrix represent? And, how is it computed? (In other words, how is it different than the various NIfTI1 image properties such as qform/sform matrices?) - What is causing the differences in the `direction` matrix between anat/seg images? - Can we somehow detect and/or avoid the change to `direction` on the SCT side of things?

EDIT: See #4025 for the conclusion.

valosekj commented 1 year ago

I have some follow-up questions:

Good questions! I don't have the answers right now. We will continue with further digging.

valosekj commented 1 year ago

Might be remotely related to https://github.com/ivadomed/model_seg_sci/issues/21. It seems that torchIO also uses SimpleITK.

valosekj commented 1 year ago

Just for the record, I encountered the same ITK direction error also for sci-zurich dataset after using sct_crop_image and sct_resample (used preprocessing script available here).

Again, I fixed the issue using change_itk_direction.py:

for sub in sub*;do cd $sub; for ses in ses*;do cd $ses/anat; python ~/code/change_itk_direction.py ${sub}_${ses}_acq-sag_T2w.nii.gz ../../../derivatives/labels/${sub}/${ses}/anat/${sub}_${ses}_acq-sag_T2w_lesion-manual.nii.gz; cd ../..; done; cd ..;echo $sub done;done

The script raised an error for sub-zh21, though:

Traceback (most recent call last):
  File "/home/GRAMES.POLYMTL.CA/p118175/code/change_itk_direction.py", line 15, in <module>
    image1  = sitk.ReadImage(sys.argv[1])
  File "/home/GRAMES.POLYMTL.CA/p118175/code/model_seg_sci/venv/lib/python3.10/site-packages/SimpleITK/extra.py", line 355, in ReadImage
    return reader.Execute()
  File "/home/GRAMES.POLYMTL.CA/p118175/code/model_seg_sci/venv/lib/python3.10/site-packages/SimpleITK/SimpleITK.py", line 8438, in Execute
    return _SimpleITK.ImageFileReader_Execute(self)
RuntimeError: Exception thrown in SimpleITK ImageFileReader_Execute: /tmp/SimpleITK-build/ITK/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx:2016:
ITK ERROR: ITK only supports orthonormal direction cosines.  No orthonormal definition found!

Thus, I removed the sub-zh21 subject entirely from the training dataset. Now nnUNetv2_plan_and_preprocess --verify_dataset_integrity is happy.

Tagging @naga-karthik to let him know.

rickymwalsh commented 1 year ago

Hi all,

I came across similar problems with discrepancies in direction matrices throwing errors when training with torchio, and also solved by using the method .SetDirection(). I dug a bit into what might cause it on SITK side and so here's my take on your questions above. I'm no expert in SITK nor C++, so take with a grain of salt :)

* What does the `direction` matrix represent? And, how is it computed? (In other words, how is it different than the various NIfTI1 image properties such as qform/sform matrices?)

It seems the direction matrix will be taken from either the qform or sform matrix based on the logic here. If sform and qform are "sufficiently" similar, then sform is used for direction matrix. Otherwise, if sform is orthonormal and sform_code==NIFTI_XFORM_SCANNER_ANAT, again sform is used. And then there are several other conditions to decide whether to use sform or qform.

Another difference is that the spacing/zooms are incorporated into the sform/qform matrices but seem not to be taken into account for the SITK direction matrix.

* What is causing the differences in the `direction` matrix between anat/seg images?

My guess is that it's because the anat image has sform_code==NIFTI_XFORM_SCANNER_ANAT and so sform will be used to get the direction matrix whereas the seg image has sform_code==NIFTI_XFORM_ALIGNED_ANAT and so the qform may be used. Besides these differences, nearly every other piece of metadata is the same in both images. 



im_anat = sitk.ReadImage('sub-001_T1w.nii.gz')
im_seg = sitk.ReadImage('sub-001_T1w_seg.nii.gz')

# Print out all metadata for both images side by side
for key in im_anat.GetMetaDataKeys():
    print(f'{key}: Anat: {im_anat.GetMetaData(key)}\t\tSeg: {im_seg.GetMetaData(key)}')
joshuacwnewton commented 1 year ago

Thank you so much for your insight!

If I'm understanding correctly, since SimpleITK's getDirection() fetches the anat image's sform (1) and the seg image's qform (2)...

If the anat image's sform (1):

> fslhd sub-001_T1w.nii.gz
sform_name  Scanner Anat
sform_code  1
sto_xyz:1   -0.351562 0.000000 -0.000166 90.746681 
sto_xyz:2   -0.000007 0.334166 0.776687 -103.529266 
sto_xyz:3   -0.000022 -0.109222 2.376290 -57.177719 
sto_xyz:4   0.000000 0.000000 0.000000 1.000000

Matched the seg image's qform (2):

> fslhd sub-001_T1w_seg.nii.gz
qform_name  Aligned Anat
qform_code  2
qto_xyz:1   -0.351562 -0.000004 -0.000004 90.746681 
qto_xyz:2   -0.000004 0.334166 0.776687 -103.529266 
qto_xyz:3   0.000001 -0.109222 2.376290 -57.177719 
qto_xyz:4   0.000000 0.000000 0.000000 1.000000

Then there would be no discrepancy between the two "direction matrices", and thus no error?

If this is the case, it seems that this is yet another issue due to discrepancies in between the qform and sform, and we've come across many:

Thankfully though, it seems that this issue doesn't stem from SCT, at least (because sct_deepseg_sc doesn't modify the sform/qform -- they are mismatched even within just the input anat image).

Actually, maybe sct_deepseg_sc is at fault here: If wasn't changing the codes of the segmented image from 1 (Scanner anat) to 2 (Aligned anat), then the sform would be fetched both times, avoiding the error. (Since the segmentation is in the same space as the anat image, shouldn't the codes be preserved?)

joshuacwnewton commented 1 year ago

Just to clean up this issue a little (since the title of this issue focuses on the dim/pixdim 1.0 issue), I've summarized the conclusions of the SITK-focused investigations in a separate issue:

And I've hidden the SITK-specific discussion as "off-topic", just so that we keep the discussion to one problem per issue.

Again, thank you so much @rickymwalsh! :hearts: