Wrong header when writing untouched NIfTI file

[fepegar]

Description

Header of the written NIfTI file is not valid if no filters are applied to a read NIfTI image. qfac, which is encoded in pixdim[0], should always be -1 or 1, but it's 0.

This was first reported on the Discourse forum.

Steps to Reproduce

1. Read NIfTI image
2. Write it
3. qfacis not valid in the header of the written file

This SimpleITK code should reproduce the issue:

import struct
import numpy as np
import nibabel as nib
import SimpleITK as sitk

def get_qfac(image_path):
    with open(image_path, 'rb') as f:
        fmt = 8 * 'f'
        size = struct.calcsize(fmt)
        f.seek(76)
        chunk = f.read(size)
        pixdim = struct.unpack(fmt, chunk)
    qfac = pixdim[0]
    return qfac

def check_qfac(image_path):
    qfac = get_qfac(image_path)
    if qfac in (-1, 1):
        print('qfac is ok:', qfac)
    else:
        print(f'qfac is {qfac} in {image_path}')

filepath_nib = '/tmp/test_nib.nii'
filepath_itk = '/tmp/test_itk.nii'
array = np.random.rand(10,10,10)
affine = np.eye(4)

print('Written with NiBabel')
nib.Nifti1Image(array, affine).to_filename(filepath_nib)
check_qfac(filepath_nib)

print('\nGetImageFromArray, written with ITK')
im = sitk.GetImageFromArray(array)
sitk.WriteImage(im, filepath_itk)
check_qfac(filepath_itk)

print('\nRead and write with ITK')
im = sitk.ReadImage(filepath_nib)
sitk.WriteImage(im, filepath_itk)
check_qfac(filepath_itk)

Expected behavior

qfac to be -1 or 1.

Actual behavior

qfac is 0.

Written with NiBabel
qfac is ok: 1.0

GetImageFromArray, written with ITK
qfac is ok: 1.0

Read and write with ITK
qfac is 0.0 in /tmp/test_itk.nii

Reproducibility

Always.

Versions

SimpleITK 1.2.4.

[stale[bot]]

This issue has been automatically marked as stale because it has not had recent activity. Thank you for your contributions.

[extragoya]

This bug just bit me, still seems to be present. Reading and writing an untouched file nifti file should work, especially given the ubiquity of the format.

[hjmjohnson]

@extragoya Is the sform information for this image correct? There are two ways that a nifti image may represent the physical space information.

[extragoya]

Hi @hjmjohnson

Here is the metadata for the image I am loading and then just saving (based off of how sitk loads it). Anything pop out at you? It seems ok to me, right?

aux_file : bitpix : 8 cal_max : 0 cal_min : 0 datatype : 2 descrip : dim[0] : 3 dim[1] : 256 dim[2] : 256 dim[3] : 293 dim[4] : 1 dim[5] : 1 dim[6] : 1 dim[7] : 1 dim_info : 0 intent_code : 0 intent_name : intent_p1 : 0 intent_p2 : 0 intent_p3 : 0 nifti_type : 1 pixdim[0] : 0 pixdim[1] : 1.9531 pixdim[2] : 1.9531 pixdim[3] : 5 pixdim[4] : 0 pixdim[5] : 0 pixdim[6] : 0 pixdim[7] : 0 qform_code : 1 qform_code_name : NIFTI_XFORM_SCANNER_ANAT qoffset_x : 249.023 qoffset_y : -249.018 qoffset_z : -1344.5 quatern_b : -0 quatern_c : 1 quatern_d : 0 scl_inter : 0 scl_slope : 1 sform_code : 1 sform_code_name : NIFTI_XFORM_SCANNER_ANAT slice_code : 0 slice_duration : 0 slice_end : 0 slice_start : 0 srow_x : -1.9531 0 0 249.023 srow_y : -0 -1.9531 -0 -249.018 srow_z : 0 -0 5 -1344.5 toffset : 0 vox_offset : 352 xyzt_units : 2

[gdevenyi]

That file has sform > 0 which means NIFTI method 3 should be used for reading, which doesn't use qfac

[extragoya]

Thanks @gdevenyi - for some reason if I load and save the image, and then load the image in a viewer like itkSnap or mitkWorkbench, the saved image is reflected in the x-direction. So I'm trying to diagnose the issue, but obviously it's not relevant to this particular issue given the above

[blowekamp]

If you directly load and then save an image ( in SimpleITK and ITK ), the meta-data dictionary is available when writing. Where as if the image if passes through a filter, the dictionary from the input file is removed.

When writing some ImageIO's can check the meta-data dictionary and change behavior. Check the dictionary in the image to write and try deleting all items, to see if that changes the behavior you are seeing.

[extragoya]

@blowekamp Thanks, I had tried that, but no difference.

I think I've narrowed down what's happening here. In the upstream process of creating the image I was working with downstream, somehow the image was saved with a negative spacing. That initial image seems to work fine as far as viewers like ITKSnap are concerned, but if that image is loaded and saved the q_form information is altered.

The following code snippet and output should make this clear:

def print_metadata(image_path):

    im = sitk.ReadImage(image_path)
    for k in im.GetMetaDataKeys():
        v = im.GetMetaData(k)
        print(f"{k} : {v}")

def create_image(output_path):

    # make some dummy data
    im_data = np.random.random((10,10,10))

    # here a negative was provided erroneously
    spacing = (2, -2, 5)
    # direction, with a negative for y direction
    direction = (1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0)

    im = sitk.GetImageFromArray(im_data)
    im.SetSpacing(spacing) # set spacing with problematic value
    im.SetDirection(direction)
    sitk.WriteImage(im, output_path)

