PR switching from qform to sform causing many downstream support issues, particularly for small voxel sizes

[gdevenyi]

The PR here, #1868

Is regularly triggering issues in various downstream projects, with users with (wrong? malformed? other?) NIFTI files which used to work with ITK-based tools and now do not.

Project-MONAI/MONAILabel#35
CoBrALab/RABIES#160
ANTsX/ANTs#1213
SimpleITK/SimpleITK#1433

I'm not sure what the solution to this is, but I'd like to propose a couple:

1. Change code to be a warning and fallthrough to qform reading
2. Provide a more informative error for the user to correct the generation of the invalid NIFTI file

[dzenanz]

People who use NIFTI a lot should weigh in. Candidates: @hjmjohnson @cookpa @vfonov @seanm @lassoan.

[cookpa]

The code does attempt to fall back to the qform if the sform cannot be used. The ANTs bug you link to was triggered by a combination of small voxels and missing qform (qform_code == 0). If I increased the voxel size, or entered a qform into the header, it worked.

I haven't had time to chase it down but I think there is some numerical test of whether to use sform that is getting confused by small voxel sizes. The data doesn't have to be oblique at all, if you have an sform of Iβ, where β is a sufficiently small number (< 0.02 seemed to trigger it for me), The sform is rejected, and the fallback to qform fails when qform_code is 0.

I think if that can be addressed, we'll be good.

EDIT: actually it's 0.02 mm, not 0.2mm.

[gdevenyi]

Thanks for the detail @cookpa. I'm going to follow up with the report I got to see what the actual state of the NIFTI header was for reference here.

[gdevenyi]

Details on the problematic scan we had:
CoBrALab/RABIES#160 (comment)

[lassoan]

I have no preference on how to change the NIFTI code or error messages to deal with this issue.

My hope is that we can build up support for OME-Zarr (or similar modern file format) in the wider image computing community and that will offer enough incentives for researchers to leave behind NIFTI and its sform/qform mess for good.

[cookpa]

I think I found it:

ITK/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx

Lines 1762 to 1767 in 02e7c44

	// Get the inverse of this matrix:
	// Make sure the matrix is invertible to begin with...
	if (fabs(vnl_determinant(mat)) <= 1e-5)
	{
	return false;
	}

This fails if the voxels are small because the determinant of the sform matrix is below the threshold (1e-5). Perhaps this could safely be made smaller.

But without eliminating the test entirely, there's always going to be a minimum spacing before this error occurs. If the entries in sform matrix are sufficiently small, qform might be the better way to get the rotation anyway. But even if the qform is used, I think very small spacing could become problematic for computing the rotation / scale matrix inverses.

Maybe users should be warned to use a more appropriate unit (eg, NIFTI_UNITS_MICRON) if the matrix determinant is below the threshold?

[cookpa]

Oh, actually other units won't help because they will get scaled to mm units

[Leengit]

@cookpa good sleuthing! The 1e-5 threshold is quite demanding; anyone know why it is there? For example, there are accurate methods for computing a matrix inverse that can handle determinants much, much smaller in absolute value than 1e-5.

[dzenanz]

PR #1868 introduced it:
34231b5#diff-fa33a46a6dc41cde5d18ee8f63cb941794512e03c6377340b8d110fe584cfdfaR1756-R1760

I guess it is fine for MRIs, which have voxels ~ 1.0 mm. @hjmjohnson is there a problem reducing this to a much lower value?

[gdevenyi]

SimpleITK/SimpleITK#1440

[dzenanz]

Gabriel, can you make a PR with your proposed changes (turn error into warning and more informative message)?

[dzenanz]

Also, can a mathematically-inclined person figure out a better threshold for determinant:

ITK/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx

Line 1764 in 02e7c44

if (fabs(vnl_determinant(mat)) <= 1e-5)

[gdevenyi]

Also, can a mathematically-inclined person figure out a better threshold for determinant:

ITK/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx

Line 1764 in 02e7c44

if (fabs(vnl_determinant(mat)) <= 1e-5)

I think solution might be to test the condition of the matrix and compare to epsilon, as suggested here, https://stackoverflow.com/a/13270760/4130016

[dzenanz]

I think solution might be to test the condition of the matrix and compare to epsilon

It would be good if this found its way into your PR 😄

[gdevenyi]

To clarify here, after working through the logic, it looks like we have two things happening together

1. ITK/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx

   Line 1764 in 02e7c44

   if (fabs(vnl_determinant(mat)) <= 1e-5)

   which is resulting in some kinds of sforms being rejected

combined with
2.

ITK/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx

