Interpolation of values on a three-dimensional surface represented by scattered points

I have got a matrix with scattered points (N x 3) representing a curved surface in three-dimensional space. Additionally there's vector (N x 1) with values (in this case these are surface temperatures) for each point from above.

My Problem: I'd like to interpolate the values on the three-dimensional surface on another set of scattered points which also represent the same surface in three-dimensional space.

I use MATLAB R2014b.

My Solutions so far:

Method 1:

Use scatteredInterpolant with 4 arguments

F = scatteredInterpolant(X, Y, Z, values)

and calculate the new values by using

values_new = F(X_new, Y_new, Z_new)

Method 2:

I've thought about the surface and came to the conclusion that for a surface in three-dimensional space the Z-coordinate is dependend on the X- and Y-coordinate, right?

Z = f(X,Y)

Therefore the dependency

value = f(X,Y,Z)

can be simplified to

value = f(X,Y)

This results in

F = scatteredInterpolant(X, Y, values)

and

values_new = F(X_new, Y_new)

Whats the mathematically correct solution for this problem?

I'really like to know which is the correct way to my interpolation problem. Maybe it's not even one of the methods I've presented.

Thanks in advance!

[Ronak Patel]

Hey Hi.. I needed some help in interpolation of data. Basically I have 3D surface coordinates (x,y,z) let say it is around 20000 and from some of this points i have a voltage value (just for 400 points) So, I want to find the voltage value for remaining points. So what should i do?

Answer 1

