Fitting a Surface to Scattered x-y-z Data Points

(1) Using Fit
(2) Using TriangularSurfacePlot
(3) Using ListSurfacePlot3D (non-uniform rectangular grid)

The problem: how to construct a surface which passes through, or nearly through, a set of (scattered) points in 3D, given a bunch of x-y-z data. For the two approaches demonstrated below, we use the same x-y-z data taken from a sample data file 3D.dat:

In[1]:=  xyz = ReadList["3D.dat", {Real, Real, Real}];

The x and y values of these data all lie between -3 and 3.

(1) Using Fit

If a smooth, ``simple'' surface which passes reasonably near the x-y-z data points is sufficient, you can try Fit. You specify the kind of simple functions of x and y to use as building blocks for the surface you seek. Fit finds the linear combination of these functions which best approximates the actual data.

Comments: Finding a good fit may be more of an art than a science. If you use a large set of functions for your basis, you may create a surface which passes through each data point but is suspiciously complicated. Using few basis functions may yield a smoother, simpler surface which only approximates the original data. The goal of getting an elementary function which is smooth and easily expressed, yet hits each data point exactly, may be unattainable. Always choose your basis according to the kinds of functions which are appropriate for your data, whether trigonometric (if periodic), polynomial, rational, exponential, logarithmic, etc.

In the example below I use a 4^th-degree polynomial in x and y to produce a surface. Fit finds the best such polynomial. There is no guarantee that any polynomial is a good choice for this purpose; perhaps i should have used trigonometric or exponentials instead, or in addition to, polynomials. (For some data there really is no good fit.)

The higher the degree of the polynomial, the better it can fit data in the middle; but at the edges of the x-y domain of the original data the polynomial goes wild. If you plot a high-degree polynomial outside its region of fit, you will get huge peaks/valleys which will dwarf the z values of the original data.

In[2]:=  polynomial = Flatten[Table[x^i y^j, {i, 0, 4}, {j, 0, 4}]]

     {1, y, y^2, y^3, y^4, x, x y, x y^2, x y^3, x y^4, x^2,
        x^2 y, x^2 y^2, x^2 y^3, x^2 y^4, x^3, x^3 y, x^3 y^2,
        x^3 y^3, x^3 y^4, x^4, x^4 y, x^4 y^2, x^4 y^3, x^4 y^4}

In[3]:=  f[x_, y_] = Fit[x yz, polynomial, {x, y}]

     -7.71905 - 0.104020 x + 0.408061 x^2 + 0.015997 x^3 - 0.067230 x^4 +
	0.591682 y - 0.357875 x y - 0.262249 x^2 y + 0.042031 x^3 y +
	0.022747 x^4 y + 0.300255 y^2 - 0.064822 x y^2 - 0.307157 x^2 y^2 +
	0.004589 x^3 y^2 + 0.031737 x^4 y^2 - 0.060495 y^3 + 0.046838 x y^3 +
	0.031218 x^2 y^3 - 0.005459 x^3 y^3 - 0.002880 x^4 y^3 - 0.031037 y^4 +
	0.008220 x y^4 + 0.026445 x^2 y^4 - 0.000590 x^3 y^4 - 0.002859 x^4 y^4

In[4]:=  Plot3D[f[x, y], {x, -3, 3}, {y, -3, 3},
	PlotRange -> {-10, -6}, BoxRatios -> {6, 6, 5},
	PlotPoints -> 30, ViewPoint -> {2, 5, 2}];

(2) Using TriangularSurfacePlot

This function uses the data points as nodes for a collection of triangles which build a surface over the convex hull which contains all the x-y values.

In[5]:=  Needs["DiscreteMath`ComputationalGeometry`"]

In[6]:=  TriangularSurfacePlot[xyz, ViewPoint -> {2, 5, 2}];

(3) Non-Uniform Rectangular Grid

If the data is not totally irregular, but the x-y values actually lie on a non-uniform rectangular grid, there is an easier method.

For example, suppose that the x-y-z data consists of z values at each node of a grid where x takes the values {1,3,5,6,7} and y takes the values {1,2,3,7,8} (for each value of x). Then, defining "griddata" as a 5×5 array of x-y-z triplets, e.g.,

In[7]:= griddata = {
{{1,1,-0.0259}, {1,2,-2.0027}, {1,3,-2.6977}, {1,7,-0.6151}, {1,8,-2.3752}},
{{3,1, 1.5260}, {3,2,-0.4508}, {3,3,-1.1458}, {3,7, 0.9368}, {3,8,-0.8234}},
{{5,1, 1.4226}, {5,2,-0.5542}, {5,3,-1.2492}, {5,7, 0.8334}, {5,8,-0.9267}},
{{6,1, 0.5533}, {6,2,-1.4235}, {6,3,-2.1185}, {6,7,-0.0360}, {6,8,-1.7961}},
{{7,1,-0.0786}, {7,2,-2.0555}, {7,3,-2.7504}, {7,7,-0.6679}, {7,8,-2.4280}}
	}

...you can plot the surface with rectangular patches using ListSurfacePlot3D from the Graphics3D package: e.g.,

In[8]:= Needs["Graphics`Graphics3D`"]

In[9]:= ListSurfacePlot3D[ griddata,
                BoxRatios->{3,3,1},ViewPoint->{0.9, 0.1, 0.4}];