Skip to content
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

[ENH] ITK IO to be LPS-aware #257

Merged
merged 29 commits into from
Oct 30, 2018
Merged

Conversation

ashgillman
Copy link
Collaborator

@ashgillman ashgillman commented Oct 4, 2018

Solution for #242
Depends #248

Still one TODO - reading orientation from DICOM. But works for now with NIfTI by assuming HFS (actually it calls it unknown, this is just treated like HFS in #248).

@ashgillman
Copy link
Collaborator Author

@rijobro perhaps you could add this PR to the Image Geometry project for bookkeeping?

@KrisThielemans
Copy link
Collaborator

@ashgillman could you rebase this on master resolving the conflict that appears after I squash-merged #248?

@ashgillman
Copy link
Collaborator Author

Done

@ashgillman
Copy link
Collaborator Author

There are some methods here that return image pointers. Should these be changes to be shared pointers?

@ashgillman ashgillman self-assigned this Oct 4, 2018
Copy link
Collaborator

@KrisThielemans KrisThielemans left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

overall looks good, but I haven't checked the detail.

Regarding convert_ITK_to_STIR and read_file_itk returning bare pointers, it's not nice, but I don't mind too much as they are internal static functions. If we both to change them, they should return unique_ptr (might be somewhat trickier than you think as we need to preserve backwards compatibility with non-C++11 compilers)

for ( it.GoToBegin(); !it.IsAtEnd(); ++it, ++stir_iter) {
itk::Point<double,3U> itk_coord;

itk_coord[0] = -double(it.Get()[0]);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why the minus signs here?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@rijobro I assume this can go now, I'll remove it. Best you double check on your data.

switch (exam_info_sptr->patient_position.get_position()) {
case PatientPosition::unknown_position:
// If unknown, assume HFS
// TODO: warning?
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes please.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, best when we set patient_position to unknown I think

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

sure. would be then when reading. I guess here.

I have no idea if ITK reads patient position for other file formats (non-DICOM). would be great if they did...

const BasicCoordinate<3,int> &min_indices,
bool is_displacement_field)
ITK_coordinates_to_STIR(const itk::ImageBase<3>::PointType &itk_coord,
const STIRImageType stir_image,
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

reference to STIRImageType?

Copy link
Collaborator Author

@ashgillman ashgillman Oct 4, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

const STIRImageType& stir_image? Sure, but curious why?

Does that prevent copying the whole image data?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes. (potentially) big objects should always be passed by (preferably const) reference. Note for stuff like int, float etc. Really should add this to our developers-guide...

@ashgillman
Copy link
Collaborator Author

ashgillman commented Oct 4, 2018

Since rebasing I am getting a numerical exception. Has anyone seen this before? Div by 0?

31: Test command: /home/gil2a4/dev/UCL@github.com/STIR/build/src/test/test_IO_DiscretisedDensity "/home/gil2a4/dev/UCL@github.com/STIR/src/test/input/test_ITKNiftiOutputFileFormat.in"
31: Test timeout computed to be: 10000000                                                                                                                  
31:                                                                                                                                                        
31: WARNING: FactoryRegistry:: overwriting previous value of key in registry.
31:      key: None
31: Testing OutputFileFormat parsing function...
31: WARNING: will overwite files called STIRtmp*
31: 
31: About to write the image to disk...
31: OK!
31: 
31: About to read the image back from disk...
5/5 Test #31: test_IO_DiscretisedDensity_test_ITKNiftiOutputFileFormat.in ..........***Exception: Numerical  0.23 sec  

Update: nevermind, having more luck with debugger

@ashgillman
Copy link
Collaborator Author

Issue was that itk_image->Clone() returned an empty image with none of the same data nor metadata as itk_image. I opted to allow conversion to modify itk_image, since it is not used elsewhere. Alternatively, this seems to be the way to duplicate an image - rather verbose in classic ITK style.

https://itk.org/Wiki/ITK/Examples/SimpleOperations/ImageDuplicator

@KrisThielemans
Copy link
Collaborator

As far as I'm concerned, this can be merged. ok?

Regarding DICOM, I've posted a question to the ITK Community list

@KrisThielemans
Copy link
Collaborator

The ITK list is now replaced by a forum. My question is here. The answer points to the section of the manual on DICOM IO. The relevant bit is at the end.

@ashgillman, want to do that before merging, or later?

@ashgillman
Copy link
Collaborator Author

Before probably makes sense, I'll have a look into it.

@ashgillman
Copy link
Collaborator Author

ashgillman commented Oct 24, 2018

Hmm, no current code is not correct. Writing is correct it seems, but if I read a NIfTI then write out again it is reversed in the SI and AP directions

The below images are the phantom acquisitions at UCL mMR acquired with @bathomas. The first is the vendor reconstruction converted to .nii with dcm2niix. This is to remove any implication that the error is in the DICOM-specific reading code, but all results are identical if I ask STIR to work directly with the DICOM. The second image is a .nii produced after reading the .nii and writing back to .nii (testing self-consistency).

Coronal view:
image
image

@ashgillman
Copy link
Collaborator Author

ashgillman commented Oct 24, 2018

Where's the error?

If I do a STIR reconstruction on raw data, and save the output to .nii, the reconstruction appear correctly oriented w.r.t. vendor reconstructions, so the write seems to be working and the error I guess must be in the reader.

Is the reorienter working?

Here I print the Orientation and Origin before and after the reorienter.

1 0 0
0 -1 0
0 0 1
[-359.216, 356.894, -283.957]
1 0 0
0 1 0
0 0 -1
[-359.216, -358.693, 128.387]

Side note
Why this LAS input axis order? The input NIfTI has affine
[[ -2.08626008, -0. , 0. , 359.21588135],
[ 0. , 2.08626008, 0. , -356.89422607],
[ 0. , 0. , 2.03125 , -283.95703125],
[ 0. , 0. , 0. , 1. ]]
NIfTI affine is defines w.r.t. RAS, so image is therefore in LAS

This is the expected behaviour. Notice though that the axes flipping is the SI and AP directions, the same as I am getting an issue with.

Does manually messing with the reorienter help?

When I manually change the reorienter for HFS from RAS to RPI, the resulting image out has shifted its origin (an expected consequence), but the axes directions are still incorrect.

The reorienter works on the meta data, what about voxel data?

So maybe the reorienter isn't actually re-ordering the voxel data. Or maybe the iterator is smart in some way and iterates not in stored order?

So I added this snippet to check the iterator data is actually changing:

  std::cerr << itk_image->GetDirection();
  std::cerr << itk_image->GetOrigin() << std::endl;
  typename ITKImageType::Pointer new_itk_image
    = orient_ITK_image<ITKImageType>(exam_info_sptr, itk_image);
  std::cerr << new_itk_image->GetDirection();
  std::cerr << new_itk_image->GetOrigin() << std::endl;

  typedef itk::ImageRegionConstIterator<ITKImageType> IteratorType;
  IteratorType it (itk_image, itk_image->GetLargestPossibleRegion());
  IteratorType nit (new_itk_image, itk_image->GetLargestPossibleRegion());
  int n_different = 0;
  for (it.GoToBegin(); !it.IsAtEnd(); ++it, ++nit)
  {
    n_different += (it.Get() != nit.Get());
  }
  std::cerr << "n different: " << n_different << std::endl;

which returns

1 0 0                                                                                           
0 -1 0                                                                       
0 0 1                               
[-359.216, 356.894, -283.957]                                                                               
1 0 0                                           
0 1 0                                     
0 0 -1                                                      
[-359.216, -358.693, 128.387]         
n different: 3465846

So the reorienter is changing the voxel data...

Test the test

check that if I craft the reorienter to do nothing that no voxels change (reorient to LAS, the same as this particular input):

1 0 0
0 -1 0
0 0 1
[-359.216, 356.894, -283.957]
1 0 0
0 -1 0
0 0 1
[-359.216, 356.894, -283.957]
n different: 0

So...

I must admit I am stumped...

@KrisThielemans
Copy link
Collaborator

@rijobro, sounds good to me. People might get fed-up seeing the warning, but we have so many repetitive warnings in STIR already that one more isn't going to hurt :-;

@rijobro
Copy link
Collaborator

rijobro commented Oct 26, 2018

@ashgillman, turns out I can make changes directly in your branch. So I did that to add the warning. I'll do it via PR in the future, but this one seemed like a tiny change.

@ashgillman
Copy link
Collaborator Author

ashgillman commented Oct 29, 2018

Alright, I was also able to reduce some of the explicit templating, however it seems many are still required when the output type is assigned to a pointer, as the output type can't be determined.

Edit: Sorry, missed some more comments made, pushing changes now.

Copy link
Collaborator

@KrisThielemans KrisThielemans left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

2 super-tiny suggestions.

Awaiting testing by @rijobro.

@@ -161,7 +162,8 @@ <h3>Changed functionality</h3>
</li>
<li>
The orientation of images read/written via ITK has changed. It should now be correct if the patient
was in HFS position. This is currently not checked (and for many file formats, never can be checked).
was in HFS position for NIfTI (this can't be verified due to lack of metadata in the format), or if
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I find this a bit awkward to read. I suggest

was in HFS position. Unfortunately, for most file formats (such as NifTI) this can't be verified due to lack of metadata in the format. If the files are in DICOM format, as long as the file contains the correct metadata, other patient positions are supported (currently HFS, FFS, HFP or FFP).

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes was a bit muddled in my head as I tried to get the words out. Much better

src/IO/ITKImageInputFileFormat.cxx Outdated Show resolved Hide resolved
KrisThielemans and others added 2 commits October 29, 2018 17:59
Co-Authored-By: ashgillman <ashley.gillman@my.jcu.edu.au>
@ashgillman
Copy link
Collaborator Author

Sorry, should have added the skip ci tag.

@rijobro
Copy link
Collaborator

rijobro commented Oct 29, 2018

Tested by performing non-rigid transformation of nifti image with nifti displacement field. Result was the same as before. Any other tests needed?

@rijobro
Copy link
Collaborator

rijobro commented Oct 29, 2018

Did a small PR. By default, you were setting is_displacement = false for both image types. We want it to be true for STIRImageMulti.

@KrisThielemans
Copy link
Collaborator

to me, setting the defaults doesn't make a lot of sense. do we need them? Anyway, whatever works really. I'm eager to press the "squash and merge" button!

@rijobro
Copy link
Collaborator

rijobro commented Oct 29, 2018

I meant to say that the default wasn't being overridden, which means that the displacement image was being treated as if it was a deformation.

Looks good to me now! (once Ash accepts my PR)

[FIX] default to displacement not deformation
@ashgillman
Copy link
Collaborator Author

ashgillman commented Oct 29, 2018

Ah I see, I hadn't gone high enough with the is_displacement_fields. Thanks.

I think the default makes sense where @rijobro left it (I have previously removed the others), since you don't want to have to choose it for single images.

I also kept the default in ITK_coordinates_to_STIR_physical_coordinates() and renamed is_displacement_field in this funtion to is_relative_coordinate, to make it explicit the actual reason. So when converting ITK_pixel_to_STIR_pixel() for multi images, we need to know if it is_displacement_field, and if so, we convert ITK_coordinates_to_STIR_physical_coordinates() such that the vector is_relative_coordinate. Again, in ITK_coordinates_to_STIR_physical_coordinates() you don't want to have to specify whether displacement if dealing with normal coordinates.

But I'm not sure you really care! Point is I also am happy now :)

@KrisThielemans
Copy link
Collaborator

done?

@ashgillman
Copy link
Collaborator Author

Yep

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants