What algorithm does the scattered interpolant 'boundary' extrapolation use when it extrapolates to a convex hull?

The documentation within the 'boundary' extrapolation mode for the scattered interpolant states that it " [evaluates] to the value of the nearest point on the convex hull." However, if you have a non-convex domain smaller than its convex hull, then that specification seems to necessitate an initial extrapolation to approximate the function values on the convex hull that were external to the original domain. I would appreciate it if someone could provide the specific details of this algorithm for me or clarify the description if I have misunderstood something fundamental.

[Matt J]

The documentation within the 'boundary' extrapolation mode for the scattered interpolant states that it " [evaluates] to the value of the nearest point on the convex hull."

I don't see that. In my documentation it says,

"Boundary extrapolation. This method extends the values on the interpolation boundary into the extrapolation domain."

Admittedly, that is equally vague.

[Torsten]

It's in the description of the "Nearest Neighbor" extrapolation method under

https://uk.mathworks.com/help/curvefit/extrapolation-methods.html

[Atharv]

Since the official documentation is vague, I opened the .m file itself and found that quoted line as a comment on the 'boundary' option. As Torsten mentioned, it is also in the description of the nearest neighbor extrapolation method. From that documentation, I would then assume PCHIP is the default for the standard scatteredInterpolant object. However, since the documentaiton is actually for the curve fitting toolbox, this is all just speculation on my part.

[Torsten]

@Atharv

If this issue is really important for you, you should directly contact the Mathworks developers team.

Answer 1

[William Rose]

@Atharv,

I don't think you have misunderstood anything fundamental. The description is unclear. Here are examples to illustrate how it works. We will check each how extrapolation method "boundary" works with each of the three interpolation methods: "nearest", "linear", and "natural".

Example 1. Method=nearest, extrapolation method=boundary. In this case, it appears that the extrapolated value is the value of the nearest of the original points.

[X,Y]=meshgrid(0:4,0:3); 
% Next: move corner point so data set is concave at corner
X(end)=3.3; Y(end)=2.3;
V=(X.^2+Y.^2).^0.5;             % V=distance from origin
% Display X,Y points
figure; plot(X,Y,'ro','MarkerFaceColor','r')

Now examine the behavior around the top right corner of the plot above (x=3 to 5, y=2 to 4), since that is where the point set is concave.

[Xq,Yq]=meshgrid(3:.05:5,2:.05:4);  % query points
% original points, near the query points (for plotting)
XL=X([15,16,19,20]); YL=Y([15,16,19,20]); VL=V([15,16,19,20]); 

% Compute scattered interpolant, method=nearest, extrap.method=boundary
F1=scatteredInterpolant(X(:),Y(:),V(:),'nearest','boundary');
Vq1=F1(Xq,Yq);                        % interpolate and extrapolate

% Plot results
figure; surf(Xq,Yq,Vq1); colorbar; hold on; view(0,90)
plot3(XL,YL,VL,'ro','MarkerFaceColor','r')
xlabel('X'); ylabel('Y'); zlabel('Z'); title('nearest, boundary');

Careful examination of the points above suggests that the extrapolated value equals the value of the nearest original point, not the value of the nearest point on the boundary of the convex hull. For example, at Xq,Yq=4.05,3, extrapolated Vq1=4.47, which is the value of the closest original point (at X,Y=4,3). The closest point on the convex hull is at about X=3.525,Y=2.475, and this point has extrapolated value V=4.02 (by nearest neighbor). So, if Vq were using the extrapolated value at the closest point on the convex hull, Vq would be 4.02 when Xq,Yq=4.05,3. But Vq1 is not that.

Example 2: Method=linear, extrap.method=boundary.

% Compute scattered interpolant, method=linear, extrap.method=boundary
F2=scatteredInterpolant(X(:),Y(:),V(:),'linear','boundary');
Vq2=F2(Xq,Yq);                        % interpolate and extrapolate

% Plot results
figure; surf(Xq,Yq,Vq2); colorbar; hold on
plot3(XL,YL,VL,'ro','MarkerFaceColor','r')
xlabel('X'); ylabel('Y'); zlabel('Z'); title('linear, boundary');

Adjust the view to further appreciate and understand the plot.

In the case above, points along the convex hull boundary are intepolated linearly between the original points that are on the convex hull, ignoring original points interior to the convex hull. Points outside the convex hull are extrapolated from the convex hull boundary, and are not affected by points interior to the convex hull.

Example 3: Method=linear, extrap.method=boundary.

% Compute scattered interpolant, method=natural, extrap.method=boundary
F3=scatteredInterpolant(X(:),Y(:),V(:),'natural','boundary');
Vq3=F3(Xq,Yq);                        % interpolate and extrapolate

% Plot results
figure; surf(Xq,Yq,Vq3); colorbar; hold on
plot3(XL,YL,VL,'ro','MarkerFaceColor','r')
xlabel('X'); ylabel('Y'); zlabel('Z'); title('natural, boundary');

Adjust the view to understand the surface better. It looks the same as the plot when method='linear'.

This case is like the linear case: points on the convex hull boundary are intepolated linearly between the original points that are on the convex hull, ignoring original points interior to the convex hull. Points outside the convex hull are esitmated from the points on the convex hull boundary. Inside the convex hull, the esitmated values are slightly different with natural method than with the linear method, consistent with natural versus linear interpolation.

