WIP: fix and refactor par / rec read / write #227

matthew-brett · 2014-02-16T00:06:52Z

par / rec read write wasn't tested and I'd broken it some time ago by
some refactoring of file IO and of the ArrayProxy class.

I fixed my breakages and also:

Added PAR / REC data to test against - thanks Michael
added some methods to PARRECHeader to make it work with the new
ArrayProxy class
Changed TR scaling from milliseconds to seconds; I think that's right.
Added tests of the script and the module

I've tested against the nifti file that Philips generated from the same
data -- from : http://psydata.ovgu.de/philips_achieva_testfiles/

A remaining issue is that our affine differs from the affine Philips
writes into its own conversion to Nifti::

In [26]: np.allclose(our_conv.get_sform(), our_conv.get_qform())
Out[26]: True

In [27]: np.allclose(philips_conv.get_sform(), np.eye(4))
Out[27]: True

In [28]: our_conv.get_qform()
Out[28]: 
array([[  -3.6499,   -0.    ,    1.8356,  123.6628],
       [  -0.    ,   -3.75  ,   -0.    ,  115.617 ],
       [   0.8605,   -0.    ,    7.7866,  -58.2506],
       [   0.    ,    0.    ,    0.    ,    1.    ]])

In [29]: philips_conv.get_qform()
Out[29]: 
array([[  -3.65  ,   -0.0016,    1.8356,  125.4881],
       [   0.0016,   -3.75  ,   -0.0004,  117.4916],
       [   0.8604,    0.0002,    7.7866,  -28.3411],
       [   0.    ,    0.    ,    0.    ,    1.    ]])

mih · 2014-02-17T19:04:35Z

Looks good. Where does the shift in the xform come from. Is that relative to the FoV origin, or scanner center?

matthew-brett · 2014-02-17T19:20:01Z

Actually - they seem to generate the same affine ('--origin=fov', '--origin=scanner')...

matthew-brett · 2014-03-04T19:36:30Z

Michael - any comments on the affine thing? Would you mind checking that?

This does not address the issue mentioned in nipygh-227 that Philips own NIfTI conversion yields a different offset column in the affine.

mih · 2014-03-15T15:23:03Z

I took a first look. For me --origin 'fov' yields a different affine than 'scanner'. I added a test for that.

I did not figure out why the offset in the Philips conversion differs. I'll try some other converters to get a second opinion before I start surgery.

Part of the problem in nipygh-227

mih · 2014-03-15T16:17:38Z

So I think I fixed the biggest issue (pushed to my fork). Here is the default ("scanner") affine that is produced now:

[[  -3.64994708    0.            1.83564171  123.66276611]
 [   0.           -3.75          0.          115.617     ]
 [   0.86045705    0.            7.78655376  -27.91161211]
 [   0.            0.            0.            1.        ]]

This is pretty close to the Philips interpretation:

array([[  -3.65  ,   -0.0016,    1.8356,  125.4881],
       [   0.0016,   -3.75  ,   -0.0004,  117.4916],
       [   0.8604,    0.0002,    7.7866,  -28.3411],
       [   0.    ,    0.    ,    0.    ,    1.    ]])

It seems like for the first two offset values they shift by a half-voxel (we have 3.75x3.75 in-plane resolution) in this file. But I have no clue where the missing 0.43 for the last axis should come from.

Any ideas?

mih · 2014-03-15T16:20:51Z

Hmm, interesting: the acquisition console screenshot for the data says slice thickness 6mm and slice gap 1mm. But exported PAR/REC, XML/REC and NII all claim 2mm gap or 8mm slice thickness. not sure who is at fault here...

screenshot here: http://psydata.ovgu.de/philips_achieva_testfiles/conversion/console_screenshot.png

parrec uses iteration over lines of a text file; implement this in the Openers wrapper and test.

Thanks Michael H for releasing it.

Test against values from converted NIfTI header

Adapt to BinOpener fileobj object as input. Change returned TR value to seconds instead of milliseconds. Add get_data_offset and get_slope_inter methods for compatibility with standard ArrayProxy object. Add option to set default scaling as input to header __init__

Use new PAR/REC data and Philips converted nifti image to test against.

Change 'transversal' to 'transverse'; 'sagital' to 'sagittal'. I think these are the standard names: http://en.wikipedia.org/wiki/Sagittal_plane

numpy strings expected to be bytes, not strings. Files must get opened with binary flag.

Make bytes / strings explicit. gzip doesn't appear to allow opening for read in text mode in Python 3.

This does not address the issue mentioned in nipygh-227 that Philips own NIfTI conversion yields a different offset column in the affine.

Part of the problem in nipygh-227

Newlines between methods; docstring formatting and adding docstring return value descriptions; expanded description of format.

