SimpleITK writing NIfTI with invalid header

Hello,

I’m using SimpleITK 1.2.4 to read and write NIfTI images. NIfTI docs say that it “should not occur” that qfac is 0. qfac is the first value of the pixdim array in the header and is expected to be -1 or 1.

When I create an image using GetImageFromArray, the written header is fine. If I read and write and image with SimpleITK, the header is incorrect. Here’s some code to show what I mean, and the output:

import gzip
import struct
from pathlib import Path
from tempfile import NamedTemporaryFile
import numpy as np
import nibabel as nib
import SimpleITK as sitk

def get_qfac(image_path):
    image_path = Path(image_path)
    with NamedTemporaryFile(suffix='.nii') as f:
        if image_path.suffix == '.gz':
            with gzip.open(image_path, 'rb') as f_in:
                with open(f.name, 'wb') as f_out:
                    shutil.copyfileobj(f_in, f_out)
            image_path = f.name
        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}')

def main():
    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)


if __name__ == "__main__":
    main()

Output:

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

Do you know what could be happening?

When you first read an image in ITK, the Image may contain additional meta data tags from the reader. I did a test with NifTi and there were quite a lot of NifTi tags:

In [1]: import SimpleITK as sitk                                                                                                                                                                                           

In [2]: img = sitk.Image([10,9], sitk.sitkUInt8)                                                                                                                                                                           

In [3]: sitk.WriteImage(img, "foo.nii")                                                                                                                                                                                    

In [4]: rimg = sitk.ReadImage("foo.nii")                                                                                                                                                                                   
...
In [6]: for k in rimg.GetMetaDataKeys(): 
   ...:     v = rimg.GetMetaData(k) 
   ...:     print("({0}) = = \"{1}\"".format(k,v)) 
   ...:                                                                                                                                                                                                                    
(ITK_FileNotes) = = ""
(aux_file) = = ""
(bitpix) = = "8"
(cal_max) = = "0"
(cal_min) = = "0"
(datatype) = = "2"
(descrip) = = ""
(dim[0]) = = "2"
(dim[1]) = = "10"
(dim[2]) = = "9"
(dim[3]) = = "1"
(dim[4]) = = "1"
(dim[5]) = = "1"
(dim[6]) = = "1"
(dim[7]) = = "1"
(dim_info) = = ""
(intent_code) = = "0"
(intent_name) = = ""
(intent_p1) = = "0"
(intent_p2) = = "0"
(intent_p3) = = "0"
(pixdim[0]) = = "1"
(pixdim[1]) = = "1"
(pixdim[2]) = = "1"
(pixdim[3]) = = "1"
(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) = = "0"
(sform_code_name) = = "NIFTI_XFORM_UNKNOWN"
(slice_code) = = ""
(slice_duration) = = "0"
(slice_end) = = "0"
(slice_start) = = "0"
(srow_x) = = "0 0 0 0"
(srow_y) = = "0 0 0 0"
(srow_z) = = "0 0 0 0"
(toffset) = = "0"
(vox_offset) = = "352"
(xyzt_units) = = ""

I’m guessing that the itk::NifTiImageIO reader is trying to do something “smart” with the NifTi tags to preserve them. Clearly there is a bug here.

With ITK, when an image if filtered the meta-data tags are not propagated, the this issue is not frequently encountered.

You can call EraseMetaData to remove the extraneous/redundant tags.

2 Likes

Thanks for your reply, @blowekamp.

Adding im.EraseMetaData('pixdim[0]') before my calls to sitk.WriteImage in my previous example didn’t help.

Try deleting all MetaData elements. The default behavior of an ITK filter is not to propagate the tags so “I’m += 0” would be the easiest was to accomplish this.

Thanks, I’ll do that. Should I create an issue on ITK/SimpleITK on GitHub?

Instead of nifti, can you use a general-purpose 3D image file format, such as nrrd or metaimage?

Nifti makes sense for neuroimaging, but for all other fields, it is so inconvenient - there are many complications due to that directions can be defined using several different approaches, header is not human-readable and not easy to extend, there are lots of irrelevant (neuroimaging specific) fields but nothing for any other specialty, etc.

1 Like

Hi Andras,

Definitely going to switch to NRRD soon, I’ve had enough issues with NIfTI over the last years.

If there is a bug it would be in the itk:NifTiImageIO class, so that is an ITK bug.

I was able to reproduce the same behavior. However, if all NifTi metadata tags are removed before rewriting the file then the error does not occur. As @lassoan says, the NifTi spec is a bit complicated and it’s not clear what the valid states are for the redundant data.

Yes, but this is a hack.
I’ll report the bug on the ITK tracker. Thanks for your help, @blowekamp.

There is something funny going on here that hopefully can be addressed. However, this is not as much of a “hack” as one might think. It is best practice to inspect the meta-data dictionary before writing an itk::Image. Since the meta-data dictionary is not propagate between filter generally this step is omitted. For this case the reader provides every detail of the NifTi header in the meta-data dictionary. These meta-data fields should not be written out as meta-data, but the fields are represented by then itk::Image’s geometric information which gets written in the NifTi header.

If you would like code that directly does the operations consider:

def strip_metadata(img):
  for k in img.GetMetaDataKeys():
    img.EraseMetaData(k)