[Pkg-exppsy-pynifti] [pymvpa] [pynifti] NiftiDatatset bug?

[Scott]

Michael Hanke wrote:
>> Either that, or the data should be handled privately (in NiftiImage) in 
>> an intelligent way that can translate the correct dimensions into the 
>> actual dimensions in the data array.  Perhaps a function like 
>> getFullData(t=True, u=False, v=False, w=False) which extracts a 4d 
>> (T,X,Y,Z) even if there's only one sample?  Otherwise you get bugs like 
>> this.  For instance, I don't know ahead of time if i'll be loading a 3D, 
>> 4D, or 4D with 1 timepoint volume (which apparently is treated 
>> differently, as my cludgy hack works that way) and every time i need to 
>> use the array I have to check for where the dimensions are.
>>     
>
> Why? The dimensions are always the reverse of x,y,z,t,u,v,w, no?
>   
yes, that's the problem - i have to check whether t is dim 0 or if I 
need to stick the data inside [] before concatenating, since singleton 
w, v, u, t are not preserved in the array.
> I wonder what your '4d with 1' image header looks like. See this:
>   In [10]: nim.header['dim']
>   Out[10]: [4, 4, 3, 2, 1, 1, 1, 1]
>
>
> This is exactly what I would expect -- and if I got you right this is
> what you want too, right?
>
>   
Yes I believe this is what happens.  The problem is that, by design, dim 
3 nii's are shaped as (x,y,z) whereas dim 4's are shaped as (t, x, y, z) 
meaning I cannot directly concatenate dim3's and dim4's together along 
the time/samples dimension.  it is not a terrible flaw, just one that 
requires care and attention to avoid error.
>> To further question this - what would happen to the Nifti data if the 
>> UVW dims are not 1?  Presumably this is rare... but would they keep 
>> adding on to the beginning of the array?
>>     
> Yes.
>
>   
>>  So the number of samples now 
>> falls in the 2nd, 3rd, or 5th position? 
>>     
> No, it is always the first.
>
>   
>> Now perhaps I extract a subset 
>> of U,V,W 4-D volumes and squeeze away the preceding singleton dims - now 
>> there's no guarantee to the number of dims I'm left with.  Or are U,V,W 
>> mapped after Z? And what about single-slice, or other 1d and 2d images 
>> encoded in an nii?  It seems like wanting to squeeze out dims stands at 
>> ends with dim compatibility in this case.
>>     
> Hmm, I did not get this.
>
>   
okay, seems like we're not on the same page.  i'm not sure why you think 
samples are in the first position if u,v,w are added to the beginning of 
the array - perhaps you're thinking of the data when handled by pymvpa?  
i mean just the nifti data itself.  if u,v,w are mapped to the 
beginning, then it must push time/samples into a subsequent position.  
Then, to illustrate my point, let's say we load an array which has 2 u 
samples, 2 timepoints, but only one slice of data in X (for argument's 
sake).  So, the dims might be:

[2, 1, 1, 2, 1, 64, 64]

At which point, if I extracted just the 'volume' at time=1, u=1, and 
squeezed out the singleton dimensions, I'd be left with

[64, 64]

Now there are two problems - perhaps my code expects v or w dims, even 
knowing i've extracted u and t.  Or, perhaps I expected a 3D volume (ie 
some library function doesn't know to expect single-slice data).  
alternatively, you could have an 'intelligent' squeeze in nifitimage 
that always preserves the last three dimensions...  But simply put, if 
the dims keep shifting, every time we use niftimage.data, we have to 
check the shape and infer which dimension is which.
>> While I'm proposing crazy ideas, it would also be simpler if the index 
>> just happened to be reversed - (Z,Y,X,T,W,V,U).  That way the proper 
>> volume is always consistent and extra dimensions can be safely dropped.  
>> I think in Matlab, nii's are typically loaded like this (I usually load 
>> using Freesurfer tools, not sure about SPM).  It would also be easy to 
>> concatenate samples or uvw coords without having to permute the array first.
>>     
>
> I don't see the advantage -- apart from the fact that this would change
> the order of information in the file and increase IO operation time.
>
>
>   
guess i don't have enough experience with array granularity to groc this 
fully.  i thought that on disk, the order of data changes fastest in the 
spatial dims, then t, then in the uvw's.  so, how does this map into 
numpy - numpy is changing fastest at the last dims?  would it be a large 
performance hit to simply permute the array after reading?  if this is a 
serious issue, then perhaps the best solution is just a method to cast 
the data into the requested space, like NiftiImage.asNd(self, N=4, t=0, 
u=0, v=0, w=0), and perhaps NiftiImage.extend(self, other, dim='t') 
which calls asNd

the advantage is just that, IMHO, you get to preserve dim consistency 
and wanting the important XYZ's up front, which frankly is the only way 
I would have thought to do this :-\
>> I understand these ideas are probably over the top complicated in the 
>> face of backwards compatibility - so perhaps they can be my 'modest 
>> proposal'.  But I hope they illustrate my point - IMHO it's crucial to 
>> easily extract and concatenate 3D volumes from a NiftiImage without 
>> having to figure out what dimension 'samples' is stored in  every time 
>> it's necessary.
>>     
>
> We should not confuse the needs of PyNIfTI with the ones of PyMVPA.
>
> But in general it is quite simple: The last three dimensions/axes are
> always the spatial ones, the one in front of them is the time and the
> rest is something else.
>
>   
>> In other words - if you are heading to a new release with pynifti, yes, 
>> i think this issue is worth addressing there :)
>>     
>
> Alright, but I still have no clear concept of what the 'issue' is in the
> first place -- let discuss some more ;-)
>
>   
to be concise - i want to always know which dim says what without having 
to check.   So, let's say i want to get a beta estimate from a nifti 
file, precomputed from some glm package.  fsl writes each beta to a 
single image, but freesurfer writes them all together, concatenated 
along time.  So, in pseudocode-ish:

def getBeta(im, range=0)
nim = niftiimage(im, load=True)
return nim.data[range]

would be a simple way to handle this.  but it is flawed, since we would 
need to check first if samples were in the first dim or not.  every time 
we want to deal with a niftiimage, we need something like

if nim.niftiheader['dim'][0]==3:#pray we don't have fewer than 3
    return N.asarray([nim.data])#to force 4d so we never do this again
return nim.data[range]

or if we expect u,v, w:

tdim = len(nim.data.shape)-4
if tdim < 0:
    return N.asarray([nim.data])
else return something fancier by subscripting range into the dimension 
tdim, using code i don't feel like looking up now :)
 
Or, let's say we just want to collect a bunch of images from disk and 
save them into a single array, concatenated by time (or any other dim).  
We would have to check each at runtime to see whether it needs upcasting 
to 4d before concatenation.

Here's a related question:  Does pynifti handle the 
niftiheader['dim'][0] spec when it writes out?  for instance, lets say 
you loaded a timecourse, did some math on it which reduces to a single 
stat, say the mean, or a beta weight/svm weight, etc, then wrote this to 
disk.  is the header set to dim[0]=3 or dim[0]=4, if it's kept the 
original header to preserve orientation, etc?  I feel like not every 
program handles this in the same way, which is perhaps one reason i've 
seen this as a problem (i've jumped through a lot of fmri software...).  
I think some of my nii's on disk actually have dim[0]=4 even when 
there's only one timepoint, and some set dim[0]=3.

Cheers
-S