Support specification of Nifti sform tolerance via environmental variable

[thewtex]

Per discussion with @effigies , it would be helpful to enable specification of a higher nifti sform tolerance via environmental variable.

There are cases where nifti headers, perhaps edited manually, that may have unintentional small shear, it would be helpful to specify an override via environmental variable a tolerance for orthogonality.

This would be used in docker images that use FSL, ANTs, etc.

Possibly related?

• https://github.com/InsightSoftwareConsortium/ITK/blob/fb61d2cadc33bf50d1f0ec5b842b70a1a7021d0e/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx#L1989-L2010

Todo:

• [ ] find a test case that demonstrates the issue
• [ ] create a patch (if needed) to enable specifying the related tolerance via environmental variable
• [ ] look into updated associated warnings to suggest setting the environmental variable if appropriate

CC @hjmjohnson

[effigies]

Just adding a note that I'm having trouble finding records of when we were hitting this. My memory was that there were small shears causing ITK to reject images, but it could also be differences between the sform column norms and pixdim[1:4], which I believe is what this particular code is handling.

Will keep my eyes out for the next one. When we've encountered it, it's generally been due to a failure to validate and update headers before passing to ANTs, so it's possible we won't trigger the case again until the next refactor leaves a hole for a messy image to slip through.

[cookpa]

I think this is already available via #4009

https://github.com/InsightSoftwareConsortium/ITK/blob/fb61d2cadc33bf50d1f0ec5b842b70a1a7021d0e/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx#L2080-L2123

This works for me:

export ITK_NIFTI_SFORM_PERMISSIVE=1
PrintHeader bad_sform.nii

[cookpa]

BTW, the expected values of ITK_NIFTI_SFORM_PERMISSIVE are true/false, on/off, yes/no. It cannot be toggled with 1/0. Is this ITK convention, or something to fix?

[cookpa]

Sorry for not reading carefully before, I see the idea here is to add a numeric tolerance rather than the boolean we currently have - this seems like a good idea.

I had a go at creating a test image. However, I noticed that the permissive mode introduces a flip.

It's very possible I'm doing something wrong here, so here's R code to create the dummy images, starting from a template

library(RNifti)

# the nifti class contains pointers to memory, so read them
# independently to be safe
mni <- readNifti('tpl-MNI152NLin2009cAsym_res-02_T1w.nii.gz')
mni_rotated <- readNifti('tpl-MNI152NLin2009cAsym_res-02_T1w.nii.gz')
mni_sheared <- readNifti('tpl-MNI152NLin2009cAsym_res-02_T1w.nii.gz')

qoffset_x <- mni$qoffset_x
qoffset_y <- mni$qoffset_y
qoffset_z <- mni$qoffset_z

# Create the translation vector using qoffset values
translation_vector <- c(qoffset_x, qoffset_y, qoffset_z)

spacing <- diag(c(2, 2, 2))

# blank out qform
qform(mni) <- structure(diag(4), code=0L)
qform(mni_rotated) <- structure(diag(4), code=0L)
qform(mni_sheared) <- structure(diag(4), code=0L)

writeNifti(mni, 'mni_sform.nii.gz')

rotation <- matrix(c(
  0.99765903, -0.06763516,  0.01009717,
  0.06808255,  0.99622887, -0.05378478,
  -0.00642135,  0.05434632,  0.9985015
), nrow = 3, byrow = TRUE)

shear_matrix <- matrix(c(
  1, 0.0001, 0.0001,
  0.0, 1, 0.0,
  0.0, 0.0, 1
), nrow = 3, byrow = TRUE)

almost_rotation <- shear_matrix %*% rotation

sform_rotated <- rbind(cbind(rotation %*% spacing, translation_vector), c(0, 0, 0, 1))

sform(mni_rotated) <- structure(sform_rotated, code=1L)

writeNifti(mni_rotated, 'mni_rotated.nii.gz')

sform_sheared <- rbind(cbind(almost_rotation %*% spacing, translation_vector), c(0, 0, 0, 1))

sform(mni_sheared) <- structure(sform_sheared, code=1L)

writeNifti(mni_sheared, 'mni_sheared.nii.gz')

The shear is small but enough to make non-permissive mode fail. ITK-SNAP 4.2.0 refuses to read mni_sheared.nii.gz, and so does ANTs PrintHeader:

PrintHeader mni_sheared.nii.gz 4
 cant read mni_sheared.nii.gz

ITK_NIFTI_SFORM_PERMISSIVE=TRUE PrintHeader mni_sheared.nii.gz 4

