[SciPy-User] how can I create B-splines of mutidimensional values?

Kay F. Jahnke _kfj@yahoo....
Mon Nov 28 11:08:50 CST 2011


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.

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 ;-)

> But ndimage is pretty zippy and it may well be that this won't turn out to be
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





More information about the SciPy-User mailing list