2D smoothing

Syntax: | `wout = BIVSMOOTH(x,y,z,mx,my)`

This function fits a smooth surface of a single-valued bivariate function z = z(x,y) to a set of input data points given at input grid points in an x-y plane. It generates a set of output grid points by equally dividing the x and y coordinates in each interval between a pair of input grid points, interpolates the z value for the x and y values of each output grid point, and generates a set of output points consisting of input data points and the interpolated points. The method is based on a piece-wise function composed of a set of bicubic polynomials in x and y. Each polynomial is applicable to a rectangle of the input grid in the x-y plane. Each polynomial is determined locally.

The first two parameters are vectors. Vector `x` contains the x-coordinates of the input grid points, in ascending or descending order. Vector `y` contains the y-coordinates of the input grid points, in ascending or descending order. Both `x` and `y` must be monotonic. The third parameter is a matrix, `z`, which contains the values of the function at the input grid points, `z[i][j]` is the data value at `(x[i],y[j])`. The last two parameters are scalars. The fourth parameter, `mx`, is the number of subintervals between each pair of input grid points in the x direction, and must be at least 2. The fifth parameter, `my`, is the number of subintervals between each pair of input grid points in the y direction, and must be at least 2. The result of the function is a matrix, `wout`, which has `my*(LEN(y)-1)+1` rows and `mx*(LEN(x)-1)+1` columns. The first row of the matrix, starting in column 2, contains the x coordinates of the output values. It can be extracted into a vector, `u`, with the following:

`u = wout[1,2:VLEN(wout)[2]]`

The first column of the matrix, starting in row 2, contains the y coordinates of the output values. It can be extracted into a vector, `v`, with the following:

`v = wout[2:VLEN(wout)[1],1]`

The smoothed values at the `(u,v)` locations can be extracted into a matrix, `w`, with the following:

`w = wout[2:VLEN(wout)[1],2:VLEN(wout)[2]]`

For example, suppose `x` and `y` are vectors and `m` is a data matrix with `LEN(y)` rows and `LEN(x)` columns:

     nx = 5 ! number of subdivisions in x
     ny = 4 ! number of subdivisions in y
     ww = bivsmooth(x,y,m,nx,ny)
     u = ww[1,2:VLEN(ww)[2]]
     v = ww[2:VLEN(ww)[1],1]
     w = ww[2:VLEN(ww)[1],2:VLEN(ww)[2]]

Algorithm derived from an article by Hiroshi Akima, Communications of the ACM, volume 17, number 1, January 1974, pp. 26-31.

Smoothing