If you vary the height of the point inside the convex hull, you will see that the results with method="nearest" still work the same way. With method="linear" or "natural", the interpolated points on the convex hull boundary, and the extrapolated points outside the convex hull, are not altered by the adjustment of the nearby point inside the convex hull. Their height is determined only by the original and interpolated points along the boundary of the convex hull.

[Atharv]

Thank you for the detailed response; it definitely cleared up a few misconceptions I had about how the choice of interpolation and extrapolation methods relate. So would it be accurate to summarize the functionality as:

1. The points on the convex hull are first interpolated via the user-selected algorithm
2. The same algorithm is then given the only the points on the boundary of the convex hull, original or recently interpolated, as its input
3. The algorithm, now basing its approximations solely on boundary points, procedes to extrapolate points external to the convex hull in the same manner that it would attempt to interpolate points

?

[William Rose]

@Atharv, I agree with your points 1 and 2.

I would not state point 3 exactly the way you did, because "extrapolate ... in the same manner that it would attempt to interpolate points" is impossible, for the linear method, when there are no point or points beyond the query point. One might expect that method linear would extrapolate using the slope that it had leading up to the edge. But it doesn't. It just flattens out. Methods nearest and natural also flatten out in the extrapolation zone.

The three examples below use the same original points as above, but here the query points extend all around and outside of the original points. For each method, the default 3D view and an overhead view are shown.

[X,Y]=meshgrid(0:4,0:3); 
% Next: move corner point so data set is concave at corner
X(end)=3.3; Y(end)=2.3;
V=(X.^2+Y.^2).^0.5;             % V=distance from origin

% Extended example going all the way around the point set.
[Xq,Yq]=meshgrid(-2:.05:6,-2:.05:5);  % query points
% Compute scattered interpolant, method=nearest, extrap.method=boundary
F1=scatteredInterpolant(X(:),Y(:),V(:),'nearest','boundary');
Vq1=F1(Xq,Yq);                        % interpolate and extrapolate
% Plot results
figure; 
h_sp1=subplot(121);
surf(Xq,Yq,Vq1,'EdgeColor','none')
hold on
plot3(X,Y,V,'ro','MarkerFaceColor','r')
xlabel('X'); ylabel('Y'); zlabel('Z'); title('nearest, boundary');
h_sp2=subplot(122);
copyobj(get(h_sp1, 'Children'), h_sp2);
xlabel('X'); ylabel('Y'); zlabel('Z'); title('nearest, boundary');
colorbar; view(0,90)

% Compute scattered interpolant, method=linear, extrap.method=boundary
F2=scatteredInterpolant(X(:),Y(:),V(:),'linear','boundary');
Vq2=F2(Xq,Yq);                        % interpolate and extrapolate
% Plot results
figure; 
h_sp1=subplot(121);
surf(Xq,Yq,Vq2,'EdgeColor','none')
hold on
plot3(X,Y,V,'ro','MarkerFaceColor','r')
xlabel('X'); ylabel('Y'); zlabel('Z'); title('linear, boundary');
h_sp2=subplot(122);
copyobj(get(h_sp1, 'Children'), h_sp2);
xlabel('X'); ylabel('Y'); zlabel('Z'); title('linear, boundary');
colorbar; view(0,90)

●

% Compute scattered interpolant, method=natural, extrap.method=boundary
F2=scatteredInterpolant(X(:),Y(:),V(:),'natural','boundary');
Vq2=F2(Xq,Yq);                        % interpolate and extrapolate
% Plot results
figure; 
h_sp1=subplot(121);
surf(Xq,Yq,Vq2,'EdgeColor','none')
hold on
plot3(X,Y,V,'ro','MarkerFaceColor','r')
xlabel('X'); ylabel('Y'); zlabel('Z'); title('natural, boundary');
h_sp2=subplot(122);
copyobj(get(h_sp1, 'Children'), h_sp2);
xlabel('X'); ylabel('Y'); zlabel('Z'); title('natural, boundary');
colorbar; view(0,90)

[Atharv]

My guess for the exact method of linear extrapolation from the boundary is something like the following:

Given a set of points on the boundary, select three sequential points p1, p2, p3. Compute (p3-p1)/||p3-p1|| as the unit tangent vector, and then apply the appropriate rotation to approximate the unit normal vector. Taking a central difference along this direction gives a fair approximation for the normal derivative at p2. This would then explain the "flattening out" we see in the diagrams, as most of that gradient is "concentrated" along the tangent direction due to the sharp turns and mostly linear nature of the convex hull of this specific shape.

This hypothesis also appears to superficially match the figures shown in the documentation page: https://www.mathworks.com/help/matlab/ref/scatteredinterpolant.html

Once again, this is all just speculation, and I appreciate the detail in your examples. I am going to mark this thread resolved, as your first answer made the mechanisms of the scatteredInterpolant more clear even if we are still uncertain about the exact details, which likely only Mathworks staff know for certain. Due to the nature of my use case, I now plan to use a more transparent algorithm or implement my own. Thank you for your responses.

[William Rose]

@Atharv, you're welcome. Good luck with your work.