I have a two-dimensional function $f(x,y)$ whose values I would like to sample. The function is very expensive to compute and it has a complex shape, so I need to find a way to get the most information about its shape using the least number of sample points.

What good methods are there to do this?

What I have so far

  • I start from an existing set of points where I have already computed the function value (this could be a square lattice of points or something else).

  • Then I compute a Delaunay triangulation of these points.

  • If two neighbouring points in the Delaunay triangulation are far enough ($ > \Delta X$) and the function value differs sufficiently in them ($> \Delta f$), then I insert a new point midway inbetween them. I do this for each neighbouring point-pair.

What's wrong with this method?

Well, it works relatively well, but on functions similar to this one it's not ideal because sample points tend to "jump over" the ridge and not even notice it's there.

Mathematica graphics

It produces results like this (if the resolution of the initial point grid is sufficiently rough):

Mathematica graphics

This plot above shows the points where the function value is calculated (actually Voronoi cells around them).

Mathematica graphics

This plot above shows the linear interpolation generated from the same points, and compares it to Mathematica's built-in sampling method (for about the same starting resolution).

How to improve it?

I think the main issue here is that my method decides whether to add a refinement point or not based on the gradient.

It would be better to take into account the curvature or at least the second derivative when adding refinement points.


What is a very simple to implement way to take into account the second derivative or curvature when the locations of my points are not constrained at all? (I don't necessarily have a square lattice of starting points, this should ideally be general.)

Or what other simple ways are there to calculate the position of refinement points in an optimal way?

I am going to implement this in Mathematica, but this question is mainly about the method. For the "easy to implement" bit it does count that I'm using Mathematica though (i.e. this was easy to do so far because it has a package for doing Delaunay triangulation)

What practical problem I am applying this to

I am calculating a phase diagram. It has a complex shape. In one region its value is 0, in another region it's between 0 and 1. There's a sharp jump between the two regions (it's discontinuous). In the region where the function is greater than zero there is both some smooth variation and a couple of discontinuities.

The function value is calculated based on a Monte Carlo simulation, so occasionally an incorrect function value or noise is to be expected (this very is rare, but for a large number of points it happens, e.g. when the steady state is not reached to due some random factor)

I have asked this already on Mathematica.SE but I can't link to it because it's still in private beta. This question here is about the method, not the implementation.

Reply to @suki

Is this the type of division you suggest, i.e. putting a new point in the middle of the triangles?

Mathematica graphics Mathematica graphics Mathematica graphics Mathematica graphics

My concern here is that it seems to require special handling at the edges of the region, otherwise it will give very long and very thin triangles, as shown above. Did you correct for this?


A problem which appear both with the method I describe and with @suki's suggestion to put subdivide based on triangles and put the subdivision points inside the triangle is that when there are discontinuities (as in my problem), recomputing the Delaunay triangulation after a step may cause triangles to change and perhaps some big triangles to appear which have different function values in the three vertices.

Here are two examples:

ex1 ex2

The first shows the end result when sampling around a straight discontinuity. The second shows the sampling point distribution for a similar case.

What simple ways are there to avoid this? Currently I am simply subdividing those egdes that disappear after a retriangulation, but this feels like a hack and needs to be done with care as in the case of symmetrical meshes (like a square grid) there are several valid Delaunay triangulations, hence edges might change randomly after retriangulation.

  • $\begingroup$ are there any new developments on this issue? $\endgroup$
    – Andrei
    Jun 20, 2014 at 22:16

2 Answers 2


I worked on a problem similar to this a while back.

I think the main difference between between our implementations is that I was choosing where to add points based on the triangles, not the edges. I also choose new points inside the triangles instead of on the edges.

I have the feeling that adding points inside the triangles would make it more efficient by giving a small increase in the average distance from old points to the new one.

Anyway another nice thing about using triangles instead of edges is that it gives an estimate of the gradient vector, instead of the slope along this particular edge.

In my matlab code I used a base class to take care of most of the machinery, with a few abstract methods:

  • weight(self) to decided the priority for which triangles to subdivide next.
  • choosePoints(self,npoints = "auto") to decide new points to evaluate based on the weight of each triangle.

I found this setup very flexible:

  • setting a sub-class's weight() function to the area of the triangle produces a constant mesh density.
  • setting weight() to calculate the average function value times the area of the triangle gives a sort of quasi-random probability sampling.
  • using var(triangle.zs) could do, for functions that have binary output, what I feel is a generalization of a bisection search to more than 1 dimension.
  • using area + var(triangle.zs) was pretty effective at putting a constant density everywhere, and an increased density along any slope (almost what you have now).

I used the variance of the z values, to approximate the importance of first order effects (slope) because the variance will never go to infinity like the slope can.

For the last example the background density was good because I was searching for discontinuous blobs of high value in a space of low value. So it would slowly fill in the entire mesh and when it would find a blob it would concentrate on following the edge of the blob all the way around because of the high weight I put on the gradient (and that it only filled in the top n triangles on each iteration). At the end I could know there were no (reasonably shaped) blobs (or holes in my blobs) of a size greater than the resulting background mesh density.

Like you I did get some bad points in my results, they weren't a problem for me because the error was such that if you re-ran nearby points they would probably give the correct answer. I'd just end up with blips of increased mesh density around my bad points.

Whatever you do I always recommend making the weights related to the triangle size so that, all else being equal, large triangles are broken up first.

Maybe a solution for you is to take my approach one step further and instead of evaluating triangles based on the contents of that triangular cell, evaluate based on that one and all three adjacent triangles.

That will contain enough information to get an estimate of the complete Hessian matrix. you can get it by doing a least squares fit of z = c1*x + C2*y c11*x^2+c12*x*y+c22*y^2 over all the vertexes in the triangles of interest (center the coordinate system on the triangle first).

I wouldn't use the gradient or Hessian (those constants) directly because they'll go to infinity at a discontinuity.

Maybe the sum-square-error of the z values relative to a planar approximation of those points would be a useful measure of how interesting second order effects are.


That looks reasonable to me.

I never actually got around to special casing the edges. It bothered me a little but for what I was doing, it was enough to just start with a lot of points around the edges.

more elegant would be to combine our two approaches, weighting edges and triangles. Then if an edge is too long, cut it in half... I like the way that concept generalizes to higher dimensions (but the numbers get big fast)...

But since you don't expect the main body of the mesh to have high aspect ratio triangles, so you could use a function like Matlab's freeboundary function to find the boundary, then run the the same algorithm in one less dimension on the boundary. If done right, on a cube for example, you could get the same mesh density on the edges, on the faces, and inside the cube. Interesting.

One thing that I never found a good solution for was the fact that my version would never explore outside the convex hull of the initial point set.

  • $\begingroup$ I also thought about using triangles first, but I had some technical problem first (which I have solved since), and later I thought it wouldn't be much better to use triangles anyway. Question: where do you put the new points? In the middle of the triangles? I didn't do this because I expected it'd create some very long and thin triangles. I'll update my post shortly with what I understand you did so you can verify if I got it right :-) thanks! $\endgroup$
    – Szabolcs
    Jan 22, 2012 at 19:29
  • $\begingroup$ Can you please see my edit and clarify? $\endgroup$
    – Szabolcs
    Jan 22, 2012 at 21:53
  • $\begingroup$ It turns out special casing the edges is unavoidable, no matter what sort of subdivision scheme I use. In my case I have a high gradient perpendicular the edge, but not parallel to it, which was making things inefficient if I didn't special case the edges. $\endgroup$
    – Szabolcs
    Feb 28, 2012 at 12:53
  • $\begingroup$ Another problem I found was that re-triangulation caused big triangles to appear occasionally where the vertices had different function values. I ended up with things like this: i.stack.imgur.com/nRPwi.png is the linearly interpolated density plot, and i.stack.imgur.com/208bP.png is the sampling points (not exactly the same). This is just a discontinuity along a straight edge. Did you hit this problem? If yes, how did you solve it? Did you retriangulate completely after each subdivision step? $\endgroup$
    – Szabolcs
    Feb 28, 2012 at 12:58
  • $\begingroup$ I am not sure the triangulation really means anything here. Each point that you assessed is the value of the function at a point, so why not do something like they use in mesh-free methods? en.wikipedia.org/wiki/Smoothed-particle_hydrodynamics You can also estimate the derivatives this way... $\endgroup$
    – meawoppl
    Mar 3, 2012 at 0:38

I think the main problem in your heuristic is that you're considering the gradient in one dimension only and thus, in regions where dfdx is small but dfdy is large (as happens in the middle of your example), you will miss points when looking in the "wrong" dimension.

One quick fix would be to consider sets of four points, taking their center of gravity, and approximating |dfdx|+|dfdy| using those four points. Another alternative is to take three points (i.e. a triangle) and take the maximum gradient of the surface over them.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.