-
-
Notifications
You must be signed in to change notification settings - Fork 665
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 failing isAffine check due to insuffcient precision in calculations #3339
BUG: Fix failing isAffine check due to insuffcient precision in calculations #3339
Conversation
This fixed my issues with my high resolution human brains and rat brains I've been working with. |
9f23b19
to
abc6b33
Compare
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.
@gdevenyi Thanks for addressing my commits, Gabriel! 👍
You might still consider adjusting the subject line of the commit text to get it at at most 72 chars, for example:
BUG: Fix failure isAffine due to insufficient precision in calculations
Because currently GitHub breaks the last word of the subject line in two: "caluc……lations" at
abc6b33
But it's not a big deal to me, of course! Approved anyway!
Depending upon where I look in ITK, I see different usage for converting from auto xd0 = static_cast<double>(xf * xf);
auto xd1 = double{ xf * xf };
double xd2 = xf * xf; Which is preferred? |
I prefer static_cast( ). It is verbose but very explicit about the desired behavior. |
Good question 🤔 I would prefer not to use |
Does the answer change for potentially lossy casts such as between |
Sure! |
I answered giving my preference under the condition of converting from |
And by its omission, am I to assume that you would not use the third form in either the lossy or lossless scenarios? |
Personally, I would not. I know that it is a lossless conversion, but novice developers benefit from having the conversion made painfully obvious. |
In this particular case, conversion from ITK/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx Lines 1756 to 1760 in abc6b33
In this case, implicit conversion |
Still, novice developers need to learn about |
While I have already approved this PR, I think it would be worth adding a minimal unit test, that checks the fix. Something like this? From @vfonov at #3333 (comment)
Would that be the most simplified (minimal) test case? Update: I see now, the ITK/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx Line 1754 in e00bad0
|
How about the following two; which do you prefer when auto xd1 = double{ xf * xf };
double xd3{ xf * xf }; |
abc6b33
to
ffd4204
Compare
I made a simple test (also showing the limits of the original code from http://callumhay.blogspot.com/2010/10/decomposing-affine-transforms.html:
testing results:
Outputs:
|
@vfonov Thanks for the test code! So I see, case2 (below here) is considered non-affine by both
Just curious. Obviously switching from |
Very minor preference for the second option. |
I don't disagree that C++ programmers need to learn and be skilled. I would just like to point out that there is a large ITK community that are primarily physicians (mathemeticians/physicists/ect), secondarily medical imaging researchers, and have a tertiary skill set of being C++ programmers. When possible, assisting with their skillset development has benefits to the ITK community. |
Just recently, I did a few contributions that should make conversions in ITK easier to read: by a new conversion function, |
d5dc985
to
e6baf23
Compare
Hopefully sorted, I had the absolute test data path instead of the relative one from the NIFTI IO directory. |
using ReaderType = itk::ImageFileReader<ImageType>; | ||
ReaderType::Pointer reader = ReaderType::New(); | ||
reader->SetFileName(argv[1]); | ||
reader->Update(); | ||
ImageType::Pointer image = reader->GetOutput(); |
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.
Please consider using itk::ReadImage
instead. Just one line of code should so 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.
Maybe some extra checks here to see that the read image is as expected...? Or at least, check that it is a valid image?
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.
- I do not have any C++ object-oriented black magic to anything beyond copying the ITK official example, https://examples.itk.org/src/io/imagebase/readanimage/documentation
- This test is to explicitly see if the NiftiIO isAffine check fails on the random image I had in my collection which was causing issues
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.
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 test is to explicitly see if the NiftiIO isAffine check fails on the random image I had in my collection which was causing issues
I guess you could still add some checks that the read image is valid, for example:
if (image == nullptr)
{
std::cerr << "Error: The read image is null!\n";
return EXIT_FAILURE;
}
if (image->GetBufferPointer() == nullptr)
{
std::cerr << "Error: the buffer of the read image is null!\n";
return EXIT_FAILURE;
}
if (image->GetBufferedRegion().GetNumberOfPixels() == 0)
{
std::cerr << "Error: The read image has no pixels!\n";
return EXIT_FAILURE;
}
Right?
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.
@gdevenyi Did you still consider adding the extra checks that I suggested here? In order to ensure that the reading does actually yield a valid image?
if (image == nullptr)
{
std::cerr << "Error: The read image is null!\n";
return EXIT_FAILURE;
}
if (image->GetBufferPointer() == nullptr)
{
std::cerr << "Error: the buffer of the read image is null!\n";
return EXIT_FAILURE;
}
if (image->GetBufferedRegion().GetNumberOfPixels() == 0)
{
std::cerr << "Error: The read image has no pixels!\n";
return EXIT_FAILURE;
}
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 can be condiserably abbreviated via itk::ReadImage
and ITK_TRY_EXPECT_NO_EXCEPTION
. We can do that in a follow-up PR. Let's merge this, as it has dragged on for quite a long time.
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.
Both done in #3448.
itk_add_test(NAME itkNiftiSmallVoxelsAffinePrecisionTest | ||
COMMAND ITKIONIFTITestDriver itkNiftiImageIOTest13 | ||
DATA{Input/SmallVoxels_AffinePrecision.nii.gz}) |
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.
What about using itkNiftiImageIOTest2 instead, like this:
itk_add_test(NAME itkNiftiSmallVoxelsAffinePrecisionTest
COMMAND ITKIONIFTITestDriver itkNiftiImageIOTest2
${ITK_TEST_OUTPUT_DIR} itkNiftiIOShouldSucceed true DATA{Input/SmallVoxels_AffinePrecision.nii.gz})
Would that also work?
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.
That test with the original code does not fail.
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.
That test with the original code does not fail.
Why not? Originally you added your SmallVoxels_AffinePrecision.nii.gz to an existing itk_add_test
. My suggestion (above here) would be to do a new itk_add_test
, referring to the old itkNiftiImageIOTest2 cxx file and your SmallVoxels_AffinePrecision.nii.gz file. But that would not work either?
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.
I don't get why itkNiftiImageIOTest2 doesn't fail, because the test file SmallVoxels_AffinePrecision.nii.gz
definitely causes an exception for me (running ANTs PrintHeader built against recent ITK).
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.
Reverting the change, and running both tests.
$ ctest -R AffinePrecision --output-on-failure
Test project /home/gdevenyi/projects/itk-build
Start 1157: itkNiftiSmallVoxelsAffinePrecisionTest
1/2 Test #1157: itkNiftiSmallVoxelsAffinePrecisionTest ............................***Failed 0.04 sec
ExceptionObject caught !
itk::ExceptionObject (0x563a086a0fa0)
Location: "unknown"
File: /home/gdevenyi/projects/ITK/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx
Line: 1988
Description: ITK ERROR: ITK only supports orthonormal direction cosines. No orthonormal definition found!
Start 1158: itkNiftiSmallVoxelsAffinePrecisionTest_fromitkNiftiImageIOTest2
2/2 Test #1158: itkNiftiSmallVoxelsAffinePrecisionTest_fromitkNiftiImageIOTest2 ... Passed 0.03 sec
50% tests passed, 1 tests failed out of 2
Label Time Summary:
ITKIONIFTI = 0.07 sec*proc (2 tests)
Total Test time (real) = 0.22 sec
The following tests FAILED:
1157 - itkNiftiSmallVoxelsAffinePrecisionTest (Failed)
Errors while running CTest
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.
So itkNiftiImageIOTest2 just cannot be used for SmallVoxels_AffinePrecision.nii.gz, for some unknown reason, right?
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.
So itkNiftiImageIOTest2 just cannot be used for SmallVoxels_AffinePrecision.nii.gz, for some unknown reason, right?
Correct.
The difference between the two sets of code appears to be
itk::NiftiImageIO::Pointer io = itk::NiftiImageIO::New();
auto imageReader = ImageReaderType::New();
imageReader->SetImageIO(io);
Where the NiftiimageIO ImageIO is explicitly set in itkNiftiImageIOTest2
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.
@gdevenyi Indeed, itkNiftiImageIOTest2 does explicitly set ImageIO, but it does also do a "read" without explicitly setting ImageIO. By calling itk::IOTestHelper::ReadImage
:
input = itk::IOTestHelper::ReadImage<ImageType>(std::string(arg2)); |
Anyway, it seems to me that the CMake code of itkNiftiIOShouldSucceed was already incorrect. It says:
itk_add_test(NAME itkNiftiIOShouldSucceed
COMMAND ITKIONIFTITestDriver itkNiftiImageIOTest2
${ITK_TEST_OUTPUT_DIR} itkNiftiIOShouldSucceed true DATA{${ITK_DATA_ROOT}/Input/LittleEndian.hdr,LittleEndian.mhd,LittleEndian.img})
ITK/Modules/IO/NIFTI/test/CMakeLists.txt
Lines 55 to 57 in 6fc6164
itk_add_test(NAME itkNiftiIOShouldSucceed | |
COMMAND ITKIONIFTITestDriver itkNiftiImageIOTest2 | |
${ITK_TEST_OUTPUT_DIR} itkNiftiIOShouldSucceed true DATA{${ITK_DATA_ROOT}/Input/LittleEndian.hdr,LittleEndian.mhd,LittleEndian.img}) |
The arguments of itk_add_test don't seem to match correctly with arg1, arg2, and prefix at
ITK/Modules/IO/NIFTI/test/itkNiftiImageIOTest2.cxx
Lines 32 to 42 in 8c2c654
if (argc != 4) | |
return EXIT_FAILURE; | |
char * arg1 = argv[1]; | |
char * arg2 = argv[2]; | |
char * prefix = argv[3]; | |
int test_success = 0; | |
using ImageType = itk::Image<short, 3>; | |
using ImagePointer = ImageType::Pointer; | |
if ((strcmp(arg1, "true") == 0) && WriteNiftiTestFiles(prefix) == -1) |
When I run it locally, arg1
== "itkNiftiIOShouldSucceed", arg2
== "true", and prefix
== "C:/X/bin/I/ITK/ExternalData/Testing/Data/Input/LittleEndian.hdr". While arg1
should have been the boolean ("true" or "false"), and arg2
should have been the input file name. Right?
I think itkNiftiSmallVoxelsAffinePrecisionTest should be more like the following, if you would still like to use the existing itkNiftiImageIOTest2:
itk_add_test(NAME itkNiftiSmallVoxelsAffinePrecisionTest
COMMAND ITKIONIFTITestDriver itkNiftiImageIOTest2
${ITK_TEST_OUTPUT_DIR} true DATA{${ITK_DATA_ROOT}/Input/SmallVoxels_AffinePrecision.nii.gz} SomeSpecificPrefix)
I don't know what exactly the prefix ("SomeSpecificPrefix") should be. But it should be the last argument to itk_add_test
.
My two cents, Niels
e6baf23
to
5735185
Compare
@gdevenyi Would a smaller test image file be possible? I just noticed that the file you added for this test, SmallVoxels_AffinePrecision.nii.gz is more than 1 MB! https://data.kitware.com/#item/627539e24acac99f42c8929a If I understand correctly, the only part of the image that is relevant to this test is the 4x4 matrix at |
I took the problematic naturalistic image I had, and zeroed all the voxels to make it maximally compressible. I don't know the NIFTI format well enough to successfully modify the file any further without accidentally breaking the important properties. |
Does this work? I just copied the header to a smaller array using RNifti |
// First make sure the bottom row meets the condition that it is (0, 0, 0, 1) | ||
{ | ||
float bottom_row_error = itk::Math::abs(mat[3][3] - 1.0f); | ||
double bottom_row_error = itk::Math::abs(mat[3][3] - 1.0); |
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.
is float epsilon several lines below (line 1769) intentional?
if (bottom_row_error > std::numeric_limits<float>::epsilon())
I thought that it might have been forgotten.
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.
is float epsilon several lines below (line 1769) intentional?
@issakomi Good question. But I wonder, isn't std::numeric_limits<T>::epsilon()
rather arbitrary anyway? I just don't know, but if mat[3][i]
must be close to zero (for i
= 0, 1, 2), then those matrix elements might as well be compared against std::numeric_limits<T>::min()
, which if much smaller than std::numeric_limits<T>::epsilon()
.
Honestly I'm just wondering what is intended as well. I have no suggestion 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.
I think the bottom row should be precisely (0,0,0,1) - but I guess it's possible there might be some epsilon error when the inverse is calculated
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.
I think it was just forgotten.
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.
there might be some epsilon error when the inverse is calculated
OK, but std::numeric_limits<T>::epsilon()
isn't just some epsilon error. It is specifically the difference between 1 and the least value greater than 1 that is representable by floating point type T
. https://www.cplusplus.com/reference/cfloat/ That's why I wonder if this specific epsilon should also be used for values that should be close to zero (rather than close to one). But I just don't know, honestly.
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.
One further idea here, currently the test is made with the sum as each entry added. Should we instead test the individual entries?
Changing to ::min()
, or testing the individual entries will be substantially stricter than the existing implementation as well.
@cookpa Thanks Philip. Your reduced nii.gz does trigger the "ITK ERROR: ITK only supports orthonormal direction cosines" when I try to read it by |
Fixes InsightSoftwareConsortium#3333 Co-authored-by: Vladimir S. FONOV <vladimir.fonov@gmail.com>
5735185
to
b8285fb
Compare
@gdevenyi That's really great, Gabriel, thanks! Can you possibly also remove the large (1.042 MB) version of SmallVoxels_AffinePrecision.nii.gz from https://data.kitware.com/#item/627539e24acac99f42c8929a ? |
As far as I understood, I did. The PR only includes the HASH of the small replacement file. |
But it's still there at https://data.kitware.com/#item/627539e24acac99f42c8929a ! Right? Isn't it possible to remove the file from https://data.kitware.com/#item/627539e24acac99f42c8929a ? |
data.kitware.com is a bit silly. Following your link there's no option for deletion, but if I navigate to my files via my personal account I can delete it. Its now gone. It would've never been fetched anyways. |
@hjmjohnson I believe we're just waiting on you here. Can you please re-review? |
Following up again here, this is ready to go. |
using ReaderType = itk::ImageFileReader<ImageType>; | ||
ReaderType::Pointer reader = ReaderType::New(); | ||
reader->SetFileName(argv[1]); | ||
reader->Update(); | ||
ImageType::Pointer image = reader->GetOutput(); |
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 can be condiserably abbreviated via itk::ReadImage
and ITK_TRY_EXPECT_NO_EXCEPTION
. We can do that in a follow-up PR. Let's merge this, as it has dragged on for quite a long time.
Fixes #3333
Co-authored-by: Vladimir S. FONOV vladimir.fonov@gmail.com
PR Checklist
No API changes were made (or the changes have been approved)
No major design changes were made (or the changes have been approved)
Updated API documentation (or API not changed)
Added license to new files (if any)
Added test
NO Added Python wrapping to new files (if any) as described in ITK Software Guide Section 9.5
NO Added ITK examples for all new major features (if any)