Finding the ranges where scatteredInterpolant is above threshold

I have been using scatteredInterpolant to reconstruct surfaces from actual measured data points. To address the slight misplacement of the object among different measurements, a reference marker was placed on the object to enable correction. Ultimately it should look for multiple points but I would like to start from the maxima.

To put it mathmatically, for function f and f', I need to know (x1,y1) and (x2,y2) where f and f' reaches its maximum.

I have tried to create a fine enough grid but it appears to reach the resolution we need it exceed the memory limit just by creating the array

ThickRBInterpolant=scatteredInterpolant(thetaThickB,zThickB,ThickRB);
>> thetalinC=linspace(0,180,180000);
>> zlinC=linspace(0,20,20000);
>> [THETAC,ZC]=meshgrid(thetalinC,zlinC);
Error using repmat
Requested 20000x180000 (26.8GB) array exceeds maximum array size preference. Creation of arrays greater than this limit may
take a long time and cause MATLAB to become unresponsive. See array size limit or preference panel for more information.

I know I can write a loop to reduce the range (and increase resolution) in steps, or just run it on a supercomputer, but am just wondering if there is already tools available?

Edit:Walter Roberson made a good comment on scatteredInterpolant not possible to exceed max sampled points, so I guess maxima is not the way to go. The question is thus changed to:

For a scatteredInterpolant f, find all (x,y) where f is larger than threshold a.

In theroy this is easy:

Thickness=ThickRBInterpolant(THETAC,ZC);
MarkerTHETA=THETAC;
MarkerZ=ZC;
MarkerTHETA(Thickness<a)=NaN;
MarkerZ(Thickness<a)=NaN;

But this would still run into the same "out of memory" wall. I am confident that the area where Thickness>a is continuous, is there anyway to find the boundary w/o using too much RAM?

[Walter Roberson]

For the revised question: is there only one such area, or could there be disjoint areas?

[Yi-xiao Liu]

There should be only one such area

Answer 1

[Walter Roberson]

scatteredInterpolant uses linear interpolation by default. The maximum possible value when you interpolate in the interior, is equal to the maximum input.

scatteredInterpolant uses linear extrapolation by default. If it is possible in your situation that the function as sampled is not uniformly surrounded by constant values smaller than the next closest interior points, then there could be an angle where the exterior points could be considered to have an upwards slope. In such a case, with linear interpolation, sampling sufficiently far away from the data boundary will give you indefinitely large values.

If you 100% confident that the maximum is somewhere on the interior, then do not bother to do the interpolation to find the maxima: just take the maxima of the input data.

Otherwise, if there is any possibility that the maxima might be along the outside, then do interpolation -- but you only need to do sampling along your outside boundary, not on the interior (the interior values cannot exceed the maxima of the input data.)

The logic behind all of this would be different if you were using a non-default interior interpolation, such as if you were using cubic spline.

If you do want to use cubic spline, then John D'Errico is much more qualified than I am to talk about possible optimization of finding maxima for cubic splines.

[Yi-xiao Liu]

I guess you are right, I will edit my question

[Walter Roberson]

Triangulate first:

https://www.mathworks.com/help/matlab/ref/delaunaytriangulation.html

Then use triangulate() https://www.mathworks.com/help/matlab/ref/triangulation.html to create a representation.

find all nodes (original points) that are above the threshold, and ask for their neighbors(): https://www.mathworks.com/help/matlab/ref/triangulation.neighbors.html

Now, all points that are above the threshold will fall strictly within the areas included in the triangles of the points above the threshold together with their neighbours. Any point that is above the threshold adjacent to a point below the threshold, you can do a simple calculation to figure out where along that line the value would equal the threshold. You do encounter complications in figuring out where to put the boundary along the face of the individual triangles... it would not surprise me if the proper line is a curve.

To save complications:

• for each isolated region as determined by points above threshold and their immediately connected neighbours:
• bounding box the region
• discretize the coordinates to your sampling grid coordinates
• generate the relevant subgrid
• interpolate at locations
• determine which ones are above the threshold

Techniques such as finding the distance along the line to where the boundary lies could potentially reduce the size of the box to search, and discretizing on a per-triangle basis instead of overall bounding box could potentially save a lot of extra calculations if the boundary is stretched out (e.g., crescent shaped)

Perhaps find the potential points within each triangle, and toss them in a bag, and unique() 'rows' them afterwards, to reduce the number of extrapolations to do.

[Yi-xiao Liu]

I am not sure I can follow you, bear with me as someone know nothing about triangulation.

ThickRBDT=delaunayTriangulation(thetaThickB,zThickB,ThickRB);
ThickRBTR=triangulation(ThickRBDT.ConnectivityList,thetaThickB,zThickB,ThickRB);

Now what's next? I assume by original points you mean measured datapoints, so the part that are above threshold is then

MarkerThetaB=thetaThickB(ThickRB>a);
MarkerZB=zThickB(ThickRB>a);

But these are just points, and the neighbours function ask fot IDs as input. How do I find the triangles whose (at least one of) nodes are among those points? (BTW it seems neighbours accept both triangulation and delaunayTriangulation as input, what's the purpose of the 2nd step?)

(As far as I understand) what I now have is a pile of triangles that are above the threshold as well as those that the boundary runs through. That's nice though still not the actual boundary. I assume you could reduce the triangles to at least one node above threshold and at least one below (though I have no idea how to do it). After that have the triangles to intercept with plane Z=a and put all the intercepts together to get the boundary (still no idea how to do it)?

PS: I don't know if it would make the problem easier by asking for the center of the boundary (point where total distance to all boundary points is minimized).