Lines 1971 to 1980 in 02e7c44

	if (prefer_sform_over_qform)
	{
	return this->m_NiftiImage->sto_xyz;
	}
	else if (this->m_NiftiImage->qform_code != NIFTI_XFORM_UNKNOWN)
	{
	return this->m_NiftiImage->qto_xyz;
	}
	itkGenericExceptionMacro("ITK only supports orthonormal direction cosines. No orthonormal definition found!");

Which, without an else condition, is falling through (and not, at least according to my reading, following the NIFTI spec, which says for #define NIFTI_XFORM_UNKNOWN 0 that this should be followed:

   METHOD 1 (the "old" way, used only when qform_code = 0):
   -------------------------------------------------------
   The coordinate mapping from (i,j,k) to (x,y,z) is the ANALYZE
   7.5 way.  This is a simple scaling relationship:

     x = pixdim[1] * i
     y = pixdim[2] * j
     z = pixdim[3] * k


   No particular spatial orientation is attached to these (x,y,z)
   coordinates.  (NIFTI-1 does not have the ANALYZE 7.5 orient field,
   which is not general and is often not set properly.)  This method
   is not recommended, and is present mainly for compatibility with
   ANALYZE 7.5 files.

For a PR for this, I will add an else fallthrough, and print a warning, but the question is how should the final mapping be handled

1. Use the analyze method (qform = 0 above)
   or
2. Use the qform anyways

[gdevenyi]

Following up on this again:

ITK/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx

Line 1764 in 02e7c44

if (fabs(vnl_determinant(mat)) <= 1e-5)

Something didn't sit right with me, so I went back to my linalg textbook and a matrix A is invertible if det(A) != 0. I wonder if this test was meant to be fabs(vnl_determinant(mat)) >= 1e-5

[cookpa]

If the header has nonzero sform or qform, and they cannot be used for reasons, I'd prefer to throw an error rather than fall back to method 1. Better no data than accidentally flipped data, IMO.

I think the logic is good as is:

• Attempt to use "method 1", Analyze backwards compatibility, only if qform_code == sform_code == 0.
• Else use sform if it's defined and orthonormal direction cosines can be recovered from it
• Else try to use qform
• Else throw an error

The only thing I would add is that after attempting and failing to use sform, the user should be warned that the attempt to use sform failed, and qform will be used as a fallback. Currently, that's a silent fallback. Arguably it should be an error too as it suggests something wrong with sform, but at least we are falling back to another transform that's actually in the header, rather than assuming ANALYZE.

[cookpa]

Regarding the threshold for the determinant, I found this in the matrix inversion code:

  //: find weights below threshold tol, zero them out, and update W_ and Winverse_
  void            zero_out_absolute(double tol = 1e-8); //sqrt(machine epsilon)

https://github.com/vxl/vxl/blob/d771c601a2072afe46af8a80b8e480d0529d8f2f/core/vnl/algo/vnl_svd.h#L89-L90

They suggest 1e-8 or sqrt(machine epsilon) as a threshold for the singular values to be zeroed out, but the class looks like it defaults to a threshold of zero. If we had spacing of 1e-8, the determinant would be 1e-24.

I'll look more at the code, I'm thinking it might be better to just get the SVD decomposition directly and check it has the correct rank.

[Leengit]

Singular values (essentially eigenvalues) get multiplied together to determine the determinant, so even when 1e-8 is the right threshold for a singular value, it is probably not the correct threshold for a determinant. Also, the smallest positive single-precision floating point number that can be represented (without resorting to subnormals representations) is 1.175494e-38, and its square root is not 1e-8.

More importantly, rigorously it is the matrix conditioning number rather than its determinant that should be checked. Is that check a possibility?

[gdevenyi]

More importantly, rigorously it is the matrix conditioning number rather than its determinant that should be checked. Is that check a possibility?

Yes, if we just compute SVD directly its the ratio of the largest to smallest eigenvalues.

[cookpa]

Using the matrix condition :

https://github.com/cookpa/ITK/blob/2a6a365b8a6ee1a715254c1c3478666c4e8b8e6f/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx#L1762-L1767

This reads small (18 micron) voxels that previously failed. I'm not sure on the best way to formally test this, as any image created with this class will contain both sform and qform.

[gdevenyi]

Using the matrix condition :

Awesome, I couldn't find something that just outputted the condition number.

This fix will solve the immediate issue, I think we can leave the logic adjustment to a seperate PR?

[cookpa]

I put in a PR for the condition number stuff. I could not find a good place to insert a warning without doing more refactoring, so I left that out. IMO the only time to warn is if we can't use sform AND sform_code == NIFTI_XFORM_SCANNER_ANAT.