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?

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

  Lines 1989 to 2010 in fb61d2c

  	// Ensure that the scales are approximately the same for spacing directions
  	bool sform_scales_ok{ true };
  	constexpr float large_value_tolerance = 1e-3; // Numerical precision of sform is not very good
  	if (itk::Math::abs(this->m_NiftiImage->dx - rotation.get_column(0).magnitude()) > large_value_tolerance)
  	{
  	sform_scales_ok = false;
  	}
  	else if (itk::Math::abs(this->m_NiftiImage->dy - rotation.get_column(1).magnitude()) >
  	large_value_tolerance)
  	{
  	sform_scales_ok = false;
  	}
  	else if (itk::Math::abs(this->m_NiftiImage->dz - rotation.get_column(2).magnitude()) >
  	large_value_tolerance)
  	{
  	sform_scales_ok = false;
  	}
  	if (!sform_scales_ok)
  	{
  	itkWarningMacro(<< this->GetFileName() << " has unexpected scales in sform");
  	}
  	}

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

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

Lines 2080 to 2123 in fb61d2c

	else if (m_SFORM_Permissive) // sform is orthonormal
	{
	// Fix to deal with non-orthogonal matrixes
	// this approach uses SVD decomposition
	// maybe it is better to use nifti_make_orthog_mat44
	// to be consistent with other software?
	itkWarningMacro(<< this->GetFileName() << " has non-orthogonal sform");
	vnl_matrix_fixed<double, 3, 3> mat;
	for (int i = 0; i < 3; ++i)
	for (int j = 0; j < 3; ++j)
	{
	mat[i][j] = double{ this->m_NiftiImage->sto_xyz.m[i][j] };
	}
	vnl_svd<double> svd(mat.as_ref(), 1e-8);
	if (svd.singularities() == 0)
	{
	mat = svd.V() * svd.U().conjugate_transpose();
	mat44 _mat;
	for (int i = 0; i < 3; ++i)
	for (int j = 0; j < 3; ++j)
	{
	_mat.m[i][j] = static_cast<float>(mat[i][j]);
	}
	// preserve origin
	for (int i = 0; i < 3; ++i)
	{
	_mat.m[i][3] = this->m_NiftiImage->sto_xyz.m[i][3];
	}
	for (int i = 0; i < 4; ++i) // should be 0 0 0 1
	{
	_mat.m[3][i] = this->m_NiftiImage->sto_xyz.m[3][i];
	}
	this->m_SFORM_Corrected = true;
	return _mat;
	}
	}
	} // sform not NIFTI_XFORM_UNKNOWN

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

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

Lines 2032 to 2123 in 9a5dade

	if (sform_decomposable_without_skew)
	{
	// Use sform for direction if qform intent is unknown, and sform intent is known and orthonormal.
	if (this->m_NiftiImage->qform_code == NIFTI_XFORM_UNKNOWN &&
	this->m_NiftiImage->sform_code != NIFTI_XFORM_UNKNOWN)
	{
	prefer_sform_over_qform = true;
	}
	// Use sform if it is labeled as SCANNER_ANAT format and is orthonormal.
	else if (this->m_NiftiImage->sform_code == NIFTI_XFORM_SCANNER_ANAT)
	{
	prefer_sform_over_qform = true;
	}
	else if (this->m_NiftiImage->qform_code != NIFTI_XFORM_UNKNOWN &&
	this->m_NiftiImage->sform_code != NIFTI_XFORM_UNKNOWN)
	{
	// If sform and qform are similar, or intent is SCANNER_ANAT, prefer sform's higher numerical precision
	const bool sform_and_qform_are_very_similar = [this]() -> bool {
	const vnl_matrix_fixed<float, 4, 4> sform_as_matrix{ &(this->m_NiftiImage->sto_xyz.m[0][0]) };
	const vnl_matrix_fixed<float, 4, 4> qform_as_matrix{ &(this->m_NiftiImage->qto_xyz.m[0][0]) };
	// extract rotation matrix
	const vnl_matrix_fixed<float, 3, 3> sform_3x3 = sform_as_matrix.extract(3, 3, 0, 0);
	const vnl_matrix_fixed<float, 3, 3> qform_3x3 = qform_as_matrix.extract(3, 3, 0, 0);
	vnl_svd<float> sform_svd(sform_3x3.as_ref());
	vnl_svd<float> qform_svd(qform_3x3.as_ref());
	// extract offset
	const vnl_vector<float> sform_offset{ sform_as_matrix.get_column(3).extract(3, 0) };
	const vnl_vector<float> qform_offset{ qform_as_matrix.get_column(3).extract(3, 0) };
	// extract perspective
	const vnl_vector<float> sform_perspective{ sform_as_matrix.get_row(3).as_vector() };
	const vnl_vector<float> qform_perspective{ qform_as_matrix.get_row(3).as_vector() };
	// if sform_3x3 * inv(qform_3x3) is approximately and identity matrix then they are very similar.
	const vnl_matrix_fixed<float, 3, 3> candidate_identity{
	sform_svd.U() * vnl_matrix_inverse<float>(qform_svd.U()).as_matrix()
	};
	const bool spacing_similar{ sform_svd.W().diagonal().is_equal(qform_svd.W().diagonal(), 1.0e-4) };
	const bool rotation_matricies_are_similar{ candidate_identity.is_identity(1.0e-4) };
	const bool offsets_are_similar{ sform_offset.is_equal(qform_offset, 1.0e-4) };
	const bool perspectives_are_similar{ sform_perspective.is_equal(qform_perspective, 1.0e-4) };
	return rotation_matricies_are_similar && offsets_are_similar && perspectives_are_similar && spacing_similar;
	}();
	prefer_sform_over_qform = sform_and_qform_are_very_similar;
	}
	}
	else if (m_SFORM_Permissive) // sform is orthonormal
	{
	// Fix to deal with non-orthogonal matrixes
	// this approach uses SVD decomposition
	// maybe it is better to use nifti_make_orthog_mat44
	// to be consistent with other software?
	itkWarningMacro(<< this->GetFileName() << " has non-orthogonal sform");
	vnl_matrix_fixed<double, 3, 3> mat;
	for (int i = 0; i < 3; ++i)
	for (int j = 0; j < 3; ++j)
	{
	mat[i][j] = double{ this->m_NiftiImage->sto_xyz.m[i][j] };
	}
	vnl_svd<double> svd(mat.as_ref(), 1e-8);
	if (svd.singularities() == 0)
	{
	mat = svd.V() * svd.U().conjugate_transpose();
	mat44 _mat;
	for (int i = 0; i < 3; ++i)
	for (int j = 0; j < 3; ++j)
	{
	_mat.m[i][j] = static_cast<float>(mat[i][j]);
	}
	// preserve origin
	for (int i = 0; i < 3; ++i)
	{
	_mat.m[i][3] = this->m_NiftiImage->sto_xyz.m[i][3];
	}
	for (int i = 0; i < 4; ++i) // should be 0 0 0 1
	{
	_mat.m[3][i] = this->m_NiftiImage->sto_xyz.m[3][i];
	}
	this->m_SFORM_Corrected = true;
	return _mat;
	}
	}
	} // sform not NIFTI_XFORM_UNKNOWN

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)

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

Line 1976 in 3b04ad1

if (!prefer_sform_over_qform || this->m_NiftiImage->sform_code != NIFTI_XFORM_UNKNOWN)

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.