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

18 views (last 30 days)
MG
MG on 25 Mar 2021
Commented: MG on 29 Apr 2021
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 27 Apr 2021
Edited: MG on 27 Apr 2021
Hi Bruno.
Thanks for your information and for the reference !!
I didnt directly find what I'm looking for in CGAL. The closest I found, which is maybe what you had in mind, is https://doc.cgal.org/latest/Surface_mesher/index.html#Chapter_3D_Surface_Mesh_Generation, however I'm not sure how to use it in my context. CGAL looks like a diverese and powerful C++ tool, so all needed tecniques might be there somewhere :)
In my case, I "only" want to have a code that adaptiively find points (x_i, y_i); such that when plotting z(x_i,y_i) and linearly interrpolating between the (x_i, y_i) points this will represent well (by some reasonable criteria) the full smooth function z(x,y). Very much like how surf actually does it (but matlab's fsurf is just a black box with no documented settings to adjust the the point sampling precission in a controllable way).
I've kept serching and found somethinig in Python, which miight be interresting to adapt to matlab: https://adaptive.readthedocs.io/en/latest/docs.html
It's a pitty that matlab is not direclty poviding tools for this -- as you pointed out, it seems almot inexistent in matlab.
MG
MG on 29 Apr 2021
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
KSSV
KSSV on 29 Mar 2021
In the above code linspace os ised to generate equal spacing array, you can try replacing that with other options like logspace etc.
MG
MG on 29 Mar 2021
Edited: MG on 30 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
MG
MG on 6 Apr 2021
Edited: MG on 12 Apr 2021
I fully agree, and this is why I tried to emphasis that "I understand that there is no perfect way (and we have to assume that the [functions] are fairly smooth and reasonably well behaved)".
Still Matlab provides both local and global minimizers, as well as 2D integration routines and dynamical 2D plotting functions. These are all "very challenging" tasks, which can in general not be solved for every function ---as you point out. Hence, such matlab routines can only "do their best" and allow the user to adjust some settings for their problem in hand.
Therefore, I initially thought that matlab might provide also routines (which do their best) to adaptively map out 2D-surfaces. Hence, I asked if some routines like that exists?
Matlab's fsurf sometimes do the job well enough, maybe it can do even better if I knew in detail what it does and the user could adjust its settings -- but I could not find any proper documentation about this (and using fsurf already seems a bit like a dirty work-around approach to me). I could also only trace fsurf's meshing routine to be hidden inside matlab.graphics.function.FunctionSurface, which is an inaccessible binary machine code.
Similarly, I could imagine to try to use Matlab's 'integral2, if I could get access to all its coordinate points used when performing its 2D-integral -- but not immediately clear to me how to do this ?
Regarding your "ponder" problem: In this case the function doesn't have a well defined limit when x and/or y approach 0. Hence, it falls outside my definition of "fairly smooth and reasonably well behaved". However, for say, x,y>0.05 this function is already challenging (and maybe not yet "unreasonable"). For example with a code like below (using fsurf as a "black-box") a ~10% precession is reached (in a few seconds), but then it's difficult to push the precision to become much better with fsurf. So here another method, or other internal fsurf settings, would be needed ---> this is why I put my question into this forum to see what options Matlab provide to adaptively map 2D-surface functions.
Thanks!
x_min = 0.05;
y_min = 0.05;
z_exact = @(x,y) sin(1./x) + cos(1./y);
subplot(1,3,1) % map function with fsurf
z_surf = fsurf(z_exact,[x_min 1 y_min 1],'MeshDensity',40,'Visible','on','EdgeColor','none');
z_surf.AdaptiveMeshDensity = 10;
%
colorbar
view(2)
set(gca, 'XScale', 'log','YScale', 'log');
xlabel('x'),ylabel('y'),zlabel('z'),title('fsurf')
subplot(1,3,2) % plot sample grid
plot(z_surf.XData',z_surf.YData','.k','MarkerSize',1)
xlabel('x'),ylabel('y'),title('Mesh points')
set(gca, 'XScale', 'log','YScale', 'log');
z_scatter = scatteredInterpolant(z_surf.XData',z_surf.YData',z_surf.ZData','natural','nearest');
f_min = @(x) -abs( z_exact(x(1,:),x(2,:)) - z_scatter(x(1,:),x(2,:)));
subplot(1,3,3) % Plot deviation at random points (in log-scale)
x = x_min*exp(log(1/x_min)*rand(1,5e4));
y = y_min*exp(log(1/y_min)*rand(1,5e4));
z = -f_min([x;y]);
hold on
zscaled = z*200;
cn = ceil(max(zscaled));
cm = colormap(jet(cn));
scatter3(x, y, z, [], cm(ceil(zscaled),:), 'filled')
cb = colorbar;
cb.Label.String = 'relative deviation';
caxis([0 max(z)])
set(gca, 'XScale', 'log','YScale', 'log');
view(2)
[a b] = max(-f_min([x;y]));
diff_max = -f_min([x(b);y(b)])
xlabel('x'),ylabel('y'),title(['Absolute deviations, max = ' num2str(diff_max,2)])
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.

Categories

Find more on Interpolation in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!