[Nipy-devel] Affine orientation ordering - please review and test

Jonathan Taylor jonathan.taylor@stanford....
Mon Jan 18 23:48:41 CST 2010


Hi Matthew,



> What I was trying to do was answer Dav, Gael's and Benjamin's need to
> be able to easily generate an image for which the _data array_  had
> been rearranged in canonical (XYZ) order.   I think they wanted it
> pull it out to an image viewer without needing to resample.   That's
> now in nipy.io.imageformats.as_closest_xyz  . That is - I'm going to
> use the orientations to change the data matrix directly by reversing
> axes and transposing - hence the apply_orientations routine.
>

I guess I didn't understand what apply_orientations routine was supposed to
do, but after this paragraph I think I see.
Effictively you want to slice arr with possibly some flips. So, if you had
an image
of shape (20,30,40) and you found out you had to swap axes 0 and 1, with
axis 1 having been flipped, then apply_orientations should return

arr[:,::-1].transpose([1,0,2])

That wasn't clear to me at first glance, my apologies. My confusion may
explain why you're frazzled
because I completely missed the ball on what this code was supposed to do.

If this is what you intended it to do, then you can ignore
my comments about the matrices because if that's the ultimate goal, then the
matrices are not necessary.



> > 3) There is really no need to assume 3-dimensions here. I know this will
> > annoy people, but io_orientation could take an (n+1)x(m+1) affine (n>m)
> and
> > return an (m+1)x(n+1) affine which is no longer necessarily square but
> has
> > only 0 or +-1 entries. Then, apply_orientation gives an (m+1)x(m+1)
> affine
> > matrix which represents the projection of the original n coordinates of
> the
> > affine onto an "optimal" (in some sense) m coordinate plane.
>
>
I think I have it extended to general dimensions (patch attached).  One test
was failing for me, namely one based on

In [16]: io_orientation(np.array([[1,1,0,0],[1,1,0,0],[0,0,1,0],[0,0,0,1]]))
Out[16]:
array([[ 1.,  1.],
       [ 0.,  1.],
       [ 2.,  1.]])

It fails because, in the general case, I took R to have the same rank as
affine.
In the original code for this case, because you didn't zero out essentially
0 eigenvalues in computing R from RS
your matrix R has rank 3 while affine has rank 2. With R as you wrote it,
you are getting a
false sense of what is happening because your matrix R is always full rank
and can be quite far from "affine". Also,
if affine is not full rank, than R as originally computed is not
mathematically unique and possibly numerically instable.
The little script "perturb_svd.py" illustrates this.

I changed the test so that all the matrices have rank 3, then all tests
pass.

-- 
Jonathan Taylor
Dept. of Statistics
Sequoia Hall, 137
390 Serra Mall
Stanford, CA 94305
Tel:   650.723.9230
Fax:   650.725.8977
Web: http://www-stat.stanford.edu/~jtaylo<http://www-stat.stanford.edu/%7Ejtaylo>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/nipy-devel/attachments/20100118/5ec56a03/attachment.html 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: orientation.patch
Type: application/octet-stream
Size: 13694 bytes
Desc: not available
Url : http://mail.scipy.org/pipermail/nipy-devel/attachments/20100118/5ec56a03/attachment.obj 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: perturb_svd.py
Type: application/octet-stream
Size: 1636 bytes
Desc: not available
Url : http://mail.scipy.org/pipermail/nipy-devel/attachments/20100118/5ec56a03/attachment-0001.obj 


More information about the Nipy-devel mailing list