Why is ScatteredInterpolant (linear) returning a value far outside of input data?

So I'm using ScatteredInterpolant to interpolate the values of a measurement on the surface of a helicopter blade. I'm showing it graphically in the figure. The input data is shown with the faint lines at a 45 degree angle. The output is the coloring of the surface. Right in the center is a bright spot. The highest input data value is ~119. The output value at that location is about 2512.

I've attached the data which generated that figure. The function used to generate the Interpolation is:

Interpolant = scatteredInterpolant(inX, inY, inZ, inV,'linear','nearest');
outV = Interpolant(outX, outY, outZ);

For reference, the bad output data point is the 2713 entry in outV. How can I possibly get an output value that is more than an order of magnitude bigger than any input value, when using linear interpolation? Is this a bug?

A quick way to visualize the error:

figure;
plot3(inX, inY, inV, 'o');
hold on
plot3(outX, outY, outV, 'o');

[Jason Whitfield]

What version of MATLAB are you running? I ran the code in R2018a, and the maximum value in outV was 118.97.

[Dan K.]

2018a

[Jason Whitfield]

I've attached the results I get by running the following code.

clear;
load('failingData.mat');
outV = Interpolant(outX, outY, outZ);

The outV data looks good to me. I'm not sure why you're getting a different result.

[Dan K.]

Yes, I see that you get good looking results, but when I execute the exact same code, I get the values of outV you see in the plot above. I've submitted a bug report to TMW, and I'll post an update if I find out more.

Answer 1

