Dead loop when inserting points to resolve an iso surface in a P3 Delaunay mesh

[gy114]

Hi,

I've been using the following code to refine a periodic 3D Delaunay mesh. This refined mesh will be later used to extract the iso surface of periodic_function(x,y,z) = 0. The idea is to insert points whenever a unique edge is detected to be "close enough" to the iso surface. The inserted points are at the center of that edge, and the info() is updated. This refinement can go through many rounds.

I found that in the middle of refinement, the code compiled in Release mode gets stuck somewhere when new points are inserted.
The Debug mode sometimes gives the following error:

CGAL error: precondition violation!
Expression : side_of_cell(tester.point(),o, c, lt_assert, i_assert, j_assert) != ON_UNBOUNDED_SIDE
File       : /opt/cgal/include/CGAL/Periodic_3_triangulation_3.h
Line       : 2690
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: side_of_cell(tester.point(),o, c, lt_assert, i_assert, j_assert) != ON_UNBOUNDED_SIDE
File: /opt/cgal/include/CGAL/Periodic_3_triangulation_3.h
Line: 2690

Sometimes the Debug mode can also be stuck when inserting the new points.

Please help.

Thanks.

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h>
#include <CGAL/Periodic_3_triangulation_ds_vertex_base_3.h>
#include <CGAL/Periodic_3_triangulation_ds_cell_base_3.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
#include <CGAL/Triangulation_cell_base_with_info_3.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
#include <CGAL/Triangulation_data_structure_3.h>

#include <CGAL/Polyhedron_3.h>
#include <CGAL/number_type_config.h> // CGAL_PI
#include <random>


typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Periodic_3_Delaunay_triangulation_traits_3<K> P3_Traits;
typedef CGAL::Periodic_3_triangulation_ds_vertex_base_3<> P3_VbDS;
typedef CGAL::Triangulation_vertex_base_3<P3_Traits, P3_VbDS >       P3_Vb;
typedef CGAL::Periodic_3_triangulation_ds_cell_base_3<> P3_CbDS;
typedef CGAL::Triangulation_cell_base_3<P3_Traits, P3_CbDS >       P3_Cb;


typedef CGAL::Triangulation_vertex_base_with_info_3<double, P3_Traits, P3_Vb > VbInfo; // Test to attach data to vertices.
typedef CGAL::Triangulation_cell_base_with_info_3<double, P3_Traits, P3_Cb > CbInfo; // Test to attach data to cells.
typedef CGAL::Triangulation_data_structure_3<VbInfo, CbInfo >                 TDS;
typedef CGAL::Periodic_3_Delaunay_triangulation_3<P3_Traits, TDS >               P3_Delaunay;
typedef P3_Delaunay::Point P3_Point;
typedef P3_Delaunay::Periodic_point PeriodicPoint;
typedef P3_Delaunay::Iso_cuboid Iso_cuboid;




double periodic_function(const P3_Point& p) {
    static unsigned long callCount = 0;
    callCount++;
    return -std::cos(2 * CGAL_PI * p[0]) - std::cos(2 * CGAL_PI * p[1]) - std::cos(2 * CGAL_PI * p[2]);
}


