I would like to find the contours of a scalar function $u(x,y)$ available as a discrete set of function values $u_i = u(x_i,y_i)$ over a scattered set of points $(x_i,y_i), i=1,...,N$. In my case, the scattered data comes from a numerical solution to a PDE using the RBF finite difference method, however interpolation of any type of scattered data is just as relevant.

In a related thread, the recommendation was to form a triangulation of the scattered points first, and then using a common contour algorithm such as marching triangles in 2D. While this is a tested approach, building the triangulation defeats one of the primary strengths of the RBF-FD method, which is exactly the fact it can avoid polygonal meshes. One of the answers in that thread suggested a RBF approach could work better.

Apart from the point coordinates and function values, suppose I also have the adjacency graph for the RBF-FD stencils of size $s$, and the RBF-FD differentation weights $w_i$:

$$ L(u(\boldsymbol{x}))|_{\boldsymbol{x}_c=\boldsymbol{x}_1} \approx \sum_{j=1}^s w_j u_j $$

which are constructed in a way, that they can provide approximations to various (linear) spatial differential operators $\partial_x, \partial_y, ...$.

Returning to the initial problem, finding a contour level $h$ can be recast to finding the level set (zero contour) of the implicit function,

$$ \phi(x,y) = u(x,y) - h = 0, $$ which can be seen as a signed distance-like function (except for the fact it doesn't really return the distances).

Given some initial guess $\boldsymbol{p} = (p_x, p_y)$ of a point near the contour, we can try and find a projection $\Delta\boldsymbol{p}$ so that

$$ \phi(\boldsymbol{p} + \Delta\boldsymbol{p}) = 0, \\ \Delta \boldsymbol{p} \,|| \, \nabla\phi(\boldsymbol{p}+\Delta \boldsymbol{p}), $$

where the projection is to be orthogonal to the gradient at the projected point. In the thesis of Per-Olof Persson (pp. 47-48) a first order approximation of the projection is found using Lagrange multipliers and a truncated Taylor expansion. The approximate projection formula is

$$ \Delta \boldsymbol{p} = \frac{\phi}{|\nabla \phi|^2} \nabla \phi. $$

An example of the approximate projection result is shown in the figure below. First, I identified the scattered points near the contour by looking for stencils, where the sign of $\phi$ changes. Using the RBF-FD spatial operators to construct the gradient $\nabla = (\partial_x, \partial_y)$, I could project the nearby points to the target contour:

Approximate projection

What I am missing now is a robust way to connect the projected points into a single contour. Is their some graph approach that could be used (construct the shortest path)?

In view of the marching squares algorithm, I am also anticipating the path-finding will have disambiguation difficulties near saddle points. Is their any way to detect such saddle regions and if yes, how to make the path-finding robust?

Approximation near saddle

  • $\begingroup$ Your introduction is quite long. Do you want to find an algorithm that connects the blue/green points with the shortest path? If yes, the problem is a NP hard problem and can be solved with the traveling salesman algorithm. $\endgroup$
    – ConvexHull
    Commented Nov 28, 2020 at 7:47
  • $\begingroup$ Precisely. I was hoping that by using the adjacency graph of the points near the contour (e.g. to zigzag between +/x symbols) it might be possible to solve this more elegantly. $\endgroup$
    – IPribec
    Commented Nov 28, 2020 at 11:14
  • $\begingroup$ I think for this you have define your problem more precisly. Most of the stuff you mentioned is not needed (my point of view) to describe your actual problem. Are the adjacency points ordered in some kind? For what kind of application do you need the exact contour line? $\endgroup$
    – ConvexHull
    Commented Nov 28, 2020 at 15:00
  • 1
    $\begingroup$ I’m not sure I understand the problem. Couldn’t you just fit a curve to the points? $\endgroup$ Commented May 1, 2021 at 13:06
  • $\begingroup$ What curve do you suggest fitting? The contour is generally an implicit function, $f(x,y)=0$. $\endgroup$
    – IPribec
    Commented Jul 29, 2023 at 21:30

2 Answers 2


Maybe this is not a full answer to your question. However, I am currently developing my own unstructured mesh generator and found this toolbox quite helpful. Perhaps there are some algorithms which will help you.

Matlab or Github

Short (not complete) overview of the genetic algorithms for the TPS problem:

  • Solves the classic Traveling Salesman Problem (TSP)
  • Solves an open variation of the TSP with no start/end constraints
  • Solves an open variation of the TSP with fixed start & end constraints
  • Solves an open variation of the TSP with only a fixed start constraint
  • Solves a variation of the TSP that maximizes the total distance
  • Nearest Neighbor (NN) solution to the TSP
  • Nearest Neighbor (NN) solution to the open variation of the TSP (no start/end constraints)
  • Farthest Neighbor (FN) solution to the TSP
  • Random Search (RS) solution to the TSP (provides interesting performance comparisons)
  • Solves the TSP with a variation of the GA that uses a cross-over operator
  • Solves the TSP with a variation of the GA that uses a hybrid set of operators
  • Variation of the TSP to find the shortest tour through clusters of points
  • Solves the TSP with a variation of the GA that increases the mutation rate on the best route
  • Solves the Open TSP with a variation of the GA that increases the mutation rate on the best route
  • Script to compare tsp_ga performance with an approach that seeds the algorithm with the NN solution
  • Script to compare tsp_ga_turbo performance with an approach that seeds the algorithm with the NN solution
  • Solves the classic Multiple Traveling Salesmen Problem (M-TSP)
  • Random Search (RS) solution to the M-TSP (provides interesting performance comparisons)
  • Solves a variation of the M-TSP where all salesmen start and end at the first city
  • ...

The following paper does a decent job of answering my question:

D. H. McLain, Drawing Contours from Arbitrary Data Points, The Computer Journal, Volume 17, Issue 4, November 1974, Pages 318–324, https://doi.org/10.1093/comjnl/17.4.318

The main ideas I found useful are 1) split the domain into small rectangles or windows, 2) apply classic grid based interpolation methods in the rectangles - this makes drawing the contours faster and easier, 3) use the data from the unstructured points to interpolate the function, and optionally it's derivatives, to the corners of the rectangles.

An important observation is that within a small rectangle (and assuming no degeneracies), the contour will generally not be far from a straight line.

In summary, the approach is to interpolate the data onto a more regular mesh.


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.