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.
It produces results like this (if the resolution of the initial point grid is sufficiently rough):
This plot above shows the points where the function value is calculated (actually Voronoi cells around them).
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?
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:
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.