bicubic spline gridded 2d interpolation
C = splin2d(x, y, z, [,spline_type])
:
This function computes a bicubic spline or sub-spline s which interpolates the (xi,yj,zij) points, ie, we have s(xi,yj)=zij for all i=1,..,nx and j=1,..,ny. The resulting spline s is defined by the triplet (x,y,C) where C is the vector (of length 16(nx-1)(ny-1)) with the coefficients of each of the (nx-1)(ny-1) bicubic patches : on [x(i) x(i+1)]x[y(j) y(j+1)], s is defined by :
The evaluation of s at some points must be done by the `interp2d`_ function. Several kind of splines may be computed by selecting the appropriate spline_type parameter. The method used to compute the bicubic spline (or sub-spline) is the old fashionned one ‘s, i.e. to compute on each grid point (xi,yj) an approximation of the first derivatives ds/dx(xi,yj) and ds/dy(xi,yj) and of the cross derivative d2s/dxdy(xi,yj). Those derivatives are computed by the mean of 1d spline schemes leading to a C2 function ( s is twice continuously differentiable) or by the mean of a local approximation scheme leading to a C1 function only. This scheme is selected with the spline_type parameter (see `splin`_ for details) :
:”not_a_knot” this is the default case. : :”periodic” to use if the underlying function is periodic : you must
have z(1,j) = z(nx,j) for all j in [1,ny] and z(i,1) = z(i,ny) for i in [1,nx] but this is not verified by the interface.
:
From an accuracy point of view use essentially the not_a_knot type or periodic type if the underlying interpolated function is periodic.
The natural, monotone, fast (or fast_periodic) type may be useful in some cases, for instance to limit oscillations ( monotone being the most powerful for that).
To get the coefficients of the bi-cubic patches in a more friendly way you can use c = hypermat([4,4,nx-1,ny-1],C) then the coefficient (k,l) of the patch (i,j) (see equation here before) is stored at c(k,l,i,j). Nevertheless the `interp2d`_ function wait for the big vector C and not for the hypermatrix c (note that one can easily retrieve C from c with C=c(:)).
// example 1 : interpolation of cos(x)cos(y)
n = 7; // a regular grid with n x n interpolation points
// will be used
x = `linspace`_(0,2*%pi,n); y = x;
z = `cos`_(x')*`cos`_(y);
C = splin2d(x, y, z, "periodic");
m = 50; // discretisation parameter of the evaluation grid
xx = `linspace`_(0,2*%pi,m); yy = xx;
[XX,YY] = `ndgrid`_(xx,yy);
zz = `interp2d`_(XX,YY, x, y, C);
emax = `max`_(`abs`_(zz - `cos`_(xx')*`cos`_(yy)));
`clf`_()
`plot3d`_(xx, yy, zz, flag=[2 4 4])
[X,Y] = `ndgrid`_(x,y);
`param3d1`_(X,Y,`list`_(z,-9*`ones`_(1,n)), flag=[0 0])
str = `msprintf`_(" with %d x %d interpolation points. ermax = %g",n,n,emax)
`xtitle`_("spline interpolation of cos(x)cos(y)"+str)
// example 2 : different interpolation functions on random data
n = 6;
x = `linspace`_(0,1,n); y = x;
z = `rand`_(n,n);
np = 50;
xp = `linspace`_(0,1,np); yp = xp;
[XP, YP] = `ndgrid`_(xp,yp);
ZP1 = `interp2d`_(XP, YP, x, y, splin2d(x, y, z, "not_a_knot"));
ZP2 = `linear_interpn`_(XP, YP, x, y, z);
ZP3 = `interp2d`_(XP, YP, x, y, splin2d(x, y, z, "natural"));
ZP4 = `interp2d`_(XP, YP, x, y, splin2d(x, y, z, "monotone"));
`xset`_("colormap", `jetcolormap`_(64))
`clf`_()
`subplot`_(2,2,1)
`plot3d1`_(xp, yp, ZP1, flag=[2 2 4])
`xtitle`_("not_a_knot")
`subplot`_(2,2,2)
`plot3d1`_(xp, yp, ZP2, flag=[2 2 4])
`xtitle`_("bilinear interpolation")
`subplot`_(2,2,3)
`plot3d1`_(xp, yp, ZP3, flag=[2 2 4])
`xtitle`_("natural")
`subplot`_(2,2,4)
`plot3d1`_(xp, yp, ZP4, flag=[2 2 4])
`xtitle`_("monotone")
`show_window`_()
// example 3 : not_a_knot spline and monotone sub-spline
// on a step function
a = 0; b = 1; c = 0.25; d = 0.75;
// create interpolation grid
n = 11;
x = `linspace`_(a,b,n);
ind = `find`_(c <= x & x <= d);
z = `zeros`_(n,n); z(ind,ind) = 1; // a step inside a square
// create evaluation grid
np = 220;
xp = `linspace`_(a,b, np);
[XP, YP] = `ndgrid`_(xp, xp);
zp1 = `interp2d`_(XP, YP, x, x, splin2d(x,x,z));
zp2 = `interp2d`_(XP, YP, x, x, splin2d(x,x,z,"monotone"));
// plot
`clf`_()
`xset`_("colormap",`jetcolormap`_(128))
`subplot`_(1,2,1)
`plot3d1`_(xp, xp, zp1, flag=[-2 6 4])
`xtitle`_("spline (not_a_knot)")
`subplot`_(1,2,2)
`plot3d1`_(xp, xp, zp2, flag=[-2 6 4])
`xtitle`_("subspline (monotone)")
Version Description 5.4.0 previously, imaginary part of input arguments were implicitly ignored. .. _interp2d: interp2d.html .. _cshep2d: cshep2d.html .. _linear_interpn: linear_interpn.html .. _splin: splin.html