Scattered Surface Fitting Tutorial

This tutorial covers bivariate spline fitting to scattered data with fitpack_surface: constructing a tensor-product B-spline \( z = s(x, y) \) from irregularly sampled \( (x_i, y_i, z_i) \) triples, controlling smoothness, and evaluating derived quantities.

Use fitpack_surface whenever your observations are not aligned on a rectangular grid. If the data do lie on a regular grid, prefer fitpack_grid_surface instead (see Grid and Parametric Surface Tutorial).

Mathematical Background

Given \( m \) scattered observations \( (x_i, y_i, z_i) \) with positive weights \( w_i \), FITPACK determines a tensor-product B-spline

\[ s(x, y) = \sum_{i=1}^{n_x - k_x - 1} \sum_{j=1}^{n_y - k_y - 1} c_{ij} \, B_i^{k_x}(x) \, B_j^{k_y}(y) \]

that minimizes the thin-plate smoothing energy subject to a residual constraint:

\[ \sum_{i=1}^{m} \bigl( w_i \, (z_i - s(x_i, y_i)) \bigr)^2 \leq S \]

where \( S \) is the user-supplied smoothing parameter. Knot positions and the number of knots \( n_x, n_y \) are chosen automatically by the surfit algorithm.

The default spline degree is bicubic ( \( k_x = k_y = 3 \)). Degrees 1 through 5 are supported.

See also

Dierckx, Ch. 5, §5.3 (pp. 85–98)

Basic Example

The code below fits a spline to 200 quasi-randomly distributed samples of \( f(x,y) = \sin(\pi x) \cos(\pi y) \) on \( [0,1]^2 \).

program surface_demo
    use fitpack, only: fitpack_surface, fp_real, fp_flag, &
                       fitpack_success, fitpack_message, zero, one, half, pi
    implicit none
 
    integer, parameter :: m = 200
    real(FP_REAL) :: x(m), y(m), z(m)
    type(fitpack_surface) :: surf
    integer(FP_FLAG) :: ierr
    integer :: i
 
    ! Quasi-random sampling (golden-ratio sequences)
    do i = 1, m
        x(i) = mod(i * 0.6180339887_fp_real, one)
        y(i) = mod(i * 0.3247179572_fp_real, one)
        z(i) = sin(pi * x(i)) * cos(pi * y(i))
    end do
 
    ! Fit with smoothing = m (a practical starting value for unit weights)
    ierr = surf%new_fit(x, y, z, smoothing=real(m, fp_real))
    if (.not. fitpack_success(ierr)) error stop fitpack_message(ierr)
 
    print '(a,i0,a,i0)', 'Knots: ', surf%knots(1), ' x ', surf%knots(2)
    print '(a,es10.3)',   'Residual: ', surf%mse()
 
    call surf%destroy()
end program

surf\new_fit(x, y, z, w, smoothing, order) stores the data, allocates workspace, and performs a smoothing fit in a single call. It returns an integer(FP_FLAG) error code.

Features

Point evaluation

Evaluate the fitted surface at one or many scattered points:

! Scalar
z_val = surf%eval(x_val, y_val, ierr)
 
! Array of scattered points
z_arr = surf%eval(x_arr, y_arr, ierr)   ! size(z_arr) == size(x_arr)

Both surf\eval variants call the unscattered spline evaluator bispeu, so the input points need not be sorted.

Grid evaluation

Evaluate on every node of a rectangular output grid:

! z_grid(j,i) = s(x_grid(i), y_grid(j))
z_grid = surf%eval_ongrid(x_grid, y_grid, ierr)

The result has shape (size(y_grid), size(x_grid)).

Partial derivatives

Compute \( \partial^{n_x+n_y} s / \partial x^{n_x} \partial y^{n_y} \) at scattered points or on a grid:

! ds/dx at a single point
dzdx = surf%dfdx(x0, y0, 1, 0, ierr)
 
! d2s/dxdy at many scattered points
d2 = surf%dfdx(x_arr, y_arr, 1, 1, ierr)
 
! ds/dy on a grid
dz_grid = surf%dfdx_ongrid(x_grid, y_grid, 0, 1, ierr)

The derivative order arguments nux, nuy must satisfy \( 0 \leq n_x < k_x \) and \( 0 \leq n_y < k_y \).

Integration

Integrate the surface over a rectangular sub-domain \( [x_1, x_2] \times [y_1, y_2] \):

vol = surf%integral([x1, y1], [x2, y2])

The two arguments are 2-element arrays lower = [x_lo, y_lo] and upper = [x_hi, y_hi].

Cross-sections

Extract a one-dimensional profile by fixing one coordinate:

type(fitpack_curve) :: profile
 
! Profile f(y) = s(x0, y)  (along_y = .true.)
profile = surf%cross_section(x0, .true., ierr)
 
! Profile g(x) = s(x, y0)  (along_y = .false.)
profile = surf%cross_section(y0, .false., ierr)

The returned fitpack_curve is a full 1D spline that can be evaluated, differentiated, and integrated independently.

Derivative spline

Build a new fitpack_surface whose evaluation gives a partial derivative of the original:

type(fitpack_surface) :: dsdx_surf
 
dsdx_surf = surf%derivative_spline(1, 0, ierr)   ! ds/dx

The resulting surface has reduced degree \( (k_x - n_x, \; k_y - n_y) \).

Smoothing Sweep

The smoothing parameter \( S \) controls the trade-off between fidelity and smoothness. Smaller \( S \) forces more knots and a closer fit; larger \( S \) yields a smoother approximation.

print '(a)', '     S       knots_x  knots_y  residual'
do i = 1, 4
    ierr = surf%new_fit(x, y, z, smoothing=10.0_fp_real**(3 - i))
    if (fitpack_success(ierr)) then
        print '(es10.1, 2i9, es12.3)', &
              surf%smoothing, surf%knots(1), surf%knots(2), surf%mse()
    end if
end do

Special cases:

• ** \( S = 0 \)** — interpolation (the spline passes through every data point). Use surf\interpolate().
• **Large \( S \)** — very smooth, few knots, possibly large residual.
• ** \( S = m \)** (number of data points, unit weights) — a practical starting value.

After any fit, surf\mse() returns the weighted sum of squared residuals \( f_p \), and surf\knots returns a 2-element integer array \( [n_x, n_y] \).

Fitting Modes

The same three strategies available for curves work for surfaces:

Method	Call	Description
Automatic knots	surf\fit(smoothing=s)	Knots chosen to satisfy the smoothing constraint
Interpolation	surf\interpolate()	Passes through all data points ( \( S = 0 \))
Least-squares	surf\least_squares()	User-supplied knots, no smoothing constraint

Error Handling

All routines return an integer(FP_FLAG) error code:

use fitpack, only: fitpack_ok, fitpack_success, fitpack_message
 
if (.not. fitpack_success(ierr)) then
    print *, 'Error: ', fitpack_message(ierr)
end if

Flag	Meaning
FITPACK_OK (0)	Success
FITPACK_INTERPOLATING_OK (-1)	Interpolating spline returned
FITPACK_S_TOO_SMALL (1)	Smoothing too small for the knot budget
FITPACK_MAXIT (2)	Maximum iterations reached
FITPACK_INPUT_ERROR (10)	Invalid input data

Complete Example

A self-contained program demonstrating all of the features above is provided in examples/example_surface.f90. Build and run it with:

fpm run --example example_surface