Delaunary Triangulation of PMP mesh

[petrasvestartas]

Hi,

I see that dual mesh is not taking boundary vertices for the dual, is there any example to take the boundary middle points into the mesh faces:

std::tuple<compas::RowMatrixXd, std::vector<std::vector<int>>>
pmp_remesh_dual(
    Eigen::Ref<compas::RowMatrixXd> vertices_a,
    Eigen::Ref<compas::RowMatrixXi> faces_a,
    double target_edge_length,
    unsigned int number_of_iterations,
    bool do_project)
{
    // Convert input matrices to CGAL mesh (primal mesh)
    compas::Mesh mesh_a = compas::mesh_from_vertices_and_faces(vertices_a, faces_a);
    // Perform isotropic remeshing
    CGAL::Polygon_mesh_processing::isotropic_remeshing(
        faces(mesh_a),
        target_edge_length,
        mesh_a,
        CGAL::Polygon_mesh_processing::parameters::number_of_iterations(number_of_iterations)
                       .do_project(do_project));

                       
    // Clean up the mesh
    mesh_a.collect_garbage();
    
    // Create the dual mesh
    typedef CGAL::Dual<compas::Mesh> DualMesh;
    DualMesh dual(mesh_a);
    
    
    // Create filtered graph to handle border edges
    typedef boost::filtered_graph<DualMesh, noborder<compas::Mesh>> FiniteDual;
    FiniteDual finite_dual(dual, noborder<compas::Mesh>(mesh_a));
    
    // Create a new mesh for the dual
    compas::Mesh dual_mesh;
    
    // Map to track face-to-vertex conversion
    std::map<typename boost::graph_traits<compas::Mesh>::face_descriptor, 
              typename boost::graph_traits<compas::Mesh>::vertex_descriptor> face_to_vertex;
    
    // 1. For each face in the primal mesh, create a vertex in the dual mesh
    for(auto face : faces(mesh_a)) {
        // Use the centroid of the face as the position for the dual vertex
        auto h = CGAL::halfedge(face, mesh_a);
        auto point = CGAL::centroid(
            mesh_a.point(target(h, mesh_a)),
            mesh_a.point(target(next(h, mesh_a), mesh_a)),
            mesh_a.point(target(next(next(h, mesh_a), mesh_a), mesh_a))
        );
        
        // Add vertex to dual mesh
        auto v = dual_mesh.add_vertex();
        dual_mesh.point(v) = point;
        face_to_vertex[face] = v;
    }
    
    // 2. For each non-border edge in the primal, create an edge in the dual
    for(auto edge : CGAL::make_range(edges(finite_dual))) {
        auto source_face = source(edge, finite_dual);
        auto target_face = target(edge, finite_dual);
        
        // Skip invalid faces
        if(source_face == typename boost::graph_traits<compas::Mesh>::face_descriptor() || 
           target_face == typename boost::graph_traits<compas::Mesh>::face_descriptor())
            continue;
        
        // Add edge connecting the corresponding vertices in the dual mesh
        dual_mesh.add_edge(face_to_vertex[source_face], face_to_vertex[target_face]);
    }
    
    // 3. Create variable-length faces in the dual mesh
    std::map<typename boost::graph_traits<compas::Mesh>::vertex_descriptor, std::vector<int>> vertex_to_face;
    
    // For each vertex in the primal, collect its incident faces to form a face in the dual
    for(auto v : vertices(mesh_a)) {
        std::vector<int> face_indices;
        
        // Collect all faces around this vertex
        for(auto h : CGAL::halfedges_around_target(v, mesh_a)) {
            auto f = face(h, mesh_a);
            // Skip invalid faces
            if(f != boost::graph_traits<compas::Mesh>::null_face()) {
                // Use the corresponding vertex in dual mesh
                if(face_to_vertex.find(f) != face_to_vertex.end()) {
                    // Get the index of the dual vertex
                    auto dual_vertex = face_to_vertex[f];
                    int idx = dual_vertex;
                    face_indices.push_back(idx);
                }
            }
        }
        
        // Only add valid faces (with 3 or more vertices)
        if(face_indices.size() >= 3) {
            vertex_to_face[v] = face_indices;
        }
    }
    
    // Print the dual mesh faces for debugging
    std::cout << "Dual mesh has " << vertex_to_face.size() << " variable-length faces" << std::endl;
    
    // Store the variable-length faces
    std::vector<std::vector<int>> variable_faces;
    for(const auto& pair : vertex_to_face) {
        variable_faces.push_back(pair.second);
    }
    
    // Convert to matrices for regular output
    auto [vertices, triangles] = compas::mesh_to_vertices_and_faces(dual_mesh);
    
    // Return the complete dual mesh information as a tuple
    return std::make_tuple(vertices, variable_faces);
}

