A way to avoid stack overflow in CVT

[daichi-ishida]

Hi Bruno, this library is very practical and I've been using this a lot!

I am currently facing an issue with a stack overflow exception when setting up a large-scale point cloud with around 650,000 points in CVT. Further investigation suggests that the NearestNeighbors and KdTree used in the Delaunay_NearestNeighbors class are consuming the stack excessively. Are there any options available to address this problem and prevent the stack overflow exception?

This is how I set now (almost same as vorpalite)

    GEO::CentroidalVoronoiTesselation CVT(&M_domain, GEO::coord_index_t(3));
    CVT.set_volumetric(true);
    CVT.delaunay()->set_vertices(points.size()/3, points.data());
    CVT.RVD()->set_exact_predicates(true);

    GEO::BuildRVDMesh callback(M_result);
    callback.set_simplify_voronoi_facets(true);

    double angle_threshold = GEO::CmdLine::get_arg_double("poly:normal_angle_threshold");
    if (angle_threshold != 0.0) callback.set_simplify_boundary_facets(true, angle_threshold);

    callback.set_tessellate_non_convex_facets(GEO::CmdLine::get_arg_bool("poly:tessellate_non_convex_facets"));
    callback.set_shrink(GEO::CmdLine::get_arg_double("poly:cells_shrink"));
    callback.set_generate_ids(true);
    CVT.RVD()->for_each_polyhedron(callback);

Many thanks

[BrunoLevy]

Hi,

• on which platform are you running ? (I know Windows has by default a smaller stack size than Linux, but it can be increased)
• Would it be possible for you to share your dataset and example programs ? Then it will be easier for us to figure out what happens. Best wishes, -- Bruno

[daichi-ishida]

The platform is Windows, but as for the RVD build, running it on Linux does not finish the process as well. Could it be related to the point cloud being too dense? (As for point cloud, there no duplicated points)

Here is my dataset and example program: https://github.com/daichi-ishida/geogram_example

The data directory contains success case and failure case (that causes stack overflow). You can switch them on lines 175 and 176 in geogram_example/src/point_check/main.cpp

Thanks in advance

[BrunoLevy]

Fantastic ! Thank you very much for this detailed issue report (the best one I have ever received !).

[BrunoLevy]

I was able to get it and compile it. I will come back to you shortly, as soon as I make progress.

[BrunoLevy]

Got it ! Pretty nice example :-)

To obtain it, I changed two things: 1) insert this line right right after all the other GEO::CmdLine(...) calls:

GEO::CmdLine::set_arg("algo:delaunay","PDEL");

2) replace the call to preprocess(M_domain) with:

M_domain.facets.triangulate();

Explanations:

• The default algorithm for RVD computation uses a KDTree, which is in general faster than computing a Delaunay triangulation. In the case of your second example, it has a very sharp variation of pointset density (around the border of the motorcycle). For this reason, the RVD algorithm needs to query a large number of nearest neighbors. Since they are allocated on the stack, it is what caused the stack overflow. The line GEO::CmdLine::set_arg("algo:delaunay","PDEL") tells the system to use a standard Delaunay triangulation algorithm instead of nearest neighbor search (and the P is for Parallel).

• It seems that the pre-processing function is bugged. Need to investigate more (it generates a kind of weird pyramid instead of keeping the original box). Anyway, we do not really need it. We just need to make sure the facets are triangulated before tetrahedralizing it.

[BrunoLevy]

Just filed a pull request to https://github.com/daichi-ishida/geogram_example with the fix.

[daichi-ishida]

Thank you for your quick response and for resolving the issue! I have also tested and confirmed that it works on my end as well, so I have merged the Pull Request.