Thin Plate Spline — Control Systems/Curve Fitting

Thin-Plate Spline

Control Systems / Curve Fitting

Fits the smoothing thin-plate spline to the last W scattered samples (x, y, z) and reads the surface back at the newest site and at fixed query points:

f(x, y) = Σ_j c_j·φ(|(x, y) − s_j|²) + a_x·x + a_y·y + a₀,   φ(r²) = r²·log(r²)

with the coefficients chosen by the smoothing system (K + λ·I)·c + P·a = z, P'·c = 0, where K_ij = φ(|s_i − s_j|²), P = [x y 1] and λ = (1 − p)/p. This is MATLAB's tpaps([x; y], z, p) read back by fnval, and p means what it means there: 1 interpolates every site, 0 is the least-squares plane, and anything between trades closeness for bending energy.

Ports

• x – the first coordinate of the current site, a scalar.
• y – the second coordinate of the current site, a scalar.
• z – the value measured at that site, a scalar. The block keeps the last W triples (x, y, z) itself.
• zf – the fitted surface at the newest site, f(x, y), [1, 1]: the smoothed version of z. At p = 1 it equals z.
• zq – the fitted surface at each column of Query Points, [Q, 1], Q the number of query points, in their order.

Parameters

• Window Length – W, how many of the most recent sites the surface is fitted to. A whole number from 4 to 24: three sites already determine a plane and leave the spline part nothing to do.
• Smoothing Parameter Choice – where p comes from:
  • Automatic (tpaps default) – recomputed every sample from the window's own sites, by tpaps's own rule when it is called without p: p = 1/(1 + mean(diag(Q₂'·K·Q₂))), where the columns of Q₂ span the vectors orthogonal to P. The default.
  • Given – Smoothing Parameter is used as it stands.
• Smoothing Parameter – p, a single number from 0 to 1, read only when the choice is Given. p = 0 fits the least-squares plane through the window and nothing else, which is what tpaps does with it. At p = 1 the surface interpolates every site – see Notes for what that costs when two sites come close.
• Query Points – where zq reads the surface: a matrix of two rows, one column per point [x; y], 1 to 8 columns – the same shape fnval takes for a bivariate spline. The default, [0; 0], is one point at the origin.
• Sampling Time (s) – zero or less inherits the solver's rate; a positive value runs the block at that period.

Code export

All ten targets: Python, MATLAB, Java, Rust, C, C++, VHDL, Verilog, SystemVerilog and PLC Structured Text. Every target keeps the window itself and, every sample, builds the kernel matrix (W² logarithms), solves the (W+3)×(W+3) bordered system by Gaussian elimination with partial pivoting, and evaluates the surface (W logarithms per point). The window length, the choice of p, a given p and the query points are all structural: they are inlined at export time rather than exposed as tunable parameters, so re-export after changing them.

The three HDL targets are simulation-only: the solve and the logarithm run in real arithmetic and only the ports are Q16.16. A dense solve whose matrix is built from the data every sample, with a logarithm in every entry, has no fixed-point form that would keep its accuracy, so these cores are offered for simulation and not for synthesis.

Simulink bridge

None (Support::None). The thin-plate smoothing spline is a Curve Fitting Toolbox function (tpaps), and that toolbox ships no Simulink library at all, so there is no path a diagram could name. The bridge reports this block rather than dropping it silently, and it therefore has no parity testbench; code export verification still covers it across all ten languages.

Notes

• Stateful, and discrete by nature (setDiscreteOnlyBlock(true)): the window advances once per sample.
• ⚠ Both outputs are 0 whenever the fit is undefined, rather than the arbitrary numbers a near-singular solve would produce:
  • for the first W − 1 samples, until the window is full. The window is not treated as zero-filled: a zero sample would be a site at the origin, and W−1 of them stacked there make a degenerate surface – measured at W = 11, a window with two real samples reaches a condition number of 6e52;
  • when the window's sites are collinear – a surface over a line is not determined, and tpaps refuses it too;
  • when the solve cannot separate the sites: p = 1 with a repeated site, or the automatic p on a window with fewer than four distinct sites, where tpaps's own rule gives p = 1.
  All three are relative tests on the elimination's pivots, so they do not depend on the units of x, y or z.
• ⚠ p = 1 interpolates, and interpolation is ill-conditioned when two sites come close with different values: the surface has to bend hard to pass through both. A small amount of smoothing (p just below 1) cures it; the automatic p is well away from 1 on ordinary data.
• Verified against MATLAB, on this block's own run. Its recorded output was replayed window by window in R2026a – 40 windows of 11 sites per setting – and agrees with fnval(tpaps(...)) at the newest site and at three query points to 2.8e−15 with the automatic p, 8.1e−15 at p = 0.83 and 1.8e−15 at p = 0; the automatic p chose 0.415 to 0.635 over those windows. The elimination here is not MATLAB's QR, but on windows like these the system's condition number stays below 200, which is what keeps the two in agreement.
• The automatic p follows the sites. It is the rule tpaps applies when it is given no p, recomputed from each window; it is not a constant to be read once and reused, and the surface can therefore change its smoothness as the sites move.
• ⚠ It is not Surface Fit or Cubic Spline. This block's surface bends to the data, with no model chosen beforehand; a polynomial surface of a given degree is Surface Fit's job, and a curve in one variable is Cubic Spline's or Spline Fit's.
• Scalar ports. One site per sample; wire one block per surface.
• No state space. The output is nonlinear in the sites, so model reduction correctly declines to merge it.