- Generated by
1.13.2
|
fitpack
Modern Fortran library for curve and surface fitting with splines
|
Fitting on non-rectangular domains—disc-shaped regions in the plane and the surface of a sphere—requires coordinate transformations that map the physical domain to a rectangular parameter space where tensor-product B-splines can be applied. These transformations introduce singularities (the origin of a disc, the poles of a sphere) where special continuity conditions must be imposed on the B-spline coefficients to ensure a smooth, single-valued result.
The presentation follows Dierckx (1993), Chapter 11.
A general polar domain is the set of points \( (x, y) \) satisfying \( x^2 + y^2 \leq R(\theta)^2 \), where \( R(\theta) \) is a positive boundary function and \( \theta = \mathrm{atan2}(y, x) \). When \( R(\theta) \) is constant the domain is a circular disc.
FITPACK maps the physical coordinates \( (x, y) \) to normalized polar coordinates \( (u, v) \) via
\[ x = u \, R(v) \cos v, \qquad y = u \, R(v) \sin v \]
with \( u \in [0, 1] \) and \( v \in [-\pi, \pi] \). The radial variable \( u \) runs from 0 (the origin) to 1 (the boundary), and \( v \) is the polar angle. In this parameter space the fitting problem becomes rectangular: find a tensor-product B-spline
\[ s(u, v) = \sum_{i} \sum_{j} c_{ij} \, N_{i,k}(u) \, M_{j,k}(v) \]
that approximates the data in a weighted least-squares sense, subject to a smoothing condition or prescribed knot placement.
| Routine | Input | Boundary |
|---|---|---|
polar (scattered) | Arbitrary \( (x_i, y_i, z_i) \) | General \( R(\theta) \) via function pointer |
pogrid (gridded) | Data on a polar grid \( (u_i, v_j) \) | Constant radius \( R \) |
At \( u = 0 \) all values of the angle \( v \) map to the same physical point—the origin. A naive tensor-product spline in \( (u, v) \) would assign potentially different values to \( s(0, v) \) for different \( v \), producing a discontinuous or non-smooth surface at the centre of the disc.
To obtain a physically meaningful result the B-spline coefficients near \( u = 0 \) must satisfy algebraic constraints that enforce continuity to the desired order.
** \( C^0 \) continuity.** The function value at the origin must be single-valued:
\[ s(0, v) = c_0 \quad \text{(constant for all } v\text{)} \]
This is achieved by constraining all coefficients associated with the first radial B-spline to be equal.
** \( C^1 \) continuity.** In addition to \( C^0 \), the first radial derivative \( \partial s / \partial u \) at \( u = 0 \) must be consistent with a well-defined gradient vector in the \( (x, y) \) plane. Without this constraint, the directional derivative at the origin could depend on the approach angle in an unphysical way. The constraint takes the form of a linear relation among the coefficients of the first two radial B-splines.
** \( C^2 \) continuity.** The second radial derivatives must likewise be consistent with a well-defined Hessian at the origin. This imposes additional linear relations on the first three rows of coefficients.
The bc_continuity_origin parameter selects the continuity order:
bc_continuity_origin | Continuity | Available for |
|---|---|---|
| 0 | \( C^0 \) | Scattered and gridded polar |
| 1 | \( C^1 \) | Scattered and gridded polar |
| 2 | \( C^2 \) (default) | Scattered polar only |
The constraints are incorporated into the least-squares system during fitting, reducing the effective number of free coefficients near \( u = 0 \).
At the outer boundary \( u = 1 \) the following options are available:
For gridded polar fitting (pogrid), additional options control the behaviour at the origin:
z0_exact = .true.) or treated as an additional data point to be approximated.z0_zero_gradient = .true. enforces \( \nabla s(0, 0) = 0 \), creating a flat spot at the centre of the disc.The unit sphere is parameterized by colatitude \( \theta \in [0, \pi] \) and longitude \( \phi \in [0, 2\pi] \). A scalar function on the sphere is represented as a tensor-product B-spline:
\[ s(\theta, \phi) = \sum_{i} \sum_{j} c_{ij} \, N_{i,k}(\theta) \, M_{j,k}(\phi) \]
Because the sphere is periodic in \( \phi \), the spline must satisfy
\[ s(\theta, 0) = s(\theta, 2\pi) \]
together with matching of all derivatives at \( \phi = 0 \) and \( \phi = 2\pi \). FITPACK enforces this by using a periodic knot vector and periodic B-spline basis in the \( \phi \) direction.
| Routine | Input |
|---|---|
sphere (scattered) | Arbitrary \( (\theta_i, \phi_i, r_i) \) on the sphere |
spgrid (gridded) | Data on a \( (\theta_i, \phi_j) \) grid |
The pole problem on the sphere is the direct analogue of the origin problem in polar coordinates. At the north pole ( \( \theta = 0 \)) and the south pole ( \( \theta = \pi \)), all longitudes \( \phi \) correspond to the same physical point. A naive tensor-product spline would allow different values of \( s(0, \phi) \) at different longitudes, producing a discontinuity at the pole.
** \( C^0 \) continuity.** The function value at each pole must be single-valued:
\[ s(0, \phi) = c_{\mathrm{NP}} \quad \text{(constant for all } \phi\text{)}, \qquad s(\pi, \phi) = c_{\mathrm{SP}} \quad \text{(constant for all } \phi\text{)} \]
** \( C^1 \) continuity.** Additionally, the angular derivatives at each pole must be consistent with a well-defined gradient on the sphere. Without this constraint the directional derivative would depend on the approach azimuth in a manner incompatible with the smoothness of the underlying surface.
FITPACK enforces these conditions by imposing linear constraints on the B-spline coefficients near \( \theta = 0 \) and \( \theta = \pi \), reducing the number of free parameters at each pole.
For gridded spherical fitting (spgrid), the north and south poles have independently configurable boundary conditions:
| Setting | Description |
|---|---|
z0 | Prescribed function value at the pole |
exact | If .true., the spline passes exactly through z0 |
continuity | Continuity order at the pole (0 = \( C^0 \), 1 = \( C^1 \)) |
zero_gradient | If .true., enforce vanishing gradient at the pole |
Each pole is configured independently via BC_north_pole and BC_south_pole, allowing asymmetric physical conditions (e.g., a known temperature at the north pole with a free south pole).
polar) accepts a general boundary function \( R(\theta) \) through a function pointer, allowing fitting on non-circular disc-shaped domains. The gridded polar fitter (pogrid) requires a constant radius \( R \).eval_many returns an output array of shape \( (n_\phi, n_\theta) \), not paired points. This matches the tensor-product structure and allows efficient grid evaluation.polar, pogrid, sphere, spgrid) support both smoothing and interpolation modes. In smoothing mode the continuity and boundary constraints are incorporated into the penalized least-squares problem; in interpolation mode they appear as hard constraints.P. Dierckx, Curve and Surface Fitting with Splines, Oxford University Press, 1993, Chapter 11.
1.13.2