[John D'Errico]

You are making a major error here in what you are trying to do.

plot3(inX,inY,inZ,'bo')
box on
grid on

Now, rotated one way, we see this surface.

Now, rotated slightly differently, we see this:

So you have essentially a nice surface in 2-dimensions. But there is NO information about what happens when you move off that surface. Now, add in the points you are trying to "interpolate".

hold on
plot3(outX,outY,outZ,'r+')

Yep. You are trying to predict values that lie HUGELY away from the 2-manifold of data that you have. Massive extrapolation here.

So, regardless of what scatteredInterpolant did, your results will be complete, meaningless garbage. You don't have three independent variables, comprising a true point cloud. You have a 2-manifold of data.

So what will scatteredInterpolant try to do here? Linear interpolation for any point will involve a tessellation of that domain. It will not be some sort of triangulation of the surface, but a set of tetrahedra spanning, dissecting that domain. At each of those points, you have a value of V that you wish to interpolate.

So scatteredInterpolant uses a tessellation of that domain. Thinking that you have a real point cloud.

tess = delaunayn([inX,inY,inZ]);

whos tess
Name          Size              Bytes  Class     Attributes
tess      43112x4             1379584  double

So there are 43112 tetrahedra. Its kind of hard to visualize, since there are so many edges. But the volume looks like this:

Most of those tetrahedra are HUGE things. Spanning across that void where no data exists. scatteredInterpolant has no clue what you did. At any poiint in that domain, it finds the tretrahedron that contains your point. Then it does LINEAR interpolation across that tetrahedron. But remember that most of the tetrahedra are huge things, connecting one "wing" of that surface with the other. (No pun intended with use of the word wing here, on what is essentially airfoil data.)

ScatteredInterpolant just does what it is told, having no idea that when you try to interpolate some point in that volume, it is creating meaningless gibberish as a result. I'm sorry, but you simply cannot use scatteredInterpolant to produce a meaningful result from this data, as you are trying to do.

Do I really need to figure out why you got something mildly strange? Garbage in, garbage out. I guess I could go a bit more deeply into what you have, if you really want it. But it would be a pure exercise in trying to debug code that I don't have, because as I recall, scatteredInterpolant is supplied to us in compiled form.

[Dan K.]

I appreciate the feedback, but I have to say there are a few reasons why I think you are mistaken.

First: the data point that comes out with a strange value is dead center in the middle of the surface, well surrounded by adjacent data points. In my actual application, I zero out the "interpolated" results of any point that lies outside the surrounding input data. Only the points which exist within the rectangular prism defined by the input coordinates are kept.

Second: No matter how you partition the space, there is no circumstance in which a linear interpolation should result in a value that exceeds the maximum of the input data. Even thinking about extrapolation wouldn't justify it, since the extrapolation method is nearest.

One other item which I've discovered since I did my original posting: this effect doesn't occur when the interpolation method is set to 'natural'.

Oh, and yes, you are correct that scatteredInterpolant is a built-in function.

[John D'Errico]

Sorry. But no. You cannot use that kind of data to interpolate.

No. It is NOT dead in the middle of the surface. It seems that way. Yes, it is close. But your surface is a 2-dimensional thing, that will be approximated by a piecewise linear interpolant, based on a tessellation that will be HIGHLY UNSTABLE, due to the fact that as I showed in the very first plot, the surface is a TWO DIMENSIONAL MANIFOLD.

No matter what, the point will be slightly off the surface, even if only by a tiny amount. So what happens is you will have a HUGE tetrahedron, spanning long distances to interpolate. The result will be garbage. How exactly they implemented the extrapolation is difficult to know, since we cannot see the code, but extrapolation in multiple dimensions is a non-trivial thing to do.

You can believe me or not. (Given that I could have written scatteredInterpolant myself, and that I have written similar code in the past, I do know what I am talking about.) Regardless, you did something silly in trying to use that code as you did. You got garbage. I'm just not even remotely surprised. But I cannot track down exactly what it did, because I cannot trace down into the code. As I said, scatteredInterpolant is compiled. I was not guessing about that fact.

Let me repeat, using scatteredInterpolant to try to interpolate the data that you have will result in garbage. That you got something strange in the form of the garbage you got will be impossible to track down. You can send this in as a bug report if it really bothers you, because the only people who can track it down are those who wrote the code. I would not waste a lot of time on it, because this tool is not designed to interpolate the kind of data you have. You cannot use a 3-dimensional interpolation tool to interpolate what is fundamentally 2-dimensional data.

Now, were I to seriously try to use your data for interpolatory purposes, in a way that would not produce random garbage? I would ...

1. Formulate the surface you have as a triangulated surface, so composed of TRIANGLES, not tetrahedra, having volume.

2. For any point that you wish to interpolate, find the closest point to that surface, and the location of that closest point on the surface, thus as a minimum distance projection onto a triangulated surface. This will result in knowing which triangle that closest point lies in, and exactly where, probably defined by a set of barycentric coordinates.

3. Then do a linear interpolation over that triangle, using those same barycentric coordinates.

All of these operations can be done in a well-posed, robust way, although step 1 may require some effort on your part to figure out. The result will be stable. It will not be subject to random garbage. It will be a better use of your time that to try to get scatteredInterpolant to work here, which is just a bad idea.

[Dan K.]

John, I don't want to give you an incorrect impression. I am well aware of your expertise, and I'm not trying to diminish it at all. In fact I have used your FEX submissions many times over the years. What is bothering me (and I have submitted a bug report on this) is that no matter how poorly I've defined the problem, any form of interpolation, using any four data points that are input, should never give a value which is an order of magnitude greater than any of the input data. That would hypothetically be possible if it were extrapolating the data, but the extrapolation method I specified was 'nearest', so even that shouldn't give a value outside the bounds of the input data. Given how close my problem output is to the center of the input data, I'm wondering if there is a rescaling that is being done in the compiled code that is resulting in a divide-by-near-zero.

Anyhow, I will look at your suggestions about re-mapping the data by unwrapping the manifold.

And again, thank you.

[John D'Errico]

Sorry.

There are two issues here.

First is that of the extrapolation not performing as expected. You could report that as a bug. There are many ways of doing extrapolation in multiple dimensions. I could guess what they are doing, but I might be wrong. And without the code, I cannot do any debugging. And on some of the codes I have written, they TRY not to extrapolate. But they can go outside the bounds slightly. In the case of a tessellation based code, if there are some highly degenerate tetrahedra created from the Delaunay tessellation, then numerical problems might arise, and that might be as good as it can do. Its just wild guessing without the code to see. So you may choose to send this in as a bug report.

The second issue is how to do interpolation on a problem like this. To make it easier, consider this simple set of points:

t = linspace(0,pi/2,10);
x = cos(t);y = sin(t);
plot(x,y,'o')
hold on
plot(sqrt(2)/2,sqrt(2)/2,'rx')

Now, you and I know that the point marked with a red 'x' is on the curve.

Suppose we have some value of z assigned to every point on that curve? We might even hope to use interpolation to interpolate a value for that point. But can we use a scatteredInterpolant like approach?

tri = delaunay(x,y);
trimesh(tri,x,y)

In fact, the x lies outside of the convex hull of the points, because the convex hull is a piecewise linear thing, approximating here an arc of a circle. So the convex hull lies always inside the circle, except for the points where it touches the perimeter. (That presumes no noise at all.)

But any point that we try to interpolate will often be inside some highly degenerate simplexes. scatteredInterpolant does not work well on this class of problem.

So my suggestion is as I said above. Formulate the 2-manifold as a triangulation floating in the 3-d domain space. You may be able to do that using an algorithm like CRUST. I think it can be made to work in 3-d, although I am not an expert on the CRUST algorithm.

A quick search suggests that is possible, based on the keywords "CRUST algorithm 3-d".

Once that is done in whatever way you can, the rest is simple. Project any point onto the nearest triangle, as a minimum distance projection. Once you have the projected point, interpolation is the easiest part of the process.

[Dan K.]

Again, thank you. The good news is that I already have the manifold geometry. That's the patch object I showed in the first graphic of my original post. The Out X,Y,Z are just the vertices of a downsampled version of the patch. What I'm trying to figure out now is how to do a coordinate transform to convince Matlab that it really is a 2D triangulation problem.

[John D'Errico]

You can't do that, at least not easily tell MATLAB that an arbitrary triangulated 2-manifold could be viewed as some sort of regularly gridded surface.

If, however, you have a triangulated surface, then the answer is easy enough. Some years ago, I wrote closestmanifoldpoint.m.

Given a triangulation and a list of points in 3-dimensions, it finds the simplex which lies closest to each point. It returns a set of barycentric coordinates inside that simplex, used for interpolation. But if you supply information about any mapping, defined at each node of the triangulation, it also does the interpolation for you.

It is part of a larger set of tools, but simply used without them.

You will need to define a struct to pass in the information about the triangulation. The second argument, sc must be a struct with three fields, sc.domain, sc.tessellation, and sc.range.

The domain field will be a list of points in R^3, one point per row, so the list of vertices of the triangulation. So sc.domain must be an NPx3 array, assuming NP nodes in the triangulation.

sc.tessellation will be an NTx3 array, listing the connectivity for each triangle in the surface, so relative references to the rows of the domain.

Finally, this code can also do interpolation as I said. That is accomplished by having a range field in the struct. The range field must be of size NPxk, where each column of this array is a variable to be interpolated. Assuming you are just using this to predict one value at each point, the range field will be a vector, of size NPx1.

Call it like this:

[x,dx,whichsimplex,bcc,rangeInterp] = closestmanifoldpoint(X,sc);

where X is an array of points to be mapped to the surface, and then interpolated.

So all you need to do is supply a triangulation. I've attached closestmanifoldpoint, as well as the only subordinate function that I noticed, rowspaces.

For example, here is a simplicial complex struct that defines a triangulation of the surface of a unit sphere.

sc
sc = 
  struct with fields:

            domain: [62×3 double]
      tessellation: [120×3 double]
             range: [62×1 double]

The color coding in that figure comes from the range field, which I generated as a uniform random variable, rand(62,1). Not very creative, I know.

[x,dx,whichsimplex,bcc,rangeInterp] = closestmanifoldpoint(rand(10,3),sc);

rangeInterp
rangeInterp =
       0.592628927726262
       0.329792850128984
       0.368604533985401
       0.633207288915207
       0.325491526196841
       0.278335770755433
         0.3232807215718
       0.469008085876214
       0.409532773231336
       0.656797929332787

So, given a set of random points in the unit cube, it found the closest point on that surface, then did an interpolation. Yeah, its random garbage here, but only because of the intensely creative way I created the vector to fill in that range field.

That means your problem just got very easy to solve, as long as you have the triangulation defining that surface manifold. And, because I wrote the code, I do know what it does.

[Dan K.]

Wow... John, thank you so much for your help. I'm still working on the details of how to implement it, since the triangulation I have is for the query points, not the data points (as you've shown in your example). But thank you.

Dan

Following up... I had an idea that I think should work, but I figured I would run it past you to see if you think I'm nuts.

What I'm ultimately trying to do is to interpolate the in-plane strain measurement on the surface of my blade. So what I'm thinking is trying to do this:

Since my blade is mostly just a uniform extrusion of a single cross section (i.e. the manifold) I should be able to map any point on the blade (except for the tapered ends) onto a 2 dimensional plane by "unrolling it". What I'm envisioning here is if I could cut the blade along its trailing edge and iron it flat, then I would have a single sheet that has all the measurements in plane. Then I do my interpolation, in-plane, using the 2D version of scattered Interpolant, and then remap it back onto the 3D geometry of the original blade.

So if I start with the original blade geometry showing all my data points.

I can do the mapping by keeping the Y value the same (i.e. distance down the blade), and for each data point I determine the equivalent XX location which is the distance around the blade from the trailing edge, using something like your arclength function. I know the cross section of the blade, so I can snap each data point's X-Z value onto the cross section and integrate distance around the perimeter to create my data map in the new XX-Y coordinate system.

I do the same thing with my query points, keeping track of each query point's mapping from X-Y-Z into XX-Y. Then do the interpolation in the XX-Y domain, and assign the value back to that query point's original X-Y-Z coordinates.

Naturally, I can't do this at the ends of the blade where my XX-Y plane projection would be distorted, but I don't actually have any data there.

Does this make sense?

Thanks, Dan