[janetournois]

The algorithm to build a Constrained Voronoi Diagram in 2D is available in Constrained_voronoi_diagram_2.h
but not directly applicable to a FaceGraph

[petrasvestartas]

It would be so cool if there would be a dual mesh algorithm that properly handles the boundaries. since it is not a 2D but a triangle mesh :)

[janetournois]

The corresponding paper is "2D Centroidal Voronoi Tessellations with Constraints"

[petrasvestartas]

It was difficult to understand, how to make this work. But here is the solution.
The most annoyint part is that you cannot get order faces around boundary vertices.

std::tuple<compas::RowMatrixXd, std::vector<std::vector<int>>>
pmp_remesh_dual(
    Eigen::Ref<compas::RowMatrixXd> vertices_a,
    Eigen::Ref<compas::RowMatrixXi> faces_a,
    double target_edge_length,
    unsigned int number_of_iterations,
    bool do_project)
{
    /////////////////////////////////////////////////////////////////////////////////
    // Initial Mesh Preparation
    /////////////////////////////////////////////////////////////////////////////////
    // This section initializes the primal mesh and performs remeshing to create
    // a high-quality input for dual mesh generation. The mesh is cleaned up
    // to ensure a consistent structure without any stray elements.
    
    compas::Mesh mesh_a = compas::mesh_from_vertices_and_faces(vertices_a, faces_a); // Convert input matrices to CGAL mesh
    CGAL::Polygon_mesh_processing::isotropic_remeshing(
        faces(mesh_a),
        target_edge_length,
        mesh_a,
        CGAL::Polygon_mesh_processing::parameters::number_of_iterations(number_of_iterations)
                       .do_project(do_project));
    mesh_a.collect_garbage(); // Clean up the mesh
    
    /////////////////////////////////////////////////////////////////////////////////
    // Dual Graph Creation
    /////////////////////////////////////////////////////////////////////////////////
    // Create the dual graph structure and a filtered version that excludes border
    // elements. This filtered graph will be used to generate the interior portion
    // of the dual mesh.
    
    typedef CGAL::Dual<compas::Mesh> DualMesh;
    DualMesh dual(mesh_a);
    
    typedef boost::filtered_graph<DualMesh, noborder<compas::Mesh>> FiniteDual;
    FiniteDual finite_dual(dual, noborder<compas::Mesh>(mesh_a));
    
    compas::Mesh dual_mesh; // Create a new mesh for the dual
    
    std::map<typename boost::graph_traits<compas::Mesh>::face_descriptor, 
              typename boost::graph_traits<compas::Mesh>::vertex_descriptor> face_to_vertex; // Map to track face-to-vertex conversion
    
    /////////////////////////////////////////////////////////////////////////////////
    // Dual Mesh Vertex Creation
    /////////////////////////////////////////////////////////////////////////////////
    // For each face in the primal mesh, create a corresponding vertex in the dual mesh.
    // The position of each dual vertex is the centroid of the corresponding primal face.
    
    for(auto face : faces(mesh_a)) {
        auto h = CGAL::halfedge(face, mesh_a);
        auto point = CGAL::centroid(
            mesh_a.point(target(h, mesh_a)),
            mesh_a.point(target(next(h, mesh_a), mesh_a)),
            mesh_a.point(target(next(next(h, mesh_a), mesh_a), mesh_a))
        );
        
        auto v = dual_mesh.add_vertex(); // Add vertex to dual mesh
        dual_mesh.point(v) = point;
        face_to_vertex[face] = v;
    }
    
    /////////////////////////////////////////////////////////////////////////////////
    // Dual Mesh Edge Creation
    /////////////////////////////////////////////////////////////////////////////////
    // Create edges in the dual mesh for each non-border edge in the filtered dual graph.
    // Each edge in the dual mesh connects two vertices that correspond to adjacent
    // faces in the primal mesh.
    
    for(auto edge : CGAL::make_range(edges(finite_dual))) {
        auto source_face = source(edge, finite_dual);
        auto target_face = target(edge, finite_dual);
        
        dual_mesh.add_edge(face_to_vertex[source_face], face_to_vertex[target_face]);
    }
    
    /////////////////////////////////////////////////////////////////////////////////
    // Dual Mesh Face Creation
    /////////////////////////////////////////////////////////////////////////////////
    // Create variable-length faces in the dual mesh. Each face corresponds to a vertex
    // in the primal mesh and consists of the dual vertices (face centroids) that
    // surround the primal vertex.
    
    std::map<typename boost::graph_traits<compas::Mesh>::vertex_descriptor, std::vector<int>> vertex_to_face;
    
    typename compas::Mesh::Property_map<typename boost::graph_traits<compas::Mesh>::face_descriptor, int> fccmap;
    fccmap = mesh_a.add_property_map<typename boost::graph_traits<compas::Mesh>::face_descriptor, int>("f:CC").first;
    int num_components = connected_components(finite_dual, fccmap);
    
    for(auto v : vertices(mesh_a)) {
        if(is_border(v, mesh_a)) {
            continue; // Skip border vertices
        }
        
        std::vector<int> face_indices;
        
        for(auto h : CGAL::halfedges_around_target(v, mesh_a)) {
            if(is_border(h, mesh_a)) continue;
            
            auto f = face(h, mesh_a);
            if(face_to_vertex.find(f) != face_to_vertex.end()) {
                auto dual_vertex = face_to_vertex[f];
                int idx = dual_vertex;
                face_indices.push_back(idx);
            }
        }
        
        if(face_indices.size() >= 3) {
            vertex_to_face[v] = face_indices; // Only add valid faces (with 3 or more vertices)
        }
    }

    std::vector<std::vector<int>> variable_faces;
    for(const auto& pair : vertex_to_face) {
        variable_faces.push_back(pair.second);
    }
    
    /////////////////////////////////////////////////////////////////////////////////
    // Boundary Processing
    /////////////////////////////////////////////////////////////////////////////////
    // The boundary processing algorithm creates special polylines/faces for boundary vertices:
    // 1. For each boundary vertex in the primal mesh, we create a sequence in the dual mesh
    // 2. This sequence consists of:
    //    a) The boundary vertex itself (same position as in primal)
    //    b) The midpoint of the first boundary edge connected to this vertex
    //    c) The centroids of all faces adjacent to this vertex, sorted by connectivity
    //       from the first boundary edge to the second boundary edge
    //    d) The midpoint of the second boundary edge connected to this vertex
    // 3. This creates a polyline that traces the boundary while connecting to the interior dual mesh
    
    /////////////////////////////////////////////////////////////////////////////////
    // Boundary Edge Midpoint Creation
    /////////////////////////////////////////////////////////////////////////////////
    // Create vertices at the midpoints of all boundary edges. These will be used
    // as the endpoints of the boundary polylines in the dual mesh.
    
    std::map<typename boost::graph_traits<compas::Mesh>::edge_descriptor,
              typename boost::graph_traits<compas::Mesh>::vertex_descriptor> edge_to_midpoint;
              
    for(auto e : edges(mesh_a)) {
        if(is_border(e, mesh_a)) {
            auto h = halfedge(e, mesh_a);
            auto source_v = source(h, mesh_a);
            auto target_v = target(h, mesh_a);
            auto midpoint = CGAL::midpoint(
                mesh_a.point(source_v),
                mesh_a.point(target_v)
            );
            
            auto midpoint_v = dual_mesh.add_vertex(); // Add midpoint vertex to dual mesh
            dual_mesh.point(midpoint_v) = midpoint;
            edge_to_midpoint[e] = midpoint_v;
        }
    }
    
    /////////////////////////////////////////////////////////////////////////////////
    // Boundary Vertex Processing
    /////////////////////////////////////////////////////////////////////////////////
    // For each boundary vertex, create a polyline in the dual mesh that connects
    // the midpoints of its adjacent boundary edges via the centroids of faces
    // adjacent to the boundary vertex.
    
    std::map<typename boost::graph_traits<compas::Mesh>::vertex_descriptor,
              typename boost::graph_traits<compas::Mesh>::vertex_descriptor> vertex_to_dual_vertex;
    
    for(auto v : vertices(mesh_a)) {
        if(is_border(v, mesh_a)) {
            auto dual_v = dual_mesh.add_vertex(); // Create a vertex in dual mesh at the same position
            dual_mesh.point(dual_v) = mesh_a.point(v);
            vertex_to_dual_vertex[v] = dual_v;
            
            std::vector<typename boost::graph_traits<compas::Mesh>::halfedge_descriptor> boundary_halfedges;
            
            for(auto h : CGAL::halfedges_around_target(v, mesh_a)) {
                if(is_border(edge(h, mesh_a), mesh_a)) {
                    boundary_halfedges.push_back(h);
                }
            }
            
            if(boundary_halfedges.size() != 2) {
                continue; // We need exactly two boundary halfedges
            }
            
            auto h1 = boundary_halfedges[0];
            auto h2 = boundary_halfedges[1];
            
            auto e1 = edge(h1, mesh_a);
            auto e2 = edge(h2, mesh_a);
            
            auto midpoint1_v = edge_to_midpoint[e1];
            auto midpoint2_v = edge_to_midpoint[e2];
            
            /////////////////////////////////////////////////////////////////////////////////
            // Face Collection and Edge Extraction
            /////////////////////////////////////////////////////////////////////////////////
            // Collect all faces adjacent to the boundary vertex and extract their edges.
            // This information will be used to determine the adjacency between faces.
            
            std::vector<typename boost::graph_traits<compas::Mesh>::face_descriptor> vertex_faces;
            std::vector<std::vector<size_t>> faces_edges;

            for(auto h : CGAL::halfedges_around_target(v, mesh_a)) {
                auto f = face(h, mesh_a);

                if(f != boost::graph_traits<compas::Mesh>::null_face() && std::find(vertex_faces.begin(), vertex_faces.end(), f) == vertex_faces.end()) {
                    vertex_faces.push_back(f);

                    auto h = halfedge(f, mesh_a);

                    std::vector<size_t> edges;
                    for(auto h : CGAL::halfedges_around_face(h, mesh_a)) {
                        auto e = edge(h, mesh_a);
                        edges.push_back(e);
                    }
                    faces_edges.push_back(edges);
                }
            }

            // Find the face containing the first boundary edge
            int e1_face = -1;
            for(int i = 0; i < faces_edges.size(); i++) {
                if(std::find(faces_edges[i].begin(), faces_edges[i].end(), e1) != faces_edges[i].end()) {
                    e1_face = i;
                    break;
                }
            }

            /////////////////////////////////////////////////////////////////////////////////
            // Face Ordering Algorithm
            /////////////////////////////////////////////////////////////////////////////////
            // Order the faces around the boundary vertex by starting from the face
            // containing the first boundary edge and traversing through adjacent faces
            // based on shared non-boundary edges.
            
            std::vector<typename boost::graph_traits<compas::Mesh>::face_descriptor> ordered_faces;
            std::vector<bool> visited(vertex_faces.size(), false);

            int current_face_idx = e1_face;
            ordered_faces.push_back(vertex_faces[current_face_idx]);
            visited[current_face_idx] = true;

            while (true) {
                bool found_next = false;
                auto& current_edges = faces_edges[current_face_idx];
                
                for (auto edge1 : current_edges) {
                    if (edge1 == e1 || edge1 == e2) continue; // Skip boundary edges
                    
                    for (size_t j = 0; j < vertex_faces.size(); j++) {
                        if (visited[j]) continue;
                        
                        if (std::find(faces_edges[j].begin(), faces_edges[j].end(), edge1) != faces_edges[j].end()) {
                            ordered_faces.push_back(vertex_faces[j]);
                            visited[j] = true;
                            current_face_idx = j;
                            found_next = true;
                            break;
                        }
                    }
                    
                    if (found_next) break;
                }
                
                if (!found_next) break; // Exit when no more adjacent faces are found
            }

            /////////////////////////////////////////////////////////////////////////////////
            // Boundary Sequence Creation
            /////////////////////////////////////////////////////////////////////////////////
            // Assemble the final boundary sequence by combining the boundary vertex,
            // edge midpoints, and face centroids in the correct order.
            
            std::vector<int> boundary_sequence;
            
            boundary_sequence.push_back(vertex_to_dual_vertex[v]); // Add boundary vertex
            boundary_sequence.push_back(midpoint1_v); // Add first edge midpoint
            
            for(auto f : ordered_faces) {
                if(face_to_vertex.find(f) != face_to_vertex.end()) {
                    boundary_sequence.push_back(face_to_vertex[f]); // Add all face centroids in order
                }
            }
            
            boundary_sequence.push_back(midpoint2_v); // Add second edge midpoint
            
            if(boundary_sequence.size() >= 3) {
                variable_faces.push_back(boundary_sequence); // Only add valid sequences
            }
        }
    }
    
    /////////////////////////////////////////////////////////////////////////////////
    // Final Output Generation
    /////////////////////////////////////////////////////////////////////////////////
    // Convert the dual mesh to matrix format and return the complete dual mesh
    // representation, consisting of vertex coordinates and variable-length faces.
    
    auto [vertices, triangles] = compas::mesh_to_vertices_and_faces(dual_mesh);
    return std::make_tuple(vertices, variable_faces);
}