# Basic 3D Interpolation Failing

I've seen a few other posts asking about the "unstructured grid" error in Interpolation. Here's a minimal (non)working example of the issue:

z = RandomReal[{-1, 1}, 50];
xy = RandomPoint[Rectangle[{0, 0}, {1, 1}], 50];


One way around this is to use InterpolatingPolynomial like this:

data = Thread[{xy, z}];
ip = InterpolatingPolynomial[data, {x, y}];
Show[Plot3D[ip, {x, 0, 1}, {y, 0, 1}],
Graphics3D[Point[Flatten /@ data]]]


However, it doesn't work well in practice as the polynomial orders are uncontrollable (too high):

Are there any viable alternatives for successfully getting quadratic InterpolationOrder in Interpolation?

• A workaround: Perhaps you can create a structured mesh using the FEMfunctionality and ElementMeshInterpolation[] for quadratuic elements? – Ulrich Neumann Dec 19 '18 at 11:01

You can use interpolation based on radial basis functions. To demonstrate the method method see the following example.

data = GeoElevationData[GeoDisk[Entity["Mountain", "MountEverest"], Quantity[3, "Miles"]]];
ListPlot3D[data, MeshFunctions -> {#3 &}, Mesh -> 40]


To make the computation faster I downsample the data

smallData = ArrayResample[QuantityMagnitude[data], {100, 100}];


Now we select 30 random positions as sample points that we are going to use for the interpolation:

randomx = RandomInteger[{1, 100}, 30];
randomy = RandomInteger[{1, 100}, 30];
stuetzstellen = Table[{randomx[[i]], randomy[[i]], QuantityMagnitude[smallData[[randomy[[i]], randomx[[i]]]]]}, {i,1, Length[randomx]}];
Show[{
ListPlot3D[smallData, MeshFunctions -> {#3 &}, Mesh -> 40],
ListPointPlot3D[stuetzstellen, PlotStyle -> PointSize -> 0.03]}]


Now we can define the RBF. There are many possible choices. You can experiment with other functions.

rbf[x_] /; (x > 0) := N[x^2 Log[x]]
rbf[x_] /; (x == 0) := 0
distances[{x_, y_}, data_] := EuclideanDistance[{x, y}, #] & /@ data[[All, {1, 2}]]


The basis for the interpolation is the distance matrix for our sample points:

distMat = DistanceMatrix[stuetzstellen[[All, {1, 2}]]];
coeffs = LinearSolve[Map[rbf, distMat, {2}], stuetzstellen[[All, 3]]]


Now we can compute the interpolated mountain:

Show[{
Plot3D[coeffs.rbf /@ distances[{x, y}, stuetzstellen], {x, 0, 100}, {y, 0, 100}],
ListPlot3D[smallData, MeshFunctions -> {#3 &}, Mesh -> 40, PlotStyle -> Directive[Red, Opacity[0.5]]],
ListPointPlot3D[stuetzstellen, PlotStyle -> PointSize -> Large]}]


The more sample points we have the better will be the approximation. This can also be expanded to higher dimensions.

• This is a nice answer. – user21 Dec 20 '18 at 6:32

The reason only first order interpolation is possible on unstructured grids is because it is not at all clear how to construct the mesh from the unstructured data given. Additional information is necessary; if you have that you can construct a higher order interpolating function.

First, we construct a Delaunay triangulation of the coordinates given.

Needs["NDSolveFEM"]
m = ToElementMesh[xy]
(*m["Wireframe"]*)


Next, we convert that mesh to a second order mesh.

m2 = MeshOrderAlteration[m, 2]


Now the mesh has coordinates at the mid side nodes and you will need data points at these positions. If you have those you can proceed with:

(* generate some random values for demonstration *)
z2 = RandomReal[{-1, 1}, Length[m2["Coordinates"]]];


Create a second order interpolating function.

if = ElementMeshInterpolation[{m2}, z2]


Then you can visualize the interpolating function and do whatever you would like to do with it.

Plot3D[if[x, y], Element[{x, y}, m2]]


To sum up, the problem with arbitrary coordinates is that no higher order mesh can be generated. If you can generate that mesh and have function values at the higher order nodes, you are good to go.