[SciPy-User] Bivariate Spline Surface Fitting

[josef.pktd at gmail.com]

On Thu, Nov 26, 2009 at 4:34 PM, Peter Combs <peter.combs at berkeley.edu> wrote:
> Hi all,
> I have localization data in 2 color channels that should agree with each other, but in practice, they don't to the level we want. I thought I'd try doing a straight polynomial least squares fit, and while that gives better registration between the two, I'm still not to the level I want. My next thought was a spline fit, so I'm trying to make two least-squares bivariate spline fits: one for taking (x,y) to x', and one for taking (x,y) to y'.
>
>
> import scipy.interpolate as interp
> ...
> def makeLSQspline(xl, yl, xr, yr):
>   """docstring for makespline"""
>
>   xmin = xr.min()-1
>   xmax = xr.max()+1
>   ymin = yr.min()-1
>   ymax = yr.max()+1
>   n = len(xl)
>
>   print "xrange: ", xmin, xmax, '\t', "yrange: ", ymin, ymax
>
>   yknots, xknots = mgrid[ymin:ymax:10j, xmin:xmax:10j]   # Makes an 11x11 regular grid of knot locations
>
>   xspline = interp.LSQBivariateSpline(xr, yr, xl, xknots.flat, yknots.flat)
>   yspline = interp.LSQBivariateSpline(xr, yr, yl, xknots.flat, yknots.flat)
>
>   def mapping(xr, yr):
>      xl = xspline.ev(xr, yr)
>      yl = yspline.ev(xr, yr)
>      return xl, yl
>   return mapping
>
>
> I have a "Registration Error" function which calculates a mapping for all but the ith point, then plugs that point into the mapping and finds the difference between the predicted value and the known value. For the 2nd order polynomial fit, I get a mean registration error around 7nm, but for the spline fitting using the function above, the mean error is more like 20,000nm. Which (along with all the random junk that gets spit out, such as
> /Library/Frameworks/Python.framework/Versions/5.1.0/lib/python2.5/site-packages/scipy/interpolate/fitpack2.py:498: UserWarning:
> Error on entry, no approximation returned. The following conditions
> must hold:
> xb<=x[i]<=xe, yb<=y[i]<=ye, w[i]>0, i=0..m-1
> If iopt==-1, then
> xb<tx[kx+1]<tx[kx+2]<...<tx[nx-kx-2]<xe
> yb<ty[ky+1]<ty[ky+2]<...<ty[ny-ky-2]<ye

(I'm out the door, so briefly)

I had this error message before when I didn't set up the knots
correctly, but I only tried the UnivariateSplines. It took me a while
to figure out how the endpoints of the knots are supposed to be
defined. If I remember correctly (?), then the endpoints where not
supposed to be included or something like this.

Josef


> warnings.warn(message)
>
> about 1 copy of this (or something similar) per call of makeLSQSpline:
> tx= 0.00000000000000 0.00000000000000 0.00000000000000 -264.022095089756 337.266858978529 3468.05790461355 6598.84895024857 9729.63999588358 12860.4310415186 15991.2220871536 19122.0131327886 22252.8041784237 25383.5952240587 28514.3862696937 337.266858978529 3468.05790461355 6598.84895024857 9729.63999588358 12860.4310415186 15991.2220871536 19122.0131327886 22252.8041784237 25383.5952240587 28514.3862696937 337.266858978529 3468.05790461355 6598.84895024857 9729.63999588358 12860.4310415186 15991.2220871536 19122.0131327886 22252.8041784237 25383.5952240587 28514.3862696937 337.266858978529 3468.05790461355 6598.84895024857 9729.63999588358 12860.4310415186 15991.2220871536 19122.0131327886 22252.8041784237 25383.5952240587 28514.3862696937 337.266858978529 3468.05790461355 6598.84895024857 9729.63999588358 12860.4310415186 15991.2220871536 19122.0131327886 22252.8041784237 25383.5952240587 28514.3862696937 337.266858978529 3468.05790461355 6598.84895024857 9729.63999588358
>  12860.4310415186 15991.2220871536 19122.0131327886 22252.8041784237 25383.5952240587 28514.3862696937 337.266858978529 3468.05790461355 6598.84895024857 9729.63999588358 12860.4310415186 15991.2220871536 19122.0131327886 22252.8041784237 25383.5952240587 28514.3862696937 337.266858978529 3468.05790461355 6598.84895024857 9729.63999588358 12860.4310415186 15991.2220871536 19122.0131327886 22252.8041784237 25383.5952240587 28514.3862696937 337.266858978529 3468.05790461355 6598.84895024857 9729.63999588358 12860.4310415186 15991.2220871536 19122.0131327886 22252.8041784237 25383.5952240587 28514.3862696937 337.266858978529 3468.05790461355 6598.84895024857 9729.63999588358 12860.4310415186 15991.2220871536 19122.0131327886 22252.8041784237 25383.5952240587 28514.3862696937 28778.4083647834 0.00000000000000 0.00000000000000 0.00000000000000
> tx= 0.00000000000000 0.00000000000000 0.00000000000000 -264.022095089756 337.266858978529 3468.05790461355 6598.84895024857 9729.63999588358 12860.4310415186 15991.2220871536 19122.0131327886 22252.8041784237 25383.5952240587 28514.3862696937 337.266858978529 3468.05790461355 6598.84895024857 9729.63999588358 12860.4310415186 15991.2220871536 19122.0131327886 22252.8041784237 25383.5952240587 28514.3862696937 337.266858978529 3468.05790461355 6598.84895024857 9729.63999588358 12860.4310415186 15991.2220871536 19122.0131327886 22252.8041784237 25383.5952240587 28514.3862696937 337.266858978529 3468.05790461355 6598.84895024857 9729.63999588358 12860.4310415186 15991.2220871536 19122.0131327886 22252.8041784237 25383.5952240587 28514.3862696937 337.266858978529 3468.05790461355 6598.84895024857 9729.63999588358 12860.4310415186 15991.2220871536 19122.0131327886 22252.8041784237 25383.5952240587 28514.3862696937 337.266858978529 3468.05790461355 6598.84895024857 9729.63999588358
>  12860.4310415186 15991.2220871536 19122.0131327886 22252.8041784237 25383.5952240587 28514.3862696937 337.266858978529 3468.05790461355 6598.84895024857 9729.63999588358 12860.4310415186 15991.2220871536 19122.0131327886 22252.8041784237 25383.5952240587 28514.3862696937 337.266858978529 3468.05790461355 6598.84895024857 9729.63999588358 12860.4310415186 15991.2220871536 19122.0131327886 22252.8041784237 25383.5952240587 28514.3862696937 337.266858978529 3468.05790461355 6598.84895024857 9729.63999588358 12860.4310415186 15991.2220871536 19122.0131327886 22252.8041784237 25383.5952240587 28514.3862696937 337.266858978529 3468.05790461355 6598.84895024857 9729.63999588358 12860.4310415186 15991.2220871536 19122.0131327886 22252.8041784237 25383.5952240587 28514.3862696937 28778.4083647834 0.00000000000000 0.00000000000000 0.00000000000000
>
> makes me think something isn't quite right. Any guesses what's going on? I have ~3400 data points, roughly evenly spread out over a 28,000nm x 23,000nm grid.
>
> On another note, how well is this going to scale up? If I end up collecting hundreds of thousands to low-millions of points, does spline fitting go as O(n^2), or more like O(n)? The error registration function runs as O(n*O(fitting)), and takes around 5 seconds now, so O(N) spline fitting is fine, about an hour run time total, but O(n^2) is very much not.
> Peter Combs
> peter.combs at berkeley.edu
>
>
>
>
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>