Search code examples
c++geometrycomputational-geometrymeshcgal

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.

Solution

  • 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