def load_and_save(image_path, output_path):

    im = sitk.ReadImage(image_path)
    sitk.WriteImage(im, output_path)

create_image('first.nii.gz')
load_and_save('first.nii.gz', 'second.nii.gz')
print('**********FIRST*********')
print_metadata('first.nii.gz')
print('*********SECOND*********')
print_metadata('second.nii.gz')

Which will give the following output

**********FIRST*********
ITK_FileNotes :
ITK_original_direction : [UNKNOWN_PRINT_CHARACTERISTICS]

ITK_original_spacing : [UNKNOWN_PRINT_CHARACTERISTICS]

aux_file :
bitpix : 64
cal_max : 0
cal_min : 0
datatype : 64
descrip :
dim[0] : 3
dim[1] : 10
dim[2] : 10
dim[3] : 10
dim[4] : 1
dim[5] : 1
dim[6] : 1
dim[7] : 1
dim_info : 0
intent_code : 0
intent_name :
intent_p1 : 0
intent_p2 : 0
intent_p3 : 0
nifti_type : 1
pixdim[0] : 0
pixdim[1] : 2
pixdim[2] : 2
pixdim[3] : 5
pixdim[4] : 0
pixdim[5] : 0
pixdim[6] : 0
pixdim[7] : 0
qform_code : 1
qform_code_name : NIFTI_XFORM_SCANNER_ANAT
qoffset_x : -0
qoffset_y : -0
qoffset_z : 0
quatern_b : -0
quatern_c : 1
quatern_d : 0
scl_inter : 0
scl_slope : 1
sform_code : 1
sform_code_name : NIFTI_XFORM_SCANNER_ANAT
slice_code : 0
slice_duration : 0
slice_end : 0
slice_start : 0
srow_x : -2 0 0 -0
srow_y : -0 -2 -0 -0
srow_z : 0 -0 5 0
toffset : 0
vox_offset : 352
xyzt_units : 2
*********SECOND*********
ITK_FileNotes :
ITK_original_direction : [UNKNOWN_PRINT_CHARACTERISTICS]

ITK_original_spacing : [UNKNOWN_PRINT_CHARACTERISTICS]

aux_file :
bitpix : 64
cal_max : 0
cal_min : 0
datatype : 64
descrip :
dim[0] : 3
dim[1] : 10
dim[2] : 10
dim[3] : 10
dim[4] : 1
dim[5] : 1
dim[6] : 1
dim[7] : 1
dim_info : 0
intent_code : 0
intent_name :
intent_p1 : 0
intent_p2 : 0
intent_p3 : 0
nifti_type : 1
pixdim[0] : 0
pixdim[1] : 2
pixdim[2] : 2
pixdim[3] : 5
pixdim[4] : 0
pixdim[5] : 0
pixdim[6] : 0
pixdim[7] : 0
qform_code : 1
qform_code_name : NIFTI_XFORM_SCANNER_ANAT
qoffset_x : -0
qoffset_y : -0
qoffset_z : 0
quatern_b : 0
quatern_c : 0
quatern_d : 1
scl_inter : 0
scl_slope : 1
sform_code : 1
sform_code_name : NIFTI_XFORM_SCANNER_ANAT
slice_code : 0
slice_duration : 0
slice_end : 0
slice_start : 0
srow_x : -2 0 0 -0
srow_y : 0 -2 0 -0
srow_z : 0 0 5 0
toffset : 0
vox_offset : 352
xyzt_units : 2

Notice the q_form information in "second.nii.gz" has been altered, despite the image only just being loaded and saved. I understand these should be ignored since s_form is 1, but it still seems to affect itkSnap and mitkWorkBench. as maybe they are not ignoring the q_forms?

Obviously, the source of the error is the user who created negative spacing, so that is not SimpleITK's fault. However, in such cases SimpleITK doesn't seem to throw an error or warn the user, and continues on silently as if there is no issue. As far as I can tell, there will be no issue with "first.nii.gz", but "second.nii.gz" will indeed have an issue as that image will be reflected in the x direction by the viewer softwares I mentioned above. It's a pernicious error for the user to trace given that "second.nii.gz" was generated only by saving and loading.

I guess ideally:

1. SimpleITK would throw and error or warn for those silly enough to set negative spacings.
2. If 1 is not possible, then at least SimpleITK should not alter the q_form data when loading and saving images, even should they have some erroneous values.

[blowekamp]

I guess ideally:

1. SimpleITK would throw and error or warn for those silly enough to set negative spacings.
2. If 1 is not possible, then at least SimpleITK should not alter the q_form data when loading and saving images, even should they have some erroneous values.

Great idea 😏

Here is the code that does it: https://github.com/InsightSoftwareConsortium/ITK/blame/cb2ff13ba02b8cbbfc94a40c3f15e00a89a16487/Modules/Core/Common/include/itkImageBase.hxx#L82-L93

It's available in the SimpleITK 2.2.0 release candidates. Please let us know if that give you the expected behavior.

[hjmjohnson]

Excellent! https://github.com/InsightSoftwareConsortium/ITK/commit/17a1266e8f92f552831497421e19cb09c72fd2f8

Thanks, @dzenanz for making this improvement.

[extragoya]

Haha well excellent timing on my part. That's very good to see :)

[dzenanz]

throw an error or warn for those silly enough to set negative spacings

I was going to ask didn't we do just that recently-ish 😄