[SciPy-User] 2D slice of transformed data

[Chris Weisiger]

I just wanted to let you all know that I have this working now, and while
it's moderately slow when it has to calculate three slices (of sizes
2x512x512, 2x512x50, and 2x50x512 in my test case), it's not so bad as to be
unusable -- and my users have more powerful computers than I do.

...hm, for that matter, I should be able to split off the three slice
calculations into separate threads...

Anyway, I put my solution online on the off-chance that someone finds it
useful:

http://paste.ubuntu.com/587110/

I can't guarantee that the comments are especially germaine, as this took a
lot of fiddling and I haven't gone through to clean everything up yet, but
the basic flow should be intelligible:

 * Figure out which axes the slice cuts across, and generate an array of the
appropriate shape to hold the results.
 * Create an array of similar size augmented with a length-4 dimension. This
array holds XYZ coordinates for each pixel in the slice; the 4th index holds
a 1 (so that we can use a 4x4 affine transformation matrix to do rotation
and offsets in the same pass). For example, an XY slice at Z = 5 would look
something like this:
[  [0, 0, 5]  [0, 1, 5]  [0, 2, 5] ...
[  [1, 0, 5]  ...
[  [2, 0, 5]
[  ...
[
 * Subtract the XYZ center off of the coordinates so that when we apply the
rotation transformation, it's done about the center of the dataset instead
of the corner.
 * Multiply the inverse transformation matrix by the coordinates.
 * Add the center back on.
 * Chop off the dummy 1 coordinate, reorder to ZYX, and prepend the time
dimension.
 * Pass the list of coordinates off to numpy.map_coordinates so it can look
up actual pixel values.

Thanks again for your help, everyone!

-Chris

On Wed, Mar 23, 2011 at 3:00 PM, Chris Weisiger <cweisiger at msg.ucsf.edu>wrote:

> In preface, I'm not remotely an expert at array manipulation here. I'm an
> experienced programmer, but not an experienced *scientific* programmer. I'm
> sure what I want to do is possible, and I'm pretty certain it's even
> possible to do efficiently, but figuring out the actual implementation is
> giving me fits.
>
> I have two four-dimensional arrays of data: time, Z, Y, X. These represent
> microscopy data taken of the same sample with two different cameras. Their
> views don't quite match up if you overlay them, so we have a
> three-dimensional transform to align one array with the other. That
> transformation consists of X, Y, and Z translations (shifts), rotation about
> the Z axis, and equal scaling in X and Y -- thus, the transformation has 5
> parameters. I can perform the transformation on the data without difficulty
> with ndimage.affine_transform, but because we typically have hundreds of
> millions of pixels in one array, it takes a moderately long time. A
> representative array would be 30x50x512x512 or thereabouts.
>
> I'm writing a program to allow users to adjust the transformation and see
> how well-aligned the data looks from several perspectives. In addition to
> the traditional XY view, we also want to show XZ and YZ views, as well as
> kymographs (e.g. TX, TY, TZ views). Thus, I need to be able to show 2D
> slices of the transformed data in a timely fashion. These slices are always
> perpendicular to two axes (e.g. an XY slice passing through T = 0, Z = 20,
> or a TZ slice passing through X = 256, Y = 256), never diagonal. It seems
> like the fast way to do this would be to take each pixel in the desired
> slice, apply the reverse transform, and figure out where in the original
> data it came from. But I'm having trouble figuring out how to efficiently do
> this.
>
> I could construct a 3D array with shape (length of axis 1), (length of axis
> 2), (4), such that each position in the array is a 4-tuple of the
> coordinates of the pixel in the desired slice. For example, if doing a YX
> slice at T = 10, Z = 20, the array would look like [[[10, 20, 0, 0], [10,
> 20, 1, 0], [10, 20, 2, 0], ...], [[10, 20, 0, 1], 10, 20, 1, 1], ...]]. Then
> perhaps there'd be some way to efficiently apply the inverse transform to
> each coordinate tuple, then using ndimage.map_coordinates to turn those into
> pixel data. But I haven't managed to figure that out yet.
>
> By any chance is this already solved? If not, any suggestions / assistance
> would be wonderful.
>
> -Chris
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20110329/8c971b07/attachment.html>