read_img_data trying to set offset of 0 even when offset is already set to 0, therefore breaking for headers that can't do ``set_data_offset``. Check whether value is already 0, don't set if so.

par header format enforces data offset of 0. Implement set_data_offset for par header enforcing this.

Replace with recommended ``dataobj.get_unscaled`` method.

matthew-brett · 2014-07-11T01:44:09Z

The PAR file obviously does think there's an 6mm slice thickness with a 2mm
gap, from the image information:

thick   gap 

6.000  2.000
6.000  2.000
6.000  2.000
6.000  2.000

There are 9 slices, and the FOV is:

.    FOV (ap,fh,rl) [mm]                :   240.000  70.000  240.000

That fits, no? (6 * 9 + 2 * 8)?

By the way, I notice that get_voxel_size doesn't include the slice gap, I
guess that was intentional. I added a comment in any case.

I played around with the affine but didn't get anywhere. I refactored a bit,
for my own understanding, but it didn't reveal anything much.

matthew-brett · 2014-07-11T01:48:12Z

Sorry - I have terrible internet here in Cuba - what does the XML file say as to the slice gap?

coveralls · 2014-07-11T01:50:29Z

Coverage increased (+0.73%) when pulling 58f7968 on matthew-brett:fixing-parrec into 9fe2e9a on nipy:master.

Test we're getting the same affine as we were (it may be the wrong affine, but at least it is our affine). Test we're getting the same-ish 3x3 component as the Philips generated NIfTI image.

A few typos and thinkos

* Rephrase the transformations in terms of LPS instead of RAS in the hope of making the code easier to read * Remove safety-net for the rotations, I'm confident the r2agui approach and euler2mat do the same thing. For clarity, unpack the order of rotations with euler2mat in case we need to fix the order later. * Use voxel FOV to get FOV offset in voxel space, then transform to get the FOV offset in mm It's possible this refactoring breaks or fixes handling of non-transverse images - we really need some of those to test against.

larsoner · 2014-08-04T18:50:28Z

Hopefully, yeah. I might be able to open a PR into your PR. Here are the cases I think we need to detect:

Mis-ordered acquisition (e.g., x, y, t, z).
Missing data (e.g., run too short, doesn't contain full dataset).
Multiple data types in same recording with different scale factors (phase contrast vs complex).
DTI b-values and vectors are recorded.

One thing to correct:

Orientation in saggital images seems incorrect.

@mrjeffs is that a brief summary of the issues we've hit thus far?

larsoner · 2014-08-04T18:57:06Z

Sorry, I meant "open a PR into your branch" not "into your PR"...

matthew-brett · 2014-08-04T19:01:19Z

:) - mi casa tu casa.

We should certainly fix orientation of sagittal images - I hope Michael's tests will help out with that.

mih · 2014-08-05T12:24:08Z

Here are the simple test scans:

http://psydata.ovgu.de/philips_achieva_testfiles/

the new ones are in the conversion2/ folder. The README in the top-folder should have the necessary bits of information.

larsoner · 2014-08-06T00:54:27Z

@matthew-brett are you planning to look at those sample files and check out the saggital orientation stuff sometime this week? In the meantime I'm working with @mrjeffs to get representative sample files from our system that satisfy the problematic conditions I outlined above. I'm going to try to look into it in the next couple days to make some progress by the end of the week.

matthew-brett · 2014-08-06T01:00:05Z

I am planning to get to this tomorrow - yes - so if you can get something done in the next few days, that would be a great help.

yarikoptic · 2014-08-06T18:41:01Z

Totally forgot about Michael's collection of phantom data and yesterday collected some on a local Philips box, aiming to come up with quite obscure setup (tilt, resolution different in all dimension, small matrix size, slice gaps) covering major acquisition protocols (T1/T2/T2*/DTI/fieldmap) and also exporting in DICOMs (as a single file per image), PAR/REC, NIFTI (non FSL one). Overall size of the .git for such a collection is just 1.6MB, while I think it could cover quite a wide range of test cases -- I would appreciate if you have a look: https://github.com/yarikoptic/nitest-balls1 , may be they would come of some use? ;)

mih · 2014-08-06T19:31:47Z

Just, quickly: I had a brief look at the affines produced by the current state of the branch for the scans with cor,sag,tra slice orientation and no further rotation or translation. All I know so far is that only axial/transverse orientiation assigns proper axis labels in FSLView. Each image should have the marker (more or less) clearly visible with min/max settings of 1000/20000. One marker for anterior, two for inferior and three for left.

That is all I got...

matthew-brett · 2014-08-13T03:24:15Z

Sorry to be slow on this one, I will review these in the next couple of days and get back to y'all.

larsoner · 2014-08-13T03:26:47Z

Cool. I assume you haven't had time to work on the transformation issues, then...? :)

