[SciPy-User] Bivariate Spline Surface Fitting

Thu Nov 26 23:03:04 EST 2009

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

knots should only specify the point of x and y not all grid points, I
added an s to play with the border values following the example in the
tests. most of it just trial and error, since I don't have a good
example of what I should get as a result

with the following knots, it finishes without warning and errors, and
I get some numbers back that might be reasonable.

  s = 1.1
  yknots = np.linspace(ymin+s,ymax-s,10)
  xknots = np.linspace(xmin+s,xmax-s,10)

Some good examples for the use of the different options in the spline
classes would be nice.

the docs are still pretty bad, but there is:

473	        Input:
474	          x,y,z  - 1-d sequences of data points (order is not
475	                   important)
476	          tx,ty  - strictly ordered 1-d sequences of knots
477	                   coordinates.
478	        Optional input:
479	          w          - positive 1-d sequence of weights
480	          bbox       - 4-sequence specifying the boundary of
481	                       the rectangular approximation domain.
482	                       By default, bbox=[min(x,tx),max(x,tx),
483	                                         min(y,ty),max(y,ty)]
484	          kx,ky=3,3  - degrees of the bivariate spline.
485	          eps        - a threshold for determining the effective rank
486	                       of an over-determined linear system of
487	                       equations. 0 < eps < 1, default is 1e-16.

Josef

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