-
-
Notifications
You must be signed in to change notification settings - Fork 662
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
BUG: Fix matrix check rejecting NIFTI sform for small voxels #2701
BUG: Fix matrix check rejecting NIFTI sform for small voxels #2701
Conversation
Fixes InsightSoftwareConsortium#2674 NIFTI sform matrices were being called non-invertible if their determinant was < 10-5. Because the determinant scales with the product of the voxel spacings, any image with voxels smaller than about 20 microns would not have its direction cosines read from the sform. This caused errors when the header sform was valid, but the qform transform was not available as a fallback. Invertibility is now tested by checking the condition number of the singular value decomposition. If this is above the float epsilon, the matrix will be considered invertible. An alternative solution would be to consider a matrix invertible whenever the determinant is > 0, which is what itkMatrix does. The condition number test as implemented is a bit more conservative, but should allow spacings down to a nanometer or less.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
LGTM, but it would be good if someone else took a look too. Even better if someone else tried 😄
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@cookpa well done. Thanks for the contribution!
A simple test with some sample images here |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is a much more comprehensive and well-thought-out condition. Thank you for your contribution!
I have only loosely followed the issue/PR. Thanks for the effort all. @cookpa wondering whether the test you wrote in the separate repository should be put into ITK/you would be willing to do that. |
I'd be willing to write the test @jhlegarreta , is there an up to date guide on how to do that? I think I would need to upload the data files from that repo, because I had to manually edit the headers to remove the qform. |
💯 Thanks @cookpa. The ITK SW Guide contains guidelines (sections 9.4 C.23) to create ITK tests. Essentially:
each one reading the corresponding
HTH. |
Thanks @jhlegarreta - I have not forgotten this. I think this existing test code would work, if it was given files with the right headers: I'll look more into how to do that. |
Are there any plans to merge this into a 5.2.x bugfix? I'm running into this in with a high-resolution 0.3mm isotropic image I'm working with. |
|
||
const float condition = vnl_matrix_inverse<float>(mat.as_matrix()).well_condition(); | ||
// Check matrix is invertible by testing condition number of inverse | ||
if (!(condition > std::numeric_limits<float>::epsilon())) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why ! > instead of <= ?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
! >
and <=
give different results for floating point NaNs. With the former, if the condition is small or is a NaN then ...
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Right! Pretty subtle though. Something like if (!isfinite(condition) || (condition < epsilon))
might make the intention clearer.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@seanm you're right, that is a nicer way to phrase it.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The ! >
approach dates back forever and its purpose is easily recognizable by us oldtimers. If we are nonetheless worried that it is opaque, that could be addressed in a comment.
If we instead go with the newer approach, instead of isfinite
we might want to be checking isnan
-- unless we are sure that a positive infinity should return false
in our current context, or unless we are sure that positive infinity is impossible in our current context.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The ! > approach dates back forever and its purpose is easily recognizable by us oldtimers
Maybe I'm an even older timer than you then, because my foggy old memory forgot about this trick. :)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The well_condition() function returns the ratio of the smallest to largest singular values. I think it was returning nan when both were 0. If there is an infinity something has gone seriously wrong, and the enclosing IsAffine function should also return false in that case.
The commit message is great, but alas my maths are too rusty... @vfonov maybe you'd care to take a look? |
Looks like there's still some corner cases here, #3333 |
Fixes #2674
NIFTI sform matrices were being called non-invertible if their determinant was <
10-5. Because the determinant scales with the product of the voxel spacings, any
image with voxels smaller than about 20 microns would not have its direction
cosines read from the sform. This caused errors when the header sform was valid,
but the qform transform was not available as a fallback.
Invertibility is now tested by checking the condition number of the singular
value decomposition. If this is above the float epsilon, the matrix will be
considered invertible.
An alternative solution would be to consider a matrix invertible whenever the
determinant is > 0, which is what itkMatrix does. The condition number test as
implemented is a bit more conservative, but should allow spacings down to a
nanometer or less.
PR Checklist
Refer to the ITK Software Guide for
further development details if necessary.