[SciPy-User] 2D slice of transformed data

[Chris Weisiger]

One thought that occurred to me shortly after leaving work: I know that the
transformation will produce a rectangle (right? There should be no way to
use the transformations specified to turn a rectangle into anything but a
different rectangle). So why not transform the four corners of the
rectangle, then step between them to generate coordinates to use with
map_coordinates? So for example, if I want to generate a slice that has
dimensions 40x512, and the transformed corners are A, B, C, and D, then I
can generate the reverse-transformed grid coordinates as
A, A + AB/40, A + 2AB/40, ...
A + AC/512, A + AC/512 + AB/40, A + AC/512 + 2AB/40, ...
A + 2AC/512, ...

That would require iterating over every pixel in the desired slice (unless
there's some way to vectorize that), but I don't see why it otherwise
wouldn't work. Then just toss those coordinates into map_coordinates and I
have my interpolated slice data.

-Chris

On Wed, Mar 23, 2011 at 5:07 PM, David Baddeley <david_baddeley at yahoo.com.au
> wrote:

> I think you are on to the right idea with using transformed coordinates
> with map_coordinates. I think the easiest way (and fastest way) to get your
> set of transformed coordinates might be something like the following:
>
> assuming Ai is the inverse of your transformation matrix ...
>
> X, Y, Z, T = mgrid[whatever your goal slices are]
>
> Xt = Ai[0,0]*X + Ai[0,1]*Y + Ai[0,2]*Z .....
> Yt = Ai[1,0]*X ........
> Tt = ...
> Tt = ...
>
> this lets you still exploit scipys vector operations on your arrays of
> coordinates, rather than having to separately multiply each tuple with the
> matrix. I'm you can probably write that as a loop to avoid all the redundant
> code. eg:
>
> def transformCoords(coords, A):
>     #coords is a list of coordinate arrays eg [X, Y, Z, T]
>     #A is the transformation matrix
>
>     ret = []
>     for i in range(len(coords)):
>          r_i = 0
>          for j in range(len(coords)):
>                r_j += A[i,j]*coords[j]
>
>          ret.append(r_i)
>
>    return ret
>
> cheers,
> David
>
> ------------------------------
> *From:* Chris Weisiger <cweisiger at msg.ucsf.edu>
> *To:* scipy-user at scipy.org
> *Sent:* Thu, 24 March, 2011 11:00:39 AM
> *Subject:* [SciPy-User] 2D slice of transformed data
>
> 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
>
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20110323/92e0e8d3/attachment.html>