int main() {
    Iso_cuboid domain(P3_Point(-0.5, -0.5, -0.5), P3_Point(0.5, 0.5, 0.5));

    auto dx = domain.xmax() - domain.xmin();
    auto dy = domain.ymax() - domain.ymin();
    auto dz = domain.zmax() - domain.zmin();

    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<> disX(domain.xmin(), domain.xmax());
    std::uniform_real_distribution<> disY(domain.ymin(), domain.ymax());
    std::uniform_real_distribution<> disZ(domain.zmin(), domain.zmax());

    P3_Delaunay T(domain);


    std::vector<P3_Point> listOfRandPoint;
    for (unsigned long long i=0; i<1000; i++)  {
        P3_Point randPoint(disX(gen), disY(gen), disZ(gen));
        listOfRandPoint.push_back(randPoint);
    }

    for (unsigned long long i=0; i<listOfRandPoint.size(); i++)  {
        P3_Delaunay::Vertex_handle vh = T.insert(listOfRandPoint[i]);
        vh->info() = periodic_function(listOfRandPoint[i]);
        std::cout << "\r" << i << " points inserted." << std::flush;
    }

    // Now do refinement for 5 times.
    int MaxRoundOfRefinement = 5;
    for (int iTimes = 0; iTimes < MaxRoundOfRefinement; iTimes++)    {
        std::vector<K::Point_3> points_to_insert;
        for (auto eit = T.unique_edges_begin(); eit != T.unique_edges_end(); ++eit) {
            auto edge = *eit;
            auto cell = edge.first;
            int i = edge.second;
            int j = edge.third;

            auto v1 = cell->vertex(i);
            auto v2 = cell->vertex(j);

            double f1 = v1->info();
            double f2 = v2->info();


            // Calculate t for the potential zero-crossing point
            double t;
            if (f2 != f1) { // Ensure denominator is not zero
                t = (0 - f1) / (f2 - f1);
            } else {
                // If f1 == f2, we cannot determine a zero-crossing in the traditional sense
                // Skip refinement for this edge, or handle as special case
                continue; // Move to the next edge
            }


            // Check if the zero-crossing is within the edge or within the extended range
            if (-2.0 < t && t < 3.0) {
                // The zero-crossing point is either directly on the edge or close enough for refinement
                PeriodicPoint pp1 = T.periodic_point(cell, i);
                PeriodicPoint pp2 = T.periodic_point(cell, j);
                P3_Delaunay::Offset offset1 = pp1.second;
                P3_Delaunay::Offset offset2 = pp2.second;

                P3_Point canonical_p1 = v1->point();
                P3_Point canonical_p2 = v2->point();

                double x1 = CGAL::to_double(canonical_p1.x()) + offset1[0] * dx;
                double y1 = CGAL::to_double(canonical_p1.y()) + offset1[1] * dy;
                double z1 = CGAL::to_double(canonical_p1.z()) + offset1[2] * dz;

                double x2 = CGAL::to_double(canonical_p2.x()) + offset2[0] * dx;
                double y2 = CGAL::to_double(canonical_p2.y()) + offset2[1] * dy;
                double z2 = CGAL::to_double(canonical_p2.z()) + offset2[2] * dz;

                K::Point_3 new_point = P3_Point( (x1 + x2)*0.5, (y1+y2)*0.5, (z1+z2)*0.5 );
                P3_Traits traits(T.domain());
                K::Point_3 canonical_point = CGAL::P3T3::internal::robust_canonicalize_point(new_point, traits);
                points_to_insert.push_back(canonical_point);
            }
        }
        std::cout << "\nRound "<< iTimes << " " << points_to_insert.size() << " new points found." << std::endl;

        for (unsigned long long i=0; i<points_to_insert.size(); i++)   {
            auto vh = T.insert(points_to_insert[i]);
            vh->info() = periodic_function(vh->point());
        }
        std::cout << "all points inserted." << std::endl;
    }
    std::cout << MaxRoundOfRefinement << " rounds of refinement finished successfully." << std::endl;
    return 0;
}

[gy114]

I found a work around by inserting new points at the centroids of unique cells instead of the mid points of unique edges. Now it seems to be working. I guess inserting new points in the middle of edges makes it a bad Delaunay mesh.

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h>
#include <CGAL/Periodic_3_triangulation_ds_vertex_base_3.h>
#include <CGAL/Periodic_3_triangulation_ds_cell_base_3.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
#include <CGAL/Triangulation_cell_base_with_info_3.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
#include <CGAL/Triangulation_data_structure_3.h>

#include <CGAL/Polyhedron_3.h>
#include <CGAL/number_type_config.h> // CGAL_PI
#include <random>


typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Periodic_3_Delaunay_triangulation_traits_3<K> P3_Traits;
typedef CGAL::Periodic_3_triangulation_ds_vertex_base_3<> P3_VbDS;
typedef CGAL::Triangulation_vertex_base_3<P3_Traits, P3_VbDS >       P3_Vb;
typedef CGAL::Periodic_3_triangulation_ds_cell_base_3<> P3_CbDS;
typedef CGAL::Triangulation_cell_base_3<P3_Traits, P3_CbDS >       P3_Cb;


typedef CGAL::Triangulation_vertex_base_with_info_3<double, P3_Traits, P3_Vb > VbInfo; // Test to attach data to vertices.
typedef CGAL::Triangulation_cell_base_with_info_3<double, P3_Traits, P3_Cb > CbInfo; // Test to attach data to cells.
typedef CGAL::Triangulation_data_structure_3<VbInfo, CbInfo >                 TDS;
typedef CGAL::Periodic_3_Delaunay_triangulation_3<P3_Traits, TDS >               P3_Delaunay;
typedef P3_Delaunay::Point P3_Point;
typedef P3_Delaunay::Periodic_point PeriodicPoint;
typedef P3_Delaunay::Iso_cuboid Iso_cuboid;




