Surface fitting

Question

I do not need a complete answer but just some advice.

I have a sparse matrix of points in a volume. I know a surface passing by these points exists and this surface is mostly flat and relatively smooth with some small harmonic contents. See the two images below (this surface of the first one is a bit too complicated and de second one is a bit too simple).

I am looking for a method (Laplacian surfaces / Splines / etc.) to fit my points to the best surface of minimal order. Actually I am expecting do describe my surface with less than 30 coefficients.

Any advice of the method I can use (with Matlab or Mathematica) ? If I can get a working example it would be awesome.

EDIT:

Ok, It seems my expectations were too high. In fact, it is very difficult to fit a 2D surface without a model. I recently tried with a pretty simple model:

$ f(x,y) = a_0 + a_1\cdot x^3 + a_2\cdot x^2y + a_3\cdot xy + a_4\cdot y^2x a_5\cdot y^3 + b_1 \sin(b_2\cdot x) + b_3 \cos(b_4\cdot x) + c_1 \sin(c_2\cdot x) + c_3 \cos(c_4 \cdot x)$

With this one I can fit my surface with a least squares approach and normal equations. The result is not bad but it would be better to find a "smart" algorithm that can build the output function with an iterative method. I still haven't find anything like this.

Answer 1

Along with Kriging using sparse matrix techniques, you could also look up literature on sparse approximation (aka sparse decomposition) methods or sparse Gaussian graphical models.

Answer 2

I would first look into the thin plate splines (TPS), which are widely available in MATLAB. They are simple to use, and were found to perform fairly good in many cases. Here is the documentation: http://www.mathworks.de/de/help/curvefit/tpaps.html?refresh=true

More detailed information is here: http://user.engineering.uiowa.edu/~aip/papers/bookstein-89.pdf

Answer 3

I guess you already tried the matlab "fit" function.

f = fit( [x, y], z, 'poly23' )

More infos on this function here: matlab fit .

Answer 4

Looks like it would be a good candidate for Gaussian Process Regression (also known as Kriging).

Edit: I'll explain a little more. Consider the space $X$ of all functions on a grid with $N$ nodes, $X := \mathbb{R}^N$. For example, given a vector $x \in X$, the first component $x_1$ is the value of the function at the first gridpoint, $x_2$ is the value at the second gridpoint, and so forth.

The basic idea is to put a gaussian prior distribution on $X$, $$\pi_\text{prior} \sim e^{-x^T C^{-1} x},$$ and consider making measurements $y_k = x_k$ of the function at some nodes of interest $k$.

We can combine the prior information and the measurement information to construct a posterior distribution on $X$ using Bayes' theorem. Some care must be taken since the likelihood is an improper distribution, but it can be done.

Then one takes the mean of the posterior to get a good function that fits the data. This mean function is equivalent to what you get from Kriging.

---

Entries of the prior covariance, $C_{ij}$, measure how similar you expect the function to be at nodes $i$ and $j$. A couple typical examples are,

• Exponential decay covariance: The similarity between the function at two points decreases as the points get further apart with decay rate $q$ (eg., $q=2$ is the most common) and characteristic length $\rho$: $$C_{ij} \propto e^{-\frac{||x_i - x_j||^q}{\rho}}.$$ This matrix is a dense $n^2$-by-$n^2$ matrix, which may present some computational difficulties.
• Laplacian covariance: The covariance falls off in the same manner as the greens function of an elliptic PDE, $\Delta^q C(x-x_0) = \delta_x$: $$C \propto \Delta^{-q},$$ where $\Delta$ is the (graph) Laplacian. This covariance is particularly nice since applying $\Delta^{-q}$ can be done efficiently with multigrid. So, you can compute the action of a dense decaying covariance matrix whenever needed (or it's square root), except without actually forming the dense matrix.