matthew-brett · 2014-08-13T03:27:28Z

Sorry - not yet.

larsoner · 2014-08-13T03:33:11Z

It almost seems like you have more responsibilities than working on nibabel or something...

Function for chaining dot product across arrays, particularly useful for dot product of affine components

Michael's new data had enough variation to test most combinations of transformations. The key seems to be to transform first to the Philips axes (anterior-posterior, inferior-superior, right-left), then apply the rotations, before finally taking into RAS output space. Add the PAR files for Michael's test data set and test against the current versions of the affines, which I checked by reslicing the data with the affines (see check_parrec_reslice.py script).

matthew-brett · 2014-08-21T18:46:01Z

Michael's new dataset was very useful in sorting this out.

The affine transforms for Michael's conversion2 dataset appear to be about
right now. I think they are OK for conversion1 and Yarick's dataset too, at
least, they are the same the Philips niftis, which look right.

The FSL nifti affines for conversion2 seemed to be a bit of a mess, so I was
checking the conversions by eye ...

We're still not dealing well with sequences with interleaved slices and so on.

Probably the way to fix that would be to make a new ImageArrayProxy object
that will unscramble the slices / echo / dynamic orders on read.

But I believe the orientation is right now, although I can't test the effect
of different rotation orders yet.

I think this current PR is ready to merge, do y'all agree? Then rebase for
the slice / volume sorting?

Now I know the format a bit better, I have a feeling we could get the sorting
etc done in time for the release - let's have a try over the next week or so? I'll have very little internet, but I'll do my best to keep up.

larsoner · 2014-08-21T18:53:49Z

Eventually an enhanced ImageArrayProxy object would be the way to go. For now, my PR into your branch should make it so that XYTZ and XYZT orders both work. But for other more complex/inconsistent orderings (if they exist) it will still fail.

I don't see a commit here -- did the rotations work already with your current changeset, or is the commit pending? I could check to see whether our data are oriented properly now.

matthew-brett · 2014-08-21T18:57:42Z

Sorry - just sending commits now...

larsoner · 2014-08-21T20:10:44Z

Shoot, it created merge conflicts with my PR. We have a few options:

You rebase your branch, then I rebase mine, and we review/merge my PR before or after release, whichever makes more sense.
I rebase my PR into your branch, then we review/merge my PR, then you rebase before release.

Let me know which you prefer...

larsoner · 2014-08-21T20:12:00Z

Oh, by the volume/slice sorting I guess you mean my PR... yeah I'm fine with you rebasing this, merging it, and the I'll rebase and open a PR into upstream/master. This would make the most sense.

matthew-brett · 2014-08-21T20:22:01Z

OK - sounds reasonable. I'll have a look at your PR today or
tomorrow, I hope I'll be able to get back to you soon - but internet
here in Cuba is very slow, and I need to go for a 20 minute drive to
get to a hotel to get internet at all.

Would you mind having a look at the latest commits for a review? With
a little review I will merge, and that will make rebasing more
straightforward...

larsoner · 2014-08-21T20:39:40Z

Edits look reasonable to me (although I am certainly no transformation expert), but perhaps more importantly the changes have fixed the orientations in our data! Awesome.

matthew-brett · 2014-08-21T20:48:45Z

OK - thanks - then I will merge - I think this is an improvement. Can
you then rebase on master? Will do merge now.

matthew-brett · 2014-08-21T20:54:14Z

Done.

larsoner · 2014-08-21T20:59:37Z

Done: #253.

This does not address the issue mentioned in nipygh-227 that Philips own NIfTI conversion yields a different offset column in the affine.

Part of the problem in nipygh-227

mih added a commit to mih/nibabel that referenced this pull request Mar 15, 2014

TST: Verify PARREC affines differ in the right places

a2ba263

This does not address the issue mentioned in nipygh-227 that Philips own NIfTI conversion yields a different offset column in the affine.

mih added a commit to mih/nibabel that referenced this pull request Mar 15, 2014

BF: parrec2nii ignored FH image offset

16aa374

Part of the problem in nipygh-227

matthew-brett mentioned this pull request Apr 10, 2014

PAR/REC failure to load ('BinOpener' object is not iterable) #237

Closed

matthew-brett and others added 14 commits July 5, 2014 16:59

NF: add iter method to Openers utility class

674fee0

parrec uses iteration over lines of a text file; implement this in the Openers wrapper and test.

NF: add PAR / REC data for testing

2e5080f

Thanks Michael H for releasing it.

NF: add skeleton for tests for PAR/REC conversion

9ccdb8c

Test against values from converted NIfTI header

TST: add test of parrec2nii script

3d8d364

Use new PAR/REC data and Philips converted nifti image to test against.

RF: change slice orientation names for par/rec

07cf98b

