Skip to content

Parrec sort fix #409

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

Merged
merged 14 commits into from
Mar 15, 2016
Merged

Parrec sort fix #409

merged 14 commits into from
Mar 15, 2016

Conversation

grlee77
Copy link
Contributor

@grlee77 grlee77 commented Feb 13, 2016

This PR is an attempt to address #408.

I added a new boolean option upon loading a .PAR/.REC file called strict_sort. If it is set to True the get_sorted_slice_indices method of PARRECHeader uses several additional keys in the call to np.lexsort. strict_sort=False corresponds to the previously existing sort mode. This option is made available via a --strict_sort command line option to parrec2nii as well.

In general, I think it is probably preferable to use strict_sort=True in all cases, but I have left the default as False for now to ensure backwards compatibility. Perhaps after a deprecation/testing period we could make True the new default. In particular, it may not be true that for any given DTI scan the b-vals/vectors would be sorted in exactly the same order as previously (If the user is using the b-vals/vectors as exported by parrec2nii this should not be an issue, but could cause problems if mixed with a previously stored table from the other approach).

There could be some debate as to the preferred order of priority for the sort keys. I chose one that seemed logical to me. As before, slices is the lowest priority so that dimension 3 of the data is still always slices and is_full is the highest so that partially acquired scans get moved to the end for possible discarding. The sorting is a bit complicated, but I tried to add comments to make the reasoning clearer. The full set of potential sort keys implemented is:

slice_nos
echos
cardiac phases
b-vectors
b-values
asl_labels
dynamics
vol_numbers()
is_full

@grlee77
Copy link
Contributor Author

grlee77 commented Feb 13, 2016

If I get a chance, I may scan a phantom to generate a couple of additional complex .PAR files to add to the tests.

@matthew-brett
Copy link
Member

That would be very good.

If you happen to get a chance to scan a phantom with (different) rotations about all three axes, that would always help a lot - we don't know what order PAR / REC uses for these at the moment.

@@ -631,10 +632,14 @@ def __init__(self, info, image_defs, permit_truncated=False):
permit_truncated : bool, optional
If True, a warning is emitted instead of an error when a truncated
recording is detected.
strict_sort : bool, optional, keyword-only
If True, a larger number of header fields are used while sorting
the REC data array.
Copy link
Member

Choose a reason for hiding this comment

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

Perhaps mention that this can give a different order of slices.

@matthew-brett
Copy link
Member

Can you also add a test for strict_sort in the script, in nibabe/test_scripts.py?

I'm afraid I still haven't got my head round the sorting, but I'll play with the code soon to see if I can understand it.

@grlee77
Copy link
Contributor Author

grlee77 commented Feb 25, 2016

@matthew-brett
As requested above (but unrelated to this PR), .PAR/.REC files corresponding to various oblique image orientations are made freely available here:
https://github.com/grlee77/parrec_oblique

I have not yet had the time to fully compare the affines produced by the scanner to those generated by parrec2nii.

@grlee77
Copy link
Contributor Author

grlee77 commented Feb 25, 2016

I tried to look at the details for the coveralls failure above, but I get the following error:

Source Not Available

The file "/home/travis/build/nipy/nibabel/venv/lib/python2.7/site-packages/nibabel/parrec.py" isn't available on github. Either it's been removed, or the repo root directory needs to be updated.

any ideas why that is?

@effigies
Copy link
Member

@grlee77 See #331.

@matthew-brett
Copy link
Member

Great - thanks for the data - I make a PR testing against it at #420

is_full = vol_is_full(slice_nos[initial_sort_order],
self.general_info['max_slices'])

# have to "unsort" is_full to match the other sort keys
Copy link
Member

Choose a reason for hiding this comment

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

This logic took me a bit to work out - is there a good test for it?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Do you mean the is_full case specifically? If so, I don't think I tested that aside from verifying that it still passes the existing tests. None of the new data sets I added have missing data, although I suppose we could truncate one of the files if desired.

I don't have an exhaustive set of tests for all available sort keys, but I did collect a few .PAR files with 4-6 dimensions to use as tests. These are the new .PAR files I added in tests/data:

  • T1_dual_echo.PAR: A 4D data set: (x, y, z, echo)
  • T1_3cho_mag_real_imag_phase.PAR: A 5D data set: (x, y, z, echo, image_type)
  • ASL_3D_Multiecho.PAR: A 6D data set: (x, y, z, echo, 'label type', dynamics)
    The tests make sure the orders are sorted as expected for these cases. I did not collect any data with cardiac phases, for example.

One other test I can think of would be to manually scramble the order of the entries in one of the .PAR files and make sure the sort keys that come out at the end still match the unscrambled case.

@grlee77
Copy link
Contributor Author

