nipy/nibabel — pull request #278 by grlee77 (2014-11-18 01:36:39+00:00)
Title: ENH:  PAR/REC data ordering and b-vector normalization

Additional PAR/REC functionality was added (building upon the recent merge of #253, #261 and #264 as implemented by @matthew-brett and @Eric89GXL)

Aside from a couple of minor fixes, there are two new features added:

**1.) ability to normalize b-values & b-vectors**

The --norm-bvs argument to parrec2nii enables the auto-scaling, otherwise the old behavior is maintained. If all b-vectors were properly normalized to 1.0 and the b-values were stored correctly, this new functionality is not be needed, so it is disabled by default.

I implemented this because I have some datasets where the .PAR header came from a custom pulse sequence where the b-value label was left at 3000 for all volumes (even the b=0 ones!), and the only way to determine the b=0 volumes was that the gradient vector is [0, 0, 0] for these. Scaling the b-values by the norm of the gradient vector gives the proper b-value array. (I did not rescale the b-values/vectors for arrays where the norm was within [1.0 +/- 1e-2] to avoid causing the b-value array to have non-unique but nearly identical b-values such as [2999.99, 3000.01, etc] due to numerical precision of the stored b-vector arrays.)

The non-normalized behavior is still required for other cases, e.g. when using a product DWI sequence with 15 different b-values and requesting only the trace volumes to be stored. In this case, the gradient vectors will all be [0, 0, 0], but the b-values are an array of 15 unique values. If the norm of the gradient vectors had been used, this would have caused the b-value array to return all zeros instead of the appropriate values.

**2.) store the labels corresponding to the fourth data dimension**
I think something along these lines is extremely important and would be interested on feedback on how best to implement it. For now, I search through a set of potentially varying image_defs keys such as "label type", "echo number", etc. For any that have non-unique values, I maintain a dictionary with the corresponding indices. This dictionary can then be used later by the user to provide custom resorting of the data or reshaping to >4D, etc. I did not implement any custom resorting options within the classes, but for now just write out the dictionary of orderings out to a .json file in parrec2nii and leave it up to the user to resort in the desired manner.

An example of where this is useful can be illustrated for a simple ASL experiment with 8 dynamics and 2 label types (dataobj.shape[3] = 16). Without labeling along the fourth dimension, we will not know which volumes correspond to `label type = 1` and which correspond to `label type = 2` (i.e. control vs. tag volumes). One logical data ordering for the 16 total volumes is:

{
"dynamic scan number": [
1,
1,
2,
2,
3,
3,
4,
4,
5,
5,
6,
6,
7,
7,
8,
8
],
"label type": [
1,
2,
1,
2,
1,
2,
1,
2,
1,
2,
1,
2,
1,
2,
1,
2
]
}

However depending on how the data was exported, the ordering within the .PAR could equally well have been:

{
"dynamic scan number": [
1,
2,
3,
4,
5,
6,
7,
8,
1,
2,
3,
4,
5,
6,
7,
8
],
"label type": [
1,
1,
1,
1,
1,
1,
1,
1,
2,
2,
2,
2,
2,
2,
2,
2
]
}

Without saving some sort of ordering info along with the .nii, there is no clean way to tell these two cases apart and properly process the resulting .nii files. I doubt json is the optimal way to do this, but it is easily readable and illustrates the desired functionality. I am certainly open to other suggestions. The .json output corresponds to the new `PARRECArrayProxy._dim_4_labels` property.

The `sorted_labels` function of the `PARRECHeader` also has the ability to include per-slice info such as slice angulation, and store this into `PARRECArrayProxy._dim_3_4_labels`. This allows .PAR/.REC files corresponding to scans such as 3-plane localizers, etc. to have the raw array data be read in along with its corresponding dimension info. These volumes with varying slice angulation still cannot be read into a PARRECImage object, though, because `PARRECArrayProxy.get_slice_orientation` only allows unique values for slice angulation so that a unique affine can be determined for the `PARRECImage`.