Change 'transversal' to 'transverse'; 'sagital' to 'sagittal'. I think these are the standard names: http://en.wikipedia.org/wiki/Sagittal_plane

RF: some python 3 fixes for par / rec stuff

245a36f

numpy strings expected to be bytes, not strings. Files must get opened with binary flag.

RF: python 3 fixes for openers test

6db3b51

Make bytes / strings explicit. gzip doesn't appear to allow opening for read in text mode in Python 3.

TST: Verify PARREC affines differ in the right places

da2cdfc

This does not address the issue mentioned in nipygh-227 that Philips own NIfTI conversion yields a different offset column in the affine.

BF: parrec2nii ignored FH image offset

c92524f

Part of the problem in nipygh-227

DOC: pep8 + expand / format docstrings for parrec

f3aa7ea

Newlines between methods; docstring formatting and adding docstring return value descriptions; expanded description of format.

RF: avoid setting data offset of 0 if already set

7c9022a

read_img_data trying to set offset of 0 even when offset is already set to 0, therefore breaking for headers that can't do ``set_data_offset``. Check whether value is already 0, don't set if so.

RF: error if trying to set par data offset !=0

a2bacc2

par header format enforces data offset of 0. Implement set_data_offset for par header enforcing this.

RF: remove use of deprecated read_img_data

58f7968

Replace with recommended ``dataobj.get_unscaled`` method.

matthew-brett added 4 commits July 11, 2014 11:26

TST: add regression tests for affine

984d173

Test we're getting the same affine as we were (it may be the wrong affine, but at least it is our affine). Test we're getting the same-ish 3x3 component as the Philips generated NIfTI image.

DOC: fix errors in top-level docstring

2965054

A few typos and thinkos

DOC: fix tiny typo in docstring

49e33b0

larsoner mentioned this pull request Aug 10, 2014

ENH: Fix PAR/REC support matthew-brett/nibabel#4

Closed

7 tasks

matthew-brett added 3 commits August 15, 2014 18:37

RF: remove unused inport in funcs

42b5105

NF: dot_reduce utility function

d7e9c3e

Function for chaining dot product across arrays, particularly useful for dot product of affine components

matthew-brett merged commit bb0973f into nipy:master Aug 21, 2014

matthew-brett deleted the fixing-parrec branch October 25, 2014 04:31

grlee77 pushed a commit to grlee77/nibabel that referenced this pull request Mar 15, 2016

TST: Verify PARREC affines differ in the right places

f581725

This does not address the issue mentioned in nipygh-227 that Philips own NIfTI conversion yields a different offset column in the affine.

grlee77 pushed a commit to grlee77/nibabel that referenced this pull request Mar 15, 2016

BF: parrec2nii ignored FH image offset

f2b7f4c

Part of the problem in nipygh-227

WIP: fix and refactor par / rec read / write #227

WIP: fix and refactor par / rec read / write #227

Uh oh!

Conversation

matthew-brett commented Feb 16, 2014

Uh oh!

mih commented Feb 17, 2014

Uh oh!

matthew-brett commented Feb 17, 2014

Uh oh!

matthew-brett commented Mar 4, 2014

Uh oh!

mih commented Mar 15, 2014

Uh oh!

mih commented Mar 15, 2014

Uh oh!

mih commented Mar 15, 2014

Uh oh!

matthew-brett commented Jul 11, 2014

Uh oh!

matthew-brett commented Jul 11, 2014

Uh oh!

coveralls commented Jul 11, 2014

Uh oh!

larsoner commented Aug 4, 2014

Uh oh!

larsoner commented Aug 4, 2014

Uh oh!

matthew-brett commented Aug 4, 2014

Uh oh!

mih commented Aug 5, 2014

Uh oh!

larsoner commented Aug 6, 2014

Uh oh!

matthew-brett commented Aug 6, 2014

Uh oh!

yarikoptic commented Aug 6, 2014

Uh oh!

mih commented Aug 6, 2014

Uh oh!

matthew-brett commented Aug 13, 2014

Uh oh!

larsoner commented Aug 13, 2014

Uh oh!

matthew-brett commented Aug 13, 2014

Uh oh!

larsoner commented Aug 13, 2014

Uh oh!

matthew-brett commented Aug 21, 2014

Uh oh!

larsoner commented Aug 21, 2014

Uh oh!

matthew-brett commented Aug 21, 2014

Uh oh!

larsoner commented Aug 21, 2014

Uh oh!

larsoner commented Aug 21, 2014

Uh oh!

matthew-brett commented Aug 21, 2014

Uh oh!

larsoner commented Aug 21, 2014

Uh oh!

matthew-brett commented Aug 21, 2014

Uh oh!

matthew-brett commented Aug 21, 2014

Uh oh!

larsoner commented Aug 21, 2014

Uh oh!

Uh oh!