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

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)

[Walter Roberson]

AdaptiveMeshDensity is undocumented; it was mentioned in https://www.mathworks.com/matlabcentral/answers/385717-how-can-i-vectorize-a-function#comment_542195 and https://www.mathworks.com/matlabcentral/answers/582839-surface-plot-with-white-patches-on-using-fsurf-but-complete-surface-using-ezsurf

[MG]

Thanks for the input Walter, these were the places where I also first saw a mentioning of the AdaptiveMeshDensity setting.

[MG]

This will not address the main question (e.g. point 3), but by adding

fig=figure('Visible', 'off')

before calling fsurf, you at least prevent fsurf's from showing its created figure (and axes).

[MG]

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]

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.

Answer 1

[KSSV]

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]

Hi KSSV,

You reply is not addressing the question. In fact, I don't even want to create a plot, which seems to be the only thing your code does.

So, your code does not automatically and adaptively adjust its sampling points depending on the functional from of the surface z(x,y). Your code is only a simple linear grid, decided manually by the the values set on m and n. Sorry, if my question was not clear from the beginning; that I'm looking for an automatic, adaptive method to sample functions. In other words, I a higer density of points should be found around potential "rapid changes" of a function and less dense in other regions -- something along what fsurf seems to do.

[KSSV]

In the above code linspace os ised to generate equal spacing array, you can try replacing that with other options like logspace etc.

[MG]

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).

Answer 2

[Bjorn Gustavsson]

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

[MG]

Thanks for your input Bjorn !

1.

I could not run your code due to: a) x and y in Iapp = F(x,y) are not defined; b) I don't have your max2D function. I therefore adjusted your code to inlclude:

[y,x] = ndgrid(1:512,1:512); % <-- Added

and

function [Mmax i1 i2] = max2D(M) % <-- Added
[M1   iM ] = max(M);
[Mmax i2 ] = max(M1);
i1 = iM(i2);

2.

The way I then understand your code is: i) You construct an interpolation function F from data points X,Y; ii) you test the goodness of that F against the true function at a set of 'reference' points x,y; iii) You then add to your set of (X,Y) the single point that gave you the maximal deviation (defined via max2D); and iv) you loop back to point i) and repeat until the maximal deviation is below a threashold value (set to 10 in your code).

It was very nice to see your code run, but the problem for me is that:

a) it's not a dynamical search of (x,y) coordinates, but it's rather a dynamical selection of pixels among an already pre-defined set of finite points (in your case the 512x512 pixels). With my continuous functions z(x,y), I do not know in advnace which x,y coordinates are optimal to test, but I need to start with some set if initial guess and then refine the set of selected x,y points (as fsurf).

b) the code seems slow, in each iteration you recalculate Iapp = F(x,y)for a large number of points (if I understood it correctly x,y are all the 512x512 pixels). The itterations can be made faster by changing from 'natural' to 'linear' interpolatiion and by commenting out all plotting commands --- but it's still much slower than matlab's fsurf. My aim is to preferable do everything within ~1 sec.

3.

Let me try to also answer your questions:

• Is the function very expensive to calculate? Moderately. The function z(x,y) takes 0.0003 seconds on my laptop for single x,y point, but it is to be evaluate about 1e10 times.
• Can you live with a regular grid or fixed number of test-points for evaluating the fit quality, are those points to be adapted? No, I can not use a regular grid, nor a pre-defined set of points to evaluate the fit quality. Intesad, good points (i.e. continuous coordinates) needs to be found adaptively.
• I'm not sure how I can best provide more needed information. Maybe I can make an analogy (but don't take it too literally -- my problem has nothing to do with picturers/images and my z(x,y) is just an analytical mathematical function!): Imagine you have a picture with infinite pixel resolution and with large range in continuous possible intensities at each pixel. You now want to extract overall information of the picture, and be able to provide the approximate intensities in any requested coordinate. That is, initially you can only request pixel by pixel to see their intensities (with the cost of time for each viewed pixel). The first step is then to efficiently re-construct the picture with good enough accuracy (say, at 1% relative accuracy in intensity). The second step is to then use your re-constructed (approximate) picture to be able to quickly find the intensity in any new specific coordinate. I understand that there is no perfect way (and we have to assume that the pictures are fairly smooth and reasonably well behaved) and, as I tried to illustrate, matlab's fsurf seems to do the job but it is not documented how it does it so I can not know how/if I can control its results.
• I will look more into the links you provided, but at a first glance it seems to be mostly only different interpolation routines -- and my guess is that I need to use something as simple as linear interpolation to not loose speed in evaluating the approximate values.

Many thanks for the feedback and the routines, I will think more about if any of your input can help.

[Bjorn Gustavsson]

You have to understand that it is "very challenging" to write a function that finds the global minima, or even local, of an arbitrary function numerically. Since the function might have a very steep minima in a very small region of the independent parameters (think about at 2-D Gaussian with very very small width) that the optimization-function never can be guaranteed to get close enough to be attracted to.

What you're asking for is one step worse - because you'll need to find all local extrema to model the variation between them. Since you might have infinite number of those this is an impossible task as expressed. Let me give you an example to ponder:

function Z = tricky_bstd(x,y)
Z = sin(1./x) + cos(1./y);
end

Think about how to approximate that one in the region [0,1]x[0,1]...

[MG]

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]

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