double periodic_function(const P3_Point& p) {
    static unsigned long callCount = 0;
    callCount++;
    return -std::cos(2 * CGAL_PI * p[0]) - std::cos(2 * CGAL_PI * p[1]) - std::cos(2 * CGAL_PI * p[2]);
}


int main() {
    Iso_cuboid domain(P3_Point(-0.5, -0.5, -0.5), P3_Point(0.5, 0.5, 0.5));

    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<> disX(domain.xmin(), domain.xmax());
    std::uniform_real_distribution<> disY(domain.ymin(), domain.ymax());
    std::uniform_real_distribution<> disZ(domain.zmin(), domain.zmax());

    P3_Delaunay T(domain);


    std::vector<P3_Point> listOfRandPoint;
    for (unsigned long long i=0; i<5000; i++)  {
        P3_Point randPoint(disX(gen), disY(gen), disZ(gen));
        listOfRandPoint.push_back(randPoint);
    }

    for (unsigned long long i=0; i<listOfRandPoint.size(); i++)  {
        P3_Delaunay::Vertex_handle vh = T.insert(listOfRandPoint[i]);
        vh->info() = periodic_function(listOfRandPoint[i]);
        std::cout << "\r" << i << " points inserted." << std::flush;
    }

    // Now do refinement for 10 times.
    int MaxRoundOfRefinement = 10;
    for (int iTimes = 0; iTimes < MaxRoundOfRefinement; iTimes++)    {
        std::vector<P3_Point> points_to_insert;

        // The following is a new strategy to refine the mesh: inserting points at the centroids of all cells if they are close to the iso surface.
        for (auto cit = T.unique_cells_begin(); cit != T.unique_cells_end(); ++cit) {
            double max_f = -std::numeric_limits<double>::infinity();
            double min_f = std::numeric_limits<double>::infinity();

            // Find max_f and min_f within the cell
            for (int i = 0; i < 4; ++i) { // Each cell has 4 vertices
                auto vh = cit->vertex(i);
                double f = vh->info();
                max_f = std::max(max_f, f);
                min_f = std::min(min_f, f);
            }

            // Calculate t using max_f and min_f
            double t;
            if (max_f != min_f) { // Ensure denominator is not zero
                t = (0 - min_f) / (max_f - min_f);
            } else {
                // If max_f == min_f, we cannot determine a zero-crossing in the traditional sense
                continue; // Move to the next cell
            }

            // Use t value for further processing, if needed
            
            if (t > -0.2 && t < 0.3) {
                // Calculate the "real" positions of the vertices considering their offsets
                double sum_x = 0, sum_y = 0, sum_z = 0;
                int count = 0;

                for (int i = 0; i < 4; ++i) { // Each cell has 4 vertices
                    PeriodicPoint pp = T.periodic_point(cit, i);
                    P3_Point real_point = T.point(pp); // Converts periodic point to its real position
                    sum_x += real_point.x();
                    sum_y += real_point.y();
                    sum_z += real_point.z();
                    count++;
                }

                // Calculate the centroid coordinates
                double centroid_x = sum_x / count;
                double centroid_y = sum_y / count;
                double centroid_z = sum_z / count;
                P3_Point centroid(centroid_x, centroid_y, centroid_z);
                P3_Traits traits(T.domain());
                P3_Point new_point = CGAL::P3T3::internal::robust_canonicalize_point(centroid, traits);
                new_point = T.construct_point(new_point);
                points_to_insert.push_back(new_point);
            }
        }

        std::cout << "\nRound "<< iTimes << " " << points_to_insert.size() << " new points found." << std::endl;

        std::cout << "inserting... " << std::flush;
        for (unsigned long long i=0; i<points_to_insert.size(); i++)   {
            try {
                auto vh = T.insert(points_to_insert[i]);
                vh->info() = periodic_function(vh->point());
            }
            catch (const std::exception& e)  {
                std::cerr << "An error occured during insertion: " << e.what() << std::endl;
            }
        }
        std::cout << "done." << std::endl;
    }
    std::cout << MaxRoundOfRefinement << " rounds of refinement finished successfully." << std::endl;
    return 0;
}