Generating set of points on a surface defined by constraint

I'm writing a differential geometry library, and one minor convenience I'd like to offer is to generate a set of points on a surface given by a constraint. For example, for a sphere,

$$x^2+y^2+z^2-r^2=0.$$

Now, if it were an equation in two variables, I could step through one of them and then solve for the other, or choose the first based on some random distribution, and solve for the other.

Here though, I'm not sure what the ideal approach is, and what can be done to make sure I've sort of 'evenly' sampled points on the surface, for triangulation later.

What algorithms would be most appropriate for this task?

• I don't think that your question have a definite answer. Do you want a evenly samples points over the parametric space or over the surface? Dec 12, 2017 at 22:03
• To add to @nicoguaro's question, will you accept points that are within floating-point error close to the surface, or must they be exact in some sense? Dec 12, 2017 at 22:46
• Bear in mind that some closely related problems to yours can be really difficult: IIRC, finding one point on a surface given by the intersection of multivariate polynomial constraints is already at least NP-hard (mathoverflow.net/questions/115608). So depending on how exactly you formulate your question, such an algorithm might very well not exist at all, in which case you'll need to rely on heuristics. Dec 12, 2017 at 23:00
• I don't have a whole answer for you, but if you're interested in ultimately meshing a surface, you might want to check out "Distmesh." It is a finite element mesh generator (for Matlab, but the algorithms are well-documented) that allows you to input a function describing your surface (the 0 level-set of the signed distance function). For more complex geometries, it will solve the Eikonal equation for a numerical approximation of the signed distance function. It may be worth looking at for inspiration, at the very least. Dec 12, 2017 at 23:07
• Maybe checking matching cubes might be useful as well. Dec 13, 2017 at 3:15

If you're interested in surface meshing, the book Delaunay Mesh Generation by Cheng, Dey, and Shewchuk is very good. Surface meshing is a pretty involved topic and I couldn't attempt to summarize it here. The constrained Delaunay triangulation and Voronoi diagram are especially important.

You also mentioned being able to find seed points that are on the surface. Part of this will involve being able to solve nonlinear equations, for which you'll want to use Newton's method. If you're dealing with algebraic surfaces, Smale's $\alpha$-theory can give you convergence guarantees that aren't available in general. My old officemate gave a talk (slides) that I think makes a great introduction. You can also take the sum of the squares of all the equations defining the surface and think of it as an optimization problem, but in that case Smale's theory applies too.

Assuming you do have some point on the surface, another important operation is to move that point following the flow of some vector field. You might also want to be able to compute geodesics. This amounts to solving a differential-algebraic system. Geometric Numerical Integration by Hairer, Lubich and Wanner has a lot of information about ODEs on manifolds.

For a real application, you can check out the manifold class in the deal.II library. deal.II uses manifolds to solve the problem of refining initially coarse meshes on surfaces. If you try to refine a mesh of, say, an annulus, a naive implementation would just put more points on the boundary segments without adjusting the curvature. With deal.II, you can say that a certain part of the boundary of a mesh is actually some manifold, and the refined mesh will adapt to the real surface much better.

Finally, you might find a lot of the engineering literature on boundary representation or constructive solid geometry modelling to be inspirational. Sounds like a cool project and good luck! Do post a link on here if you make any of it open source.

I think what you're looking for is a surface-meshing algorithm. There are a multitude of different approaches and algorithms that can be used for these type of problems - I'll try to (briefly) outline a few options here:

• Parametric space approaches: One approach is to build a mesh in a two-dimensional parametric space associated with the surface. While the meshing problem here becomes a 2D task (an obvious upside!), the construction of robust parametric mappings can be difficult. Such methods typically also require that the parametric space meshes be anisotropic, so that they are 'nicely shaped' when mapped back onto the 3D surface.

• Delaunay-refinement methods: It's also possible to build meshes in 3D space directly; dodging around the parametric mapping problem entirely. One way of doing this is to incrementally build a restricted Delaunay triangulation of the surface - a subset of the triangular faces embedded in a 3D Delaunay tetrahedralisaton that form a good approximation to the surface. Typically, vertices are inserted one-by-one according to a given refinement-strategy until the surface mesh is 'good enough'. One thing that (I think) is nice about these methods is that they come with various provable guarantees on mesh quality, surface approximation error, and etc.

• Variational techniques: Meshing can also be cast as an optimisation problem: find the set of vertex positions and mesh topology that maximise a given mesh-quality functional. Given some initial distribution of vertices, such methods essentially 'slide' vertices over the surface to improve the shape of the mesh; re-computing mesh topology (usually via Delaunay-type methods) as they go. Popular variational methods include Centroidal Voronoi Tessellation (CVT) and Optimal Delaunay Triangulation (ODT).

There are a number of high quality meshing packages that can handle surface problems. To list a few: CGAL, GMSH, GEOGRAM, NETGEN, and (my own library) JIGSAW.