Developing a meshfree contouring algorithm

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 $Deltaboldsymbol{p}$ so that

$$
phi(boldsymbol{p} + Deltaboldsymbol{p}) = 0,
Delta boldsymbol{p} ,|| , nablaphi(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:

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?