# Is there a way to adaptively sample 2-D surfaces ?

49 views (last 30 days)
MG on 25 Mar 2021
Commented: MG on 29 Apr 2021 at 15:46
Hi,
I have 2D-surface functions z(x,y) that I want to sample efficiently. Does matlab include any such addaptive mesh-refinement functions and, if so, how to use it?
Backgrund: I want to find sample points on my functions, so that I later can intepolate (using e.g. linear 'scatteredInterpolant') between those sampled points and by this effectively recover my "full" surface function. I will probably need around 1e6 sample points on each surface -- in order to precisely enough capture the full shape of my 2D surfaces z(x,y) -- but still these points must be mainly concentrated in those regions where z have peaks or edges.
Let me give a simple example (my functions typically also have narrow edges, and higher peaks):
z = @(x,y) 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2)- 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) - 1/3*exp(-(x+1).^2 - y.^2);
Now, I want to adaptaively sample x, y points such that it allows to capture the full behaviour of z.
The only function I could find that does something like this was 'fsurf':
h = fsurf(z, [-5 5 -5 5],'MeshDensity',10,'Visible','off');
I can then access fsurf's sampled points via:
x=h.XData; y=h.YData; z = h.ZData;
However, I want to control the adaptive refinments. fsurf gives no such options (excpet the initial MeshDensity). I accidentally noticed that I can also set a parameter called h.AdaptiveMeshDensity to larger values (e.g. to 9) after fsurf is executed (if I keep fplot's 'empty' window open!). This setting seems to allow to futher refine the sampling in the most relevant regions. Let me illustrate it by this code:
figure()
plot3(h.XData,h.YData,h.ZData,'.b')
hold on
plot3(h.XData,h.YData,h.ZData,'.r','MarkerSize',10) So I do obtain sampling points adaptively in this way, but:
1. I can not find any documentation about the setting`AdaptiveMeshDensity', nor is this atribute listed when I run get(h), but still this setting influences the adaptive sampling in fsurf?
2. I don't want fsurf to open a figure, but I don't understand how to prevent it, if using fsurf?
3. I imagine that there exists a more fundamental matlab function/operation (which fsurf might even be using) to adaptively find good mesh/sampling points. However, I could not find any such adaptive-mesh-refinement routine(s) within matlab?
Thanks for any useful help and inputs!
(p.s. I noticed there is a function generateMesh, but unclear if this can be of any use here)
MG on 29 Apr 2021 at 15:46
Thanks for the reply Bruno. As far as I understand that code you used does the process of reducing the number of faces in a surface mesh while still keeping the overall shape. However, we need the "reverse" process -- to start with a limited number of points on a function surface and then infromatively add more and more points until the functional surface is recovered precise enough.

KSSV on 29 Mar 2021
z = @(x,y) 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2)- 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) - 1/3*exp(-(x+1).^2 - y.^2);
m = 100 ; n = 100 ; % change these values
xi = linspace(-5,5,m) ;
yi = linspace(-5,5,n) ;
[x,y] = meshgrid(xi,yi) ;
z = z(x,y) ;
surf(x,y,z)
MG on 29 Mar 2021
Manually you can of course choose any grid you want, but that is not an automatic adaptive approach. I'm also not interested in a Cartesian rectangular grid (as meshgrid provide).

Bjorn Gustavsson on 5 Apr 2021
Edited: Bjorn Gustavsson on 5 Apr 2021
It is not entirely clear what your "problem situation" or "design desires" are. These might be very important for what type of solution to suggest. There is a sparse-grid-interpolation tool on the file exchange that is rather powerful and clever. However, you indicate that you have a general 2-D function you want to approximate efficiently on a completely unstructured grid - which sounds like quite a challenge. Once upon a time a coleague and I tried to make a "sparse image representer" by iteratively adding points to a Delaunay-triangulation where the absolute difference were largest:
I = double(I);
X = [1 1 512 512];
Y = [1 512 512 1];
F = scatteredInterpolant(X(:),Y(:),(diag(I(Y,X))));
Iapp = F(x,y);
[Imax,i1,i2] = max2D(abs(double(I)-Iapp)); % My own max2D-version
IMax = Imax;
subplot(1,3,1)
imagesc(abs(double(I))),colorbar
while Imax > 10
subplot(1,3,2)
imagesc(abs(double(Iapp))),colorbar
subplot(1,3,3)
imagesc(abs(double(I)-Iapp)),colorbar,hold on,plot(X,Y,'r.'),hold off,
X = [X,i2];
Y = [Y,i1];
F = scatteredInterpolant(X(:),Y(:),(diag(I(Y,X))),'natural');
% 'natural' in this case gives vibes of Turner, 'linear' reminds me of cubist-period Picasso
Iapp = F(x,y);
title(numel(X))
[Imax,i1,i2] = max2D(abs(double(I)-Iapp));
IMax(end+1) = Imax;
drawnow
end
This was mainly for laughs and giggles. One could also turn to a quad-tree-type parameterization (the attached function uses a general polynomial model-fiting-function orignating from John d'Ericco similar to: polyfitn).
Additional tools that might be of interest on the file exchange might be: griddatalsc, gaussian-interpolation-with-successive-corrections and barnes-interpolation-barnes-objective-analysis.
In the end it all comes down to how you plan to calculate the fit between your function and approximant. Is the function very expensive to calculate? Can you live with a regular grid or fixed number of test-points for evaluating the fit quality, are those points to be adapted? For more detailed advice more details about the problem is necessary.
HTH
Bjorn Gustavsson on 7 Apr 2021
If you'd be happy to get the grid from the quadrature-functions then it sounds like you could just have a look at the source-code of quadgk and integral2 and see what you should extract from those functions. They are normal .m-files at least in 2020a.
HTH