grlee77 commented Mar 15, 2016

@matthew-brett
I think this is now ready for review again. The primary changes since last time you reviewed are:

  • I removed the previously untested _direction_numbers routine since it was based on a 'diffusion' field in image_defs that isn't even present in the V4 .PAR files that would have been the only case to call this routine (see FIX: properly handle V4 .PAR files with diffusion info #426)
  • In the tests for strict sort I made sure they pass even when the image_defs array is randomly sorted via np.shuffle.
  • Since V4 .PAR files have no gradient orientation number key, I added vol_numbers back in as a key in the strict sort routine so the diffusion directions would get assigned unique volume numbers even for V4 .PAR files. As was already the case for vol_is_full this has to be done after the initial stage of sorting. In this way the sorting is robust even to random shuffling of the image info lines from the .PAR file.

bvecs = bvecs[0]
# rotate bvecs to match stored image orientation
permute_to_psl = ACQ_TO_PSL[self.get_slice_orientation()]
bvecs = apply_affine(np.linalg.inv(permute_to_psl), bvecs)
return bvals, bvecs

def get_def(self, name):
""" Return a single image definition field. """
Copy link
Member

Choose a reason for hiding this comment

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

Add "or None if field not present" ?

@matthew-brett
Copy link
Member

Great - thanks for the changes - nearly there.

It's good to avoid backslashes for line continuation if possible - https://www.python.org/dev/peps/pep-0008/#id17

@grlee77
Copy link
Contributor Author

grlee77 commented Mar 15, 2016

I made the suggested changes above.

I also simplified the final lines of _strict_sort_order. I was able to remove the "unsorting" entirely. Now it is an initial sort as before, followed by a second sort that is just:

initial_sort_order[np.lexsort((vol_nos, is_full))]

The docstring was expanded to (hopefully!) better describe what is going on in the function`.

@nibotmi
Copy link
Contributor

nibotmi commented Mar 15, 2016

☔ The latest upstream changes (presumably #426) made this pull request unmergeable. Please resolve the merge conflicts.

@grlee77
Copy link
Contributor Author

grlee77 commented Mar 15, 2016

rebased to fix merge conflicts with #426

@matthew-brett
Copy link
Member

Thanks for your patience - it looks good - merging.

matthew-brett added a commit that referenced this pull request Mar 15, 2016
MRG: Add option for PARREC sorting with more header fields

This PR is an attempt to address #408.

I added a new boolean option upon loading a .PAR/.REC file called `strict_sort`.  If it is set to `True` the `get_sorted_slice_indices` method of `PARRECHeader` uses several additional keys in the call to `np.lexsort`.  `strict_sort=False` corresponds to the previously existing sort mode.  This option is made available via a `--strict_sort` command line option to `parrec2nii` as well.

In general, I think it is probably preferable to use `strict_sort=True` in all cases, but I have left the default as `False` for now to ensure backwards compatibility.  Perhaps after a deprecation/testing period we could make `True` the new default.  In particular, it may not be true that for any given DTI scan the b-vals/vectors would be sorted in exactly the same order as previously (If the user is using the b-vals/vectors as exported by `parrec2nii` this should not be an issue, but could cause problems if mixed with a previously stored table from the other approach).

There could be some debate as to the preferred order of priority for the sort keys.  I chose one that seemed logical to me.  As before, slices is the lowest priority so that dimension 3 of the data is still always slices and `is_full` is the highest so that partially acquired scans get moved to the end for possible discarding.  The sorting is a bit complicated, but I tried to add comments to make the reasoning clearer.  The full set of potential sort keys implemented is:
```
slice_nos
echos
cardiac phases
b-vectors
b-values
asl_labels
dynamics
vol_numbers()
is_full
```
@matthew-brett matthew-brett merged commit 23bc068 into nipy:master Mar 15, 2016
@grlee77
Copy link
Contributor Author

grlee77 commented Mar 15, 2016

Great! This will be very useful for some of the data we acquire.

matthew-brett added a commit that referenced this pull request Mar 18, 2016
MRG: add a function to retrieve .PAR data labels for the 4th dimension

This was pulled out of the previously submitted #278 to make reviewing easier.  This complements the functionality in the recently merged #409.

The primary enhancement is a new method in the PARREC header called `get_volume_labels` that lists which .PAR fields the data stored along the 4th dimension correspond to.  (i.e. all data >4D is still collapsed into 4 dimensions exactly as before, but labels can be retrieved to make the data ordering unambiguous).

The script interface to parrec can write these labels out as a CSV file with one row per volume.
matthew-brett added a commit that referenced this pull request May 4, 2016
MRG: expand API docs for PAR/REC

Add API documentation for the PAR/REC enhancements introduced by #409 and #427
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.

4 participants