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

[Michael Hanke]

On Mon, Jan 12, 2009 at 03:52:31PM -0500, Scott wrote:
> 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.
Indeed, _leading_ singletons would not be preserved -- actually, they were
never there in the first place.


>> 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)  
No. dim 3 images are shaped (z,y,x) !!

> 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.
You can perfectly do that.

>>> 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?   
No, I mean PyNIfTI, actually there is not much difference between both.

> 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.   
Hmm, not sure what you mean.

> 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]
No. They would be (1, 1, 2, 2, 64, 64, 1) !!
>
> 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]
No, you'd get (64, 64, 1).

> 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.
I guess we have found the main problem. There is no squeezing in
pynifti. The order of dimensions is simply reversed wrt the dim vector
in the header, but nothing else is done to them.

Can you perhaps provide one of the images that cause problem? Maybe I am
still not getting it...


>>> 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
Shoudl not be necessary if my above speculations are correct.


> 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.
You always know: the last is always X, the one before is always Y, ....

>   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]
That should work perfectly.


> 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
No.

> 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]
Don't get this. Above returns you a 3d volume, but this would return a
4d image -- is that intended?

> 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?
Sure.

>  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.
Let the code speak: ;-)

Python 2.5.2 (r252:60911, Nov 14 2008, 19:46:32) 
[GCC 4.3.2] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> from nifti import *
>>> nim=NiftiImage('example4d')
>>> nim.header['dim']
[4, 128, 96, 24, 2, 1, 1, 1]
>>> m=nim.data.mean(axis=0)
>>> m.shape
(24, 96, 128)
>>> nim2=NiftiImage(m, nim.header)
>>> nim2.header['dim']
[3, 128, 96, 24, 1, 1, 1, 1]
>>> 


That is how it should go (IMHO).


Michael


-- 
GPG key:  1024D/3144BE0F Michael Hanke
http://apsy.gse.uni-magdeburg.de/hanke
ICQ: 48230050