<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
  <meta content="text/html;charset=ISO-8859-1" http-equiv="Content-Type">
</head>
<body bgcolor="#ffffff" text="#000000">
<br>
<br>
Michael Hanke wrote:
<blockquote cite="mid:20090112211117.GA16553@runkel" type="cite"><br>
  <pre wrap="">No. dim 3 images are shaped (z,y,x) !!

  </pre>
</blockquote>
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.<br>
<blockquote cite="mid:20090112211117.GA16553@runkel" type="cite">
  <pre wrap=""></pre>
  <blockquote type="cite">
    <pre wrap="">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.
    </pre>
  </blockquote>
  <pre wrap=""><!---->You can perfectly do that.
  </pre>
</blockquote>
if and only if i manually configure them all to be the right data
shape, every time i want to handle the array.&nbsp; ie lets say i want to
get all my beta estimates from a series of FSL files, i should be able
to just:<br>
<br>
d = numpy.asarray([NiftiImage(f, load=True).data for f in file_list])<br>
<br>
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.&nbsp; If they
had unequal #tp this would throw an exception.&nbsp; Or, i could do
something like:<br>
<br>
d = None<br>
for f in file_list:<br>
&nbsp;&nbsp;&nbsp; nim = NiftiImage(f)<br>
&nbsp;&nbsp;&nbsp; imd = nim.data.reshape(nim.timepoints, nim.volextents*)<br>
&nbsp;&nbsp;&nbsp; if d is None:<br>
&nbsp;&nbsp;&nbsp; &nbsp;&nbsp; d = imd<br>
&nbsp;&nbsp;&nbsp; else:<br>
&nbsp;&nbsp;&nbsp; &nbsp;&nbsp; d.extend(imd)<br>
<br>
&nbsp;&nbsp;&nbsp; or,<br>
&nbsp;&nbsp;&nbsp; imd = nim.data<br>
&nbsp;&nbsp;&nbsp; if len(imd.shape) &lt; 4:<br>
&nbsp;&nbsp;&nbsp; &nbsp;&nbsp; imd = [imd]<br>
&nbsp;&nbsp;&nbsp; d.extend(imd)<br>
<br>
which is obviously manageable, but more work than should be necessary
(especially for nifti noobs).&nbsp; apologies for being liberal with the
list&lt;-&gt;numpy.array syntax, i wish they were more consistent...<br>
<br>
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.&nbsp; For instance, if I want
a function to operate on every sample, the elegant way would be:<br>
<br>
for vol in nim.data:<br>
&nbsp;&nbsp;&nbsp; some code...<br>
<br>
but I instead must:<br>
<br>
for s in range(nim.timepoints):<br>
&nbsp;&nbsp;&nbsp; if nim.data.ndim==4: # ignoring u,v,w for simplicity here, but this
should be managed too<br>
&nbsp;&nbsp;&nbsp; &nbsp;&nbsp; vol = nim.data[s]<br>
&nbsp;&nbsp;&nbsp; else:<br>
&nbsp;&nbsp;&nbsp; &nbsp;&nbsp; vol = nim.data&nbsp;&nbsp; #Indexing data[s] here would extract a z slice!<br>
<br>
&nbsp;&nbsp;&nbsp; some code with vol...<br>
<br>
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.&nbsp; 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).&nbsp; Perhaps even an
implementation of extend would be useful too so that NiftiImage objects
can be directly concatenated<br>
<blockquote cite="mid:20090112211117.GA16553@runkel" type="cite"><br>
  <pre wrap=""><!---->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.

  </pre>
</blockquote>
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.&nbsp; Even though you've never called
squeeze, the result is the same as if they had been squeezed out of the
reference 7d array.&nbsp; 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.
<blockquote cite="mid:20090112211117.GA16553@runkel" type="cite">
  <pre wrap="">You always know: the last is always X, the one before is always Y, ....

  </pre>
</blockquote>
I've never needed to know what's in shape[-1], i always need to know
what's in shape[0].&nbsp; 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.&nbsp;
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.<br>
<blockquote cite="mid:20090112211117.GA16553@runkel" type="cite">
  <pre wrap=""></pre>
  <blockquote type="cite">
    <pre wrap="">  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]
    </pre>
  </blockquote>
  <pre wrap=""><!---->That should work perfectly.


  </pre>