WARNING: In /Users/pcook/tmp/NOT_BACKED_UP/antsMasterBuild/build/ITKv5/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx, line 2111
NiftiImageIO (0x7fccec328ce0): mni_sheared.nii.gz has non-orthogonal sform

WARNING: In /Users/pcook/tmp/NOT_BACKED_UP/antsMasterBuild/build/ITKv5/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx, line 2111
NiftiImageIO (0x7fccec328440): mni_sheared.nii.gz has non-orthogonal sform

-0.997662x-0.0680327x0.00647123x0.0675826x-0.996232x-0.0543497x0.0101444x-0.0537853x0.998501

To try to visualize what happens in permissive mode, I used antsApplyTransforms. This should read the image using the permissive mode code, and then resample it into the space of mni_rotated, which does not trigger the permissive mode sform formulation

ITK_NIFTI_SFORM_PERMISSIVE=TRUE antsApplyTransforms \
  -d 3 -i mni_sheared.nii.gz -r mni_rotated.nii.gz -o mni_shear_permissive.nii.gz  \
  -t Identity --verbose --float

But loading all these into ITK-SNAP, mni_shear_permissive.nii.gz appears rotated. I figured this was because it was going back to the qform, so I removed that

library(RNifti)
permissive <- readNifti('mni_shear_permissive.nii.gz')
qform(permissive) <- structure(diag(4), code=0L)
c<- structure(diag(4), code=0L)
writeNifti(permissive, 'mni_shear_permissive_sform.nii.gz')

So I think the qform is not the issue, the sform created in permissive mode seems to be flipped somehow. Perhaps some sign indeterminancy in SVD?

itksnap -g mni_sform.nii.gz -o mni_rotated.nii.gz mni_shear_permissive_sform.nii.gz 

cc: @vfonov

[seanm]

Why an environment variable and not an API to the NIfTI reader class?

[dzenanz]

Both would be ideal. The environment variable is usable by end-users of the application, whereas the API only by the programmer.

[cookpa]

There's also the CMake option

option(ITK_NIFTI_IO_SFORM_PERMISSIVE_DEFAULT "Allow non-orthogonal rotation matrix in NIFTI sform by default" OFF)

Do we want a compile time numerical threshold as well?

[seanm]

Both would be ideal. The environment variable is usable by end-users of the application, whereas the API only by the programmer.

Well, the programmer can use the API to provide some control to the end user by some user interface.

Environment variables are global state, and like a secret API, that application developers might not even know about, they could result in surprising changes in end results. With NIfTI, we are talking medical-adjacent software, which can often be highly regulated, you might explicitly want for such behaviour changes to not be possible except for via the UI you have designed and tested.

[cookpa]

One other thing I think we may want to correct is that in permissive mode, a non-orthonormal sform is always used even if a valid qform exists

https://github.com/InsightSoftwareConsortium/ITK/blob/9a5dade7f7e13fcef307f8224a01c2ce2b0ab13b/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx#L2032-L2123

IMO, we should only use a non-orthonormal sform if no qform exists.

[effigies]

By the NIfTI standard, if sform_code > 0, then the sform should be used, and the qform only consulted if sform_code == 0 && qform_code > 0. The preference of ITK for qform over sform has been a pain point of interoperability for years. IMO it would be a shame to intentionally move back in that direction.

[cookpa]

I made a flow chart of what I think the current code does.

The qform parameters are used to make a qto 4x4 matrix similar to the sform sto. These are then checked element-by-element for similarity. But I don't think this check is of much consequence because there's a logical OR here (I think it should be an AND)

https://github.com/InsightSoftwareConsortium/ITK/blob/3b04ad1b377c6d5d77a2c1f454c363eab8197340/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx#L1976

Next, the sto and qto are decomposed with SVD. The sto matrix is checked for skew (qto by design can't have any) and the axes are checked for similarity between qto and sto.

I think the "sform decomposable without skew" threshold is the best target for a customizable threshold. Passing this test will guarantee the use of the sform in cases where

1. sform_code == NIFTI_XFORM_SCANNER_ANAT

2. sform_code != NIFTI_XFORM_UNKNOWN && qform_code == NIFTI_XFORM_UNKNOWN. This case leads to an unhandled exception if the skew test fails.

BTW, even when sform is used, spacing is still derived from the pixdims. Rotation and translation come from the sform. If other scaling is present in the sform, a warning is printed, but that scaling is disregarded. This is contrary to the NIFTI standard, which says pixdims should be ignored when using the sform.

[hjmjohnson]

PR request with sample desired implementation would be greatly appreciated!