Computing curvatures with CGAL

I am trying to understand how I can compute the principal curvatures of a SurfaceMesh (and so the gaussian one) using CGAL.

First of all, since I didn't see it, is there any function that computes a curvature for all vertices in a mesh?

There seems to be an example in Jet_fitting_3, but it's really convoluted, with several additional classes...

As far as I understand, I have to, for each vertex:

1. gather points around it
2. use monge_fit on them obtaining a monge form

But the example then modifies the normal with comply_wrt_given_normal which I do not understand completely why it would be needed.

I might also be skipping over a quite obvious function that computes this, as I said before... after all it's really an important property.

The code I have succeeded to write is this:

for (auto v : result.vertices())
{
    adj.clear();
    
    auto h = result.halfedge(v);
    
    // find all adjacent vertices
    for (HalfedgeDescriptor he : result.halfedges_around_target(h))
        adj.push_back(result.point(result.source(he)));
            
    CGAL::Monge_via_jet_fitting<Kernel> monge_fit;
    CGAL::Monge_via_jet_fitting<Kernel>::Monge_form monge_form;
    
    monge_form = monge_fit(adj.begin(), adj.end(), 2, 1);
    monge_form.comply_wrt_given_normal(vn.first[v]);
    
    // get the principal curvatures
    double k1 = monge_form.principal_curvatures(0);
    double k2 = monge_form.principal_curvatures(1);
}

However, it raises an exception in calling monge_fit, with the following text:

CGAL error: precondition violation!
Expression : nb_input_pts >= nb_d_jet_coeff
File       : /opt/homebrew/include/CGAL/Monge_via_jet_fitting.h
Line       : 297
Explanation: 
Refer to the bug-reporting instructions at https://www.cgal.org/bug_report.html
libc++abi: terminating due to uncaught exception of type CGAL::Precondition_exception: CGAL ERROR: precondition violation!
Expr: nb_input_pts >= nb_d_jet_coeff
File: /opt/homebrew/include/CGAL/Monge_via_jet_fitting.h
Line: 297

So the number of points are not sufficient, I guess?

I don't know, however, how to correct this... any hints are more than welcome.

The method comply_wrt_given_normal 'orients' the fitted jet.
Fitting a jet by only giving points gives curvatures values and a normal vector that are valid up to a sign.
If you care about the sign of the Gaussian curvature then you must call comply_wrt_given_normal, else you can avoid calling it and use the absolute value of the curvature.
Internally, comply_wrt_given_normal just swap and multiply by -1 some coefficients of the jet.

Fitting a jet requires to solve a system of equations where the dimensions depend on the number of points and the order of the jet.
With a mesh, you have very few adjacent vertices so the system of equation is under-determined and cannot be solved, hence the error you get.

I see 4 solutions

• lower the degree of the jet: in the example https://doc.cgal.org/latest/Jet_fitting_3/index.html#Jet_fitting_3SingleEstimationaboutaPoint, d_fitting and d_monge are set to 4, but since you are only interested in principal curvatures (2nd order derivatives of the surface), then you can set them to 2. Order 2 may be sufficient for a few points, but it will depends on your mesh connectivity.
• increase the number of points: add other points to the vertices you already give to the jet fitting, like the center point of each adjacent triangles or incident edges. It's a little bit 'artificial' and it may produces different curvatures values, but that can solve your issue. A complex but accurate solution would be to uniformly sample each adjacent triangles with a sufficiently high number of points.
• use another library if that is possible:
  • libigl https://libigl.github.io/tutorial/#gaussian-curvature
  • ponca https://poncateam.github.io/ponca/example_cxx_basic_page.html
  • libpcl https://pcl.readthedocs.io/en/latest/normal_estimation.html
  • vcglib (Meshlab) https://github.com/cnr-isti-vclab/vcglib (https://www.meshlab.net)
• approximate yourself the Gaussian curvature 'manually' on each vertex. For a given vertex, compute the sum of angles between pair of incident and successive edges. Then the Gaussian curvature can be approximated by the difference between 2*pi and this sum divided by the sum of areas of adjacent triangles (this is a formula I remembered from memory so it should be double checked).