</blockquote>
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!<br>
<br>
Also, to address your confusion below, I was unclear when I wrote
nim.data[range] above.&nbsp; 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.&nbsp; For clarity, I meant
something more like:<br>
<br>
def getBetas(im, timepoints=None):<br>
&nbsp;&nbsp;&nbsp; nim = niftiimage(im, load=True)&nbsp;&nbsp; # Let's assume 4d, even if only 1
timepoint<br>
&nbsp;&nbsp;&nbsp; if timepoints is None:<br>
&nbsp;&nbsp;&nbsp; &nbsp;&nbsp; return nim.data<br>
&nbsp;&nbsp;&nbsp; if not operater.isSequenceType(timepoints):<br>
&nbsp;&nbsp;&nbsp; &nbsp;&nbsp; timepoints=[timepoints]&nbsp;&nbsp; # 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<br>
&nbsp;&nbsp;&nbsp; return nim.data[timepoints]<br>
<blockquote cite="mid:20090112211117.GA16553@runkel" type="cite">
  <pre wrap=""></pre>
  <blockquote type="cite">
    <pre wrap="">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
    </pre>
  </blockquote>
  <pre wrap=""><!---->No.

  </pre>
  <blockquote type="cite">
    <pre wrap="">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]
    </pre>
  </blockquote>
  <pre wrap=""><!---->Don't get this. Above returns you a 3d volume, but this would return a
4d image -- is that intended?

  </pre>
</blockquote>
yes, because we need to concatenate samples together along the 4th (or
rather, 1st of four) dimension, so everything needs to have that dim.&nbsp;
(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!)<br>
<br>
<blockquote cite="mid:20090112211117.GA16553@runkel" type="cite">
  <blockquote type="cite">
    <pre wrap=""> 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.
    </pre>
  </blockquote>
  <pre wrap=""><!---->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.
  </pre>
  <blockquote type="cite">
    <blockquote type="cite">
      <blockquote type="cite">
        <pre wrap="">from nifti import *
nim=NiftiImage('example4d')
nim.header['dim']
        </pre>
      </blockquote>
    </blockquote>
  </blockquote>
  <pre wrap=""><!---->[4, 128, 96, 24, 2, 1, 1, 1]
  </pre>
  <blockquote type="cite">
    <blockquote type="cite">
      <blockquote type="cite">
        <pre wrap="">m=nim.data.mean(axis=0)
m.shape
        </pre>
      </blockquote>
    </blockquote>
  </blockquote>
  <pre wrap=""><!---->(24, 96, 128)
  </pre>
  <blockquote type="cite">
    <blockquote type="cite">
      <blockquote type="cite">
        <pre wrap="">nim2=NiftiImage(m, nim.header)
nim2.header['dim']
        </pre>
      </blockquote>
    </blockquote>
  </blockquote>
  <pre wrap=""><!---->[3, 128, 96, 24, 1, 1, 1, 1]
  </pre>
  <pre wrap=""><!---->

That is how it should go (IMHO).


Michael

  </pre>
</blockquote>
Yes - but this works because m itself was reduced to 3 dimensions, when
mean squeezed out the singleton on axis=0.&nbsp; 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].&nbsp; The
behavior would be more like:<br>
<br>
&gt;&gt;&gt; r = numpy.random.rand<br>
&gt;&gt;&gt; r = numpy.random.rand(2, 3, 4)<br>
&gt;&gt;&gt; r.shape<br>
(2, 3, 4)<br>
<br>
&gt;&gt;&gt; nim = NiftiImage(r) # As above<br>
&gt;&gt;&gt; nim.data.shape<br>
(2, 3, 4)<br>
&gt;&gt;&gt; nim.header['dim']<br>
[3, 4, 3, 2, 1, 1, 1, 1]<br>
<br>
&gt;&gt;&gt; nim2 = NiftiImage(numpy.asarray([r]))&nbsp;&nbsp;&nbsp; # 4d with 1
timepoint<br>
&gt;&gt;&gt; nim2.header['dim']<br>
[4, 4, 3, 2, 1, 1, 1, 1]<br>
&gt;&gt;&gt; nim2.data.shape<br>
(1, 2, 3, 4)<br>
<br>
&gt;&gt;&gt; nim3 = NiftiImage(numpy.asarray([r]), header=nim.header) #
Fed 3d header<br>
&gt;&gt;&gt; nim3.data.shape<br>
(1, 2, 3, 4)<br>
&gt;&gt;&gt; nim3.header['dim']&nbsp;&nbsp; # Oops! The dims are still
inconsistent with nim and r<br>
[4, 4, 3, 2, 1, 1, 1, 1]<br>
<br>
Now we have a 4d image with 1 timepoint, even if we fed it the correct
header info!&nbsp; 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.&nbsp; 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. <br>
<br>
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.<br>
</body>
</html>