...

scipy.ndimage.map_coordinates() performs b-spline interpolation of regularly-spaced data

my data are pairs of numbers, like complex numbers. I can't see a way of processing them.

My apologies; I misread your email. The traditional way of interpolating moltivariate data is to do multiple univariate interpolations, as far as I can tell. (E.g. the "thin plate spline" and related literature for defining/manipulating image deformations deals in sparse transforms of (x_old, y_old) -> (x_new, y_new), not unlike what you describe. But all the operations are defined separately for (x_old, y_old) -> x_new and (x_old, y_old) -> y_new, simplifying matters.) Are you thinking that doing the operations "together" (getting spline coefficients simultaneously for the x and y mapping) would/should somehow yield different coefficients than doing them separately? I've certainly never seen anything like that, but I'm far from an expert on the matter. But, as above, from everything I've seen, you can just do the interpolations separately for x and y, and then knit the results together at the end. Zach

...

(spline order 0-5, with several options for boundary conditions). The syntax can seem a bit tricky at first, and you need to watch out for ringing artifacts at sharp transitions (as these are interpolating splines), but it should do the trick.

Zach

Thanks, Zach, but I tried all the routines in interpolate, ndimage and signal, and all of these only seem to use use one-dimensional values. Using the routines in ndimage, I can easily have a 3D array of floats and interpolate at arbitrary 3D coordinates, but this is not what I want. I want multidimensional values, not coordinates. My coordinates are plain 2D x,y coordinates, but the values defined over them are pairs of numbers.

My data would look something like:

(V1,V2) (V1,V2) .... (V1,V2) (V1,V2) (V1,V2) .... (V1,V2) ... (V1,V2) (V1,V2) .... (V1,V2)

(a 2D matrix of pairs)

I'd expect a spline coefficient matrix of the same shape

(C1,C2) (C1,C2) ... (C1,C2) (C1,C2) (C1,C2) ... (C1,C2) ... (C1,C2) (C1,C2) ... (C1,C2)

and, when interpolating at (x,y) I'd like a result

(I0,I1)

(a single pair)

Kay

_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user

Zachary Pincus <zachary.pincus <at> yale.edu> writes:

...

Are you thinking that doing the operations "together" (getting spline

coefficients simultaneously for

the x and y mapping) would/should somehow yield different coefficients than doing them separately? I've certainly never seen anything like that, but I'm far from an expert on the matter. But, as above, from everything I've seen, you can just do the interpolations separately for x and y, and then knit the results together at the end.

I certainly hope that doing the interpolations together and separately should yield precisely the same values :) I really only have performance concerns. I can explain this with a very simple example. Let's suppose I just have two data points P0 and P1 with the values (A0,B0) and (A1,B1). Instead of using splines, let's do a linear interpolation, and let the location of Pi be determined by a single coordinate. The underlying function would be f(x), and let f(x0)=P0 and f(x1)=P1. If we interpolate f(x) at some arbitrary point between x0 and x1, we'd have to calculate ( A0 * (x1-x)/(x2-x1) + A1 * (x-x0)/(x2-x1) , B0 * (x1-x)/(x2-x1) + B1 * (x-x0)/(x2-x1) ) obviously, the only difference between the calculations in the first and second line are the A and B values; the remainder, what could be called the interpolation infrastructure, is precisely the same. So calculating it twice, as would be done if the problem were separated, would be inefficient. It would be a better strategy to calculate the weights first W0 = (x1-x)/(x2-x1) W1 = (x-x0)/(x2-x1) and arrive at the result by ( A0 * W0 + A1 * W1 , B0 * W0 , B1 * W1 ) in short, the code lends itself to vectorization - along the vectors making up individual values. If the values are N-tuples, the larger N is the more time can be saved. Now I do not know if these concerns are relevant in my 'real-world' case, but I assume they are. Naively I would assume that calculating a B-spline at a given location would require determining weights for the coefficients relevant to the calculation. This calculation would be the same for each component value, just as in my simple 1D example. So there should be saving potential here.

Oh and PS, as far as performance concerns about doing things this way? Don't worry about it until it becomes a bottleneck! Profile the code to determine whether it is, and then if so, you might consider re-writing a parallel coefficient-evaluation-loop in Cython or something.

But ndimage is pretty zippy and it may well be that this won't turn out to be

maybe my worries are indeed needless. But I'm doing image analysis and I've got real-time stuff at the back of my head - every milisecond I can save is a good milisecond ;-) the performance hit you fear. For now I'll separate the problem, but I might have to dig deeper. I just had high hopes I could do without separation; after all everything in numpy/scipy is made to work together orthogonally, and when the software informed me it won't even handle complex values (this would be all I need currently) I couldn't quite believe it :( Thanks again for your reply. Kay

Zachary Pincus <zachary.pincus <at> yale.edu> writes:

You could then, as I mentioned, write a Cython module to erect all the "infrastructure" in parallel to do the interpolation,

(gulp) I may need another few years until I'm THAT advaced... I've just learnt a bit of swig and I'm not really so keen on learning yet another interface generator

or perhaps you could even figure out how to vectorize it in pure numpy. (Perhaps that's not hard to do...)

My suspicion.

If you do either, definitely drop a line to the list because I'm sure folks might be interested in the code.

sure will. It should make a difference to all applications along these lines - I mean, how do they go on about RGB images? do each channel separately? Someone must have thought of exploiting the obvious redundancies. Kay

Zachary Pincus <zachary.pincus <at> yale.edu> writes:

...

...

You could then, as I mentioned, write a Cython module

Haha... cython is actually less scary than all that. It's more a

"compiled" to C, and which transparently interfaces with normal

python-like language that can be python code. After adding a few type

annotations, it's really easy to write loops that operate on numpy arrays at C-ish speed.

I had a good look at cython before I opted to use swig for my project - I had a large body of extant C++ code to deal with, and swig is better at dealing with existing stuff (you basically include the original C++ headers and add a bit of collateral code). Cython is more for writing new code. But I'd have to write the cython code to operate on numpy data and fit into the scipy environment. That's where the gulp really came from. I'd have to look at what the current code does and figure out how to write a new piece of scipy to fit in.

...

how do they go on about RGB images? do each channel separately? Someone must have thought of exploiting the obvious redundancies.

Probably, but they might not have been writing in Python! After all, 3x slower doesn't change the big-O complexity, and for most tasks that's often going to be fine. I mean, everything equal it's always better to have faster, more general tools, but there's an API complexity cost, and someone needs to write the code, etc. etc. etc...

Still I feel going that way is right. But maybe I can get lucky elsewhere - I suppose if I look around among the various image processing packages I might find just what I want. In hugin, which is the project I work with, they use vigra, and since vigra uses 'generic' programming I wouldn't be too surprised if the type of vectorization I anticipate is already there. I was just hoping to be able to do it all inside scipy rather than having to pull in yet another dependency. Kay