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

[Scott]

Michael Hanke wrote:
>
> No. dim 3 images are shaped (z,y,x) !!
>
>   
sorry, i knew i was gonna mix up the x,y,z's - for all below, of course 
you are right about the dim names, but it's the dim sizes that I posit 
are the problem, so let's assume those are correct.
>> 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.
>   
if and only if i manually configure them all to be the right data shape, 
every time i want to handle the array.  ie lets say i want to get all my 
beta estimates from a series of FSL files, i should be able to just:

d = numpy.asarray([NiftiImage(f, load=True).data for f in file_list])

But if one of them had more than one time point (or just had dim[0]=4) 
this would fail, and if all of them had equal # of timepoints, the array 
would be shape=(#files, #times, z,y,x) instead of 4d.  If they had 
unequal #tp this would throw an exception.  Or, i could do something like:

d = None
for f in file_list:
    nim = NiftiImage(f)
    imd = nim.data.reshape(nim.timepoints, nim.volextents*)
    if d is None:
       d = imd
    else:
       d.extend(imd)

    or,
    imd = nim.data
    if len(imd.shape) < 4:
       imd = [imd]
    d.extend(imd)

which is obviously manageable, but more work than should be necessary 
(especially for nifti noobs).  apologies for being liberal with the 
list<->numpy.array syntax, i wish they were more consistent...

i should also note that just keeping an array of [niftiimages] simply 
postpones this problem, as any code will have to deal with the dimension 
variability when it views the data.  For instance, if I want a function 
to operate on every sample, the elegant way would be:

for vol in nim.data:
    some code...

but I instead must:

for s in range(nim.timepoints):
    if nim.data.ndim==4: # ignoring u,v,w for simplicity here, but this 
should be managed too
       vol = nim.data[s]
    else:
       vol = nim.data   #Indexing data[s] here would extract a z slice!

    some code with vol...

hence my proposal of a method like NiftiImage.asNd(self, N=4, t=0, u=0, 
v=0, w=0) which can map the data array as necessary, specifically, 
either to a 3d volume from any dim size, or to any dim size as required 
from a single 3d volume.  The implementation should be simple, just an 
extrapolation of what I've put above, but robust for all possible dim 
sizes (perhaps just reshape into 7d, index at the kwargs provided if 
necessary, and reshape to the last N dims).  Perhaps even an 
implementation of extend would be useful too so that NiftiImage objects 
can be directly concatenated
>
> 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.
>
>   
well, yes, reversed but not including the leading singletons, meaning 
that every dim's position in the shape of the array is dependent on how 
many leading singletons were ignored.  Even though you've never called 
squeeze, the result is the same as if they had been squeezed out of the 
reference 7d array.  That is, the order is not what's important here (i 
know it's constant), it's which dimension is living in shape[0], whether 
it's a (z,y,x), time, or another.
> You always know: the last is always X, the one before is always Y, ....
>
>   
I've never needed to know what's in shape[-1], i always need to know 
what's in shape[0].  the point is not that i can't figure it out, it's 
that I *must* figure it out every time I want to manipulate the array.  
IMHO it's too easy to make a mistake doing this (or not doing it!), and 
end up with mismatched arrays errors or mistaking t and z, like the 
niftidataset bug which set slices as the number of samples.
>>   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.
>
>
>   
no, this only works if the first dimension is time, but not if it's a 3d 
array - because then it would return the range in z instead of t!

Also, to address your confusion below, I was unclear when I wrote 
nim.data[range] above.  My intention was that range could be a number, 
or a list, so that we can extract one or more samples from nim, and they 
would all be returned together in 4d.  For clarity, I meant something 
more like:

def getBetas(im, timepoints=None):
    nim = niftiimage(im, load=True)   # Let's assume 4d, even if only 1 
timepoint
    if timepoints is None:
       return nim.data
    if not operater.isSequenceType(timepoints):
       timepoints=[timepoints]   # If index is always an array, we never 
loose a dim in indexing, and will preserve 4d so we can always expect a 
(samples x Z x X x Y) array returned
    return nim.data[timepoints]
>> 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?
>
>   
yes, because we need to concatenate samples together along the 4th (or 
rather, 1st of four) dimension, so everything needs to have that dim.  
(alternatively, we can force everything to be 3-d, but then we must 
first split any array with more than 1 sample in it before recombining 
them all!)

>>  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
>
>   
Yes - but this works because m itself was reduced to 3 dimensions, when 
mean squeezed out the singleton on axis=0.  perhaps a function doesn't 
reduce the dims, like if we're iterating over sets and we don't know how 
many there are, something like [fun(vol) for vol in nim.data].  The 
behavior would be more like:

 >>> r = numpy.random.rand
 >>> r = numpy.random.rand(2, 3, 4)
 >>> r.shape
(2, 3, 4)

 >>> nim = NiftiImage(r) # As above
 >>> nim.data.shape
(2, 3, 4)
 >>> nim.header['dim']
[3, 4, 3, 2, 1, 1, 1, 1]

 >>> nim2 = NiftiImage(numpy.asarray([r]))    # 4d with 1 timepoint
 >>> nim2.header['dim']
[4, 4, 3, 2, 1, 1, 1, 1]
 >>> nim2.data.shape
(1, 2, 3, 4)

 >>> nim3 = NiftiImage(numpy.asarray([r]), header=nim.header) # Fed 3d 
header
 >>> nim3.data.shape
(1, 2, 3, 4)
 >>> nim3.header['dim']   # Oops! The dims are still inconsistent with 
nim and r
[4, 4, 3, 2, 1, 1, 1, 1]

Now we have a 4d image with 1 timepoint, even if we fed it the correct 
header info!  I know this is by design, not a bug - but I think it's 
important to be clear that this will cause other bugs if a programmer is 
not careful.  If you're still unclear what the problem is, try writing 
nim and nim3 to disk, load them, and then try to work with them together 
without examining the header or shape of the array, or your code 
assuming which one is in which shape.

Even if pynifti handles all of this perfectly (should you decide there 
is a perfect way to handle this...), there is still no guarantee that 
other software will handle singleton dims in the same manner - hence why 
I think pynifti should have a proactive way to deal with data in the 
shape required by any given function.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.alioth.debian.org/pipermail/pkg-exppsy-pynifti/attachments/20090112/d7edc62a/attachment-0001.htm