[John D'Errico]

Probably none of the answers you have posed is right.

You have a curved surface in R^3, represented as a set of points (x,y,z). As such, they represent a sampling from some general curved manifold, embedded in that 3-d space. In addition, you have a function, w(x,y,z), that you wish to interpolate onto another disjoint set of points that also hopefully lie on the same manifold.

What you have not said is that this manifold can be represented as a single valued function, z(x,y). For example, the surface of a sphere is a curved manifold, a general surface, but it fails this last item, so Method 2 as you outlined it may not be any kind of option. Or, if that function may in theory be single valued, but it has locations where there would effectively be a derivative singularity, so an infinite gradient. Again, z(x,y) will be useless as an idea.

Method 1 must also fail in general. In fact, method 1 will be useless ALWAYS. That scatterred interpolant will try to form a tessellation, in THREE DIMENSIONS. Since your data lies not scattered in 3-d, but on the surface of a smooth manifold, those simplexes that are generated will be terrible in terms of their value for an interpolant. That is, they will either be terribly flat, or they will use points from far away on the surface, spanning across big folds of your manifold. This is not what scatteredInterpolant was designed to solve as a problem. Sorry, but NOT.

You have a manifold. Depending on how that manifold is parameterized, a new point need not even lie on the manifold. For example, given a triangulated surface that represents APPROXIMATELY the surface of a sphere, a new point from another set need not even lie exactly in that approximately triangulated manifold. So interpolation is meaningless until you manage to map that new point onto the triangulated surface. Then a simple linear interpolation will suffice.

So far, I don't know enough about your problem. You have posed it too generally, not given sufficient information.

[bash0r]

Okay, let me describe the situation I'm in.

I measure the surface temperature of a measurement specimen using infrared thermography. The surface is three-dimensional, curved and cannot be described by any algebraic function. But each combination of X and Y has an unique Z-value (--> no undercuts). So singularities shouldn't be a problem, I guess. Let's assume the surface looks like the one on the image. The real surface does not have a pattern like the one on the image.

The surface of that specimen is represented by a point cloud which has been created by a mesh generator (I guess it was ANSYS ICEM) based on the CAD model. Therefore I have the triangulation (points and the according triangles) of that continuous surface.

After reconstructing the measured temperature field in MATLAB, I'm in possesion of the temperature for each point of the triangulation. I use this temperature field (X,Y,Z,T) as a boundary condition in ABAQUS. The temperature gets mapped to the mesh in ABAQUS (differs from the one in MATLAB) on the surface of the specimen. I use ABAQUS to calculate the heat flux (HFL) into the surface of the specimen. Unfortunately the mesh in ABAQUS is different and I have another set of points which represents the surface of the specimen.

For further processing I need the heat flux values at the same positions (X,Y,Z) like the temperature values. That's why I'd like to interpolate the heat flux data (X*,Y*,Z*, HFL).

I hope I could provide a sufficient insight into my problem. Thanks in advance!

[John D'Errico]

Good. It is easier than it might have been.

You have a problem that can be viewed as a pair of mappings, as I would describe it, of the form z(x,y) and T(x,y). The triangulation of the surface can be taken in the (x,y) plane. So at any vertex of that triangulation, you have two values, z and T. At any point (x,y) in that triangulated domain, it is simplest to do a linear interpolation within the corresponding triangle.

So, there are two options. You should be able to use tsearchn to identify a triangle that contains any point. Tsearchn also returns a set of barycentric coefficients that allow linear interpolation within that triangle, as a linear combination of the function values T at each vertex of the triangle. This is all done fairly easily, and I can show you how to do it, in a vectorized form if you need that.

If you cannot deal with that, I have written a complete toolbox that does all of this, although I tend not to give it out too often, as it takes some effort to learn to use that toolbox. My experience has been that most people then need some tutoring to use the tools, so I just avoid giving it away.

[bash0r]

Thanks for your answer! I've implemented the method you've recommended really quickly, but it seems to work (even vectorized). The execution of the tsearchn command takes a few minutes, because I work with around 1 million points. Some of the new points don't lie inside any triangle of the old points, but those are just a few (~100).

I've added the code I wrote, in case you want to have a look at it. To compare the results I used scatteredInterpolant(X, Y, value, 'linear', 'none'). Except for 77 points the results of both methods are identical. And even those 77 points have only a slight different value with no physical influence on my results.

Is it right that both methods do the same?

%%Description
% points_FEM: points representing the 3d curved surface
% values: scalar value for each point in 'points_FEM'
% points: new points (3d) for which the values have to be interpolated

%%Method 1: Mapping
% create a triangulation of the FEM coordinates
tri_FEM = delaunay(points_FEM(:,1:2));
% search for the enclosing simplexes
[T, P] = tsearchn(points_FEM(:,1:2), tri_FEM, points(:,1:2));
% create a nanmask (because some values lie outside)
nanmask = isnan(T);
% read the corresponding triangle to each point
tri_matched = tri_FEM(T(~nanmask),:);
% get the 3 values for each point of each triangle
values_matched = reshape(values(reshape(tri_matched,[],1), 1), [],3);
% use the barycentric coordinates to determine the value
values_new(~nanmask,:) = sum(values_matched .* P(~nanmask,:),2);
values_new(nanmask,:) = nan;

%%Method 2: ScatteredInterpolant
interp_value = scatteredInterpolant(points_FEM(:,1:2), values, 'linear', 'none');
values_new2 = interp_value(points(:,1:2));

[John D'Errico]

I distinctly saw you write that you HAD a triangulation already! So why would you use delaunay here? I did say to use tsearchn directly.

As for your question, yes. The scatteredInterpolant function will first do a delaunay triangulation, as the very first thing it does. So all you did was to effectively reproduce what scatteredInterpolant does.

Again, this is why you use tsearchn directly on the triangulation that you said you already had. Why do the work twice? As importantly, since a delaunay triangulation may often choose a different way to triangulate a given set of points than the one you have, it can yield subtly inconsistent predictions, since they will be based on a different triangulation than the one you were given.

For example, here is the triangulation generated by delaunay, on a regular rectangular, 20x20 mesh of points.

As you can see, the triangles chosen are rather random in their orientation. (Although it probably would make a nice random maze generator.)

[bash0r]

I'm in possession of the triangulation of the point cloud (named 'points' in the code above) which I use in MATLAB. I'd like to interpolate the values I get from ABAQUS on those points. The results from ABAQUS only contain the coordinates X,Y,Z (named 'points_FEM' in the code above) and the values. Unfortunately I dont have the triangulation for those points, that's why I used delaunay.

[bash0r]

I have managed to extract the points and the triangulation of the ABAQUS mesh. Now I have the points and the triangulation for both point clouds.

I have tried to use

 [T,P] = tsearchn(points_ABAQUS(:,1:2), tri_ABAQUS, points(:,1:2))

but half the entries in T (or P) are NaN. When I used delaunay(points_ABAQUS(:,1:2)) to generate tri_ABAQUS, tsearchn delivered a result for almost all the points.

The help of tsearchn says: " tsearchn requires a triangulation TRI of the points X obtained from delaunayn."

Do you know what the problem could be? Maybe I can't use my own triangulation?

[John D'Errico]

Sorry. I guess tsearchn requires a delaunay triagulation. Its been a while since I played with it. My own tools do not have any requirement at all in that context, but they use a more robust algorithm that can handle a totally general triangulation. A quick test shows that my code is slower though.