MATLAB Answers

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

49 views (last 30 days)
MG
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()
h.AdaptiveMeshDensity=9;
plot3(h.XData,h.YData,h.ZData,'.b')
hold on
h.AdaptiveMeshDensity=2;
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)
  5 Comments
MG
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.

Sign in to comment.

Answers (2)

KSSV
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)
  3 Comments
MG
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).

Sign in to comment.


Bjorn Gustavsson
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 = imread('liftingbody.png');
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
  4 Comments
Bjorn Gustavsson
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

Sign in to comment.

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by