I am looking for a way to fit a surface given a set of measured data $(x, y) \mapsto z$. A typical example would consist of anywhere between $10$ and $30$ measurements spread evenly over a disc. Until now it was sufficient to fit a plane through these points with a least square error fit. This model is not sufficient anymore since the accuracy has to be increased. The first thing we tried was to move to more general polynomial fitting. However, low order polynomials ($2$, $3$, $4$) cannot model the surface well enough in general and with higher orders we quickly seem to run into the "Runge phenomenon" that the global fitting error is hard to control.
Therefore I'm currently thinking in the direction of triangulation and piecewise fitting. I can imagine that even piecewise linear fitting could be an option. I'm not sure how to proceed from here however. Are there ways to construct a simple (low count) triangulation and piecewise fit? (The triangulation vertices may be located anywhere and need not coincide with the locations of the measurements.)
Any (other) suggestions are welcome!