How to use Boost's Kruskal MST algorithm on CGAL's Triangulation_3?

Question

I've been trying to puzzle out how to form edge descriptors for a CGAL Triangulation_3 such that I can use Boost's implementation of Kruskal's Minimum Spanning Tree on that Triangulation.

I have been reading through the reference documentation for a Triangulation_2 (provided here), but have observed that no implementation exists for boost::graph_traits<Triangulation_3>. While puzzling it out, I found that I could possibly provide my own implementation for edge descriptors through an adjacency list as shown in Boost's example for a Kruskal MST, but got lost and confused at this step, and didn't know if that would be a sufficient approach.

Ultimately, it seems that what I need to do is create a boost Graph implementation, but am lost at what resources I need to accomplish this step. From there, the desired use is to be able to traverse this MST to perform graph-based min-cuts at specific edges matching a predicate.

// EDIT :>

My current attempt revolves around creating the EMST via pushing simplex edges defined as a pair of vertex iterate indices, with weights defined as euclidean distance between vertices (a Point_3 with data attached), using the Graph construction shown in the Boost example.

My hope is to have BGL vertices (as a Point_3 with color information attached) be connected by BGL edges (as a simplex[!] edge after the triangulation). My ultimate use just requires that I traverse some sort of contiguous spatial ordering of my Point_3's (with RGB info), and split estimated planes into "patches" which meet a max-distance (complete linkage?) threshold, or a within-patch distance variance. It's not exactly segmentation, but similar.

// some defns...
using RGBA = std::array<uint16_t, 3>;
using PointData = boost::tuple<
    Point_3,    // Point location; Easting-Altitude-Northing
    Vector_3,   // Estimated Normal Vector at Point
    RGBA,       // Photo Color
    RGBA,       // RANSAC Shape Colorization
    size_t,     // Estimated Patch ID
    RGBA>;      // Estimated Patch Colorization

//
//  Some operations on the points and RANSAC estimation occurs here
//

// iterate through shapes
while (it != shapes.end()) {
        boost::shared_ptr<EfficientRANSAC::Shape> shape = *it;
        std::cout << (*it)->info() << std::endl;

        // generate a random color code for this shape
        RGBA rgb;
        for (int i=0; i<3; i++) {
            rgb[i] = rand()%256;
        }

        // Form triangulation to later convert into Graph representation
        using VertexInfoBase = CGAL::Triangulation_vertex_base_with_info_3<
                                    PointData,
                                    Kernel
                                >;
        using TriTraits = CGAL::Triangulation_data_structure_3<
                                VertexInfoBase,
                                CGAL::Delaunay_triangulation_cell_base_3<Kernel>,
                                CGAL::Parallel_tag
                            >;
        using Triangulation_3 = CGAL::Delaunay_triangulation_3<Kernel, TriTraits>;

        Triangulation_3 tr;

        // Iterate through point indices assigned to each detected shape. 
        std::vector<std::size_t>::const_iterator 
            index_it = (*it)->indices_of_assigned_points().begin();

        while (index_it != (*it)->indices_of_assigned_points().end()) {
            PointData& p = *(points.begin() + (*index_it));

            // assign shape diagnostic color info
            boost::get<3>(p) = rgb;

            // insert Point_3 data for triangulation and attach PointData info
            auto vertex = tr.insert(boost::get<0>(p));
            vertex->info() = p;

            index_it++; // next assigned point
        }
        
        std::cout << "Found triangulation with: \n\t" << 
            tr.number_of_vertices() << "\tvertices\n\t" <<
            tr.number_of_edges() << "\tedges\n\t" <<
            tr.number_of_facets() << "\tfacets" << std::endl;

        // build a Graph out of the triangulation that we can do a Minimum-Spanning-Tree on
        using Graph = boost::adjacency_list<
                        boost::vecS,        // OutEdgeList
                        boost::vecS,        // VertexList
                        boost::undirectedS, // Directed
                        boost::no_property, // VertexProperties
                        boost::property< boost::edge_weight_t, int >,  // EdgeProperties
                        boost::no_property, // GraphProperties
                        boost::listS        // EdgeList
                        >;
        using Edge = boost::graph_traits<Graph>::edge_descriptor;
        using E = std::pair< size_t, size_t >; // <: TODO - should be iterator index of vertex in Triangulation_3 instead of size_t?

        std::vector<E> edge_array; // edges should be between Point_3's with attached RGBA photocolor info. 
                                   // It is necessary to later access both the Point_3 and RGBA info for vertices after operations are performed on the EMST
        std::vector<float> weights; // weights are `std::sqrt(CGAL::squared_distance(...))` between these Point_3's

        // Question(?) :> Should be iterating over "finite" edges here? 
        for (auto edge : tr.all_edges()) {
            // insert simplex (!!) edge (between-vertices) here 
            edge_array.push_back(...);
            // generate weight using std::sqrt(CGAL::squared_distance(...))
            weights.push_back(...);
        }

        // build Graph from `edge_array` and `weights`
        Graph g(...);
        // build Euclidean-Minimum-Spanning-Tree (EMST) as list of simplex edges between vertices
        std::list<E> emst;
        boost::kruskal_minimum_spanning_tree(...);

        // - traverse EMST from start of list, performing "cuts" into "patches" when we have hit
        // max patch distance (euclidean) from current "first" vertex of "patch". 
        // - have to be able to access Triangulation_3 vertex info (via `locate`?) here
        // - foreach collection of PointData in patch, assign `patch_id` and diagnostic color info,
        //   then commit individual serialized "patches" collections of Point_3 and RGBA photocolor to database 
        todo!();

        it++; // next shape
    }

The end goal of traversing each of the shapes using a Minimum Spanning Tree via Triangulation is to break each of the RANSAC estimated shapes into chunks, for other purposes. Picture example:

Answer 1

Do you want the graph of vertices and edges, or the graph of the dual, that is the tetrahedra would be BGL vertices and the faces between tetrahedra would be BGL edges?
For both it is not that hard to write the specialization of the graph traits class and the some free functions to navigate. Get inspired by the code for the 2D version for the graph_traits

Answer 2

Ultimately, it seems that what I need to do is create a boost Graph implementation, but am lost at what resources I need to accomplish this step.

The algorithm documents the concept requirements:

You can zoom in on the implications here: VertexListGraph and EdgeListGraph.

I found that I could possibly provide my own implementation for edge descriptors through an adjacency list as shown in Boost's example for a Kruskal MST, but got lost and confused at this step, and didn't know if that would be a sufficient approach.

It would be fine to show your attempt as a question, because it would help us know where you are stuck. Right now there is really no code to "go at", so I'll happily await a newer, more concrete question.

Answer 3

I was able to find an attempt at an answer. I added another property to my Point collection implementation (that included the index of that point in an array), and used this to iterate over edges in the triangulation to build the Graph, before running the EMST algorithm on it.

However, the real answer is don't do this -- it still is not working correctly (incorrect number of edges, including infinite vertices in the edge list, and other problems).

// Form triangulation to later convert into Graph representation
        using VertexInfoBase = CGAL::Triangulation_vertex_base_with_info_3<
                                    PointData,
                                    Kernel
                                >;
        using TriTraits = CGAL::Triangulation_data_structure_3<
                                VertexInfoBase,
                                CGAL::Delaunay_triangulation_cell_base_3<Kernel>,
                                CGAL::Parallel_tag
                            >;
        using Triangulation_3 = CGAL::Delaunay_triangulation_3<Kernel, TriTraits>;

        Triangulation_3 tr;

        // Iterate through point indices assigned to each detected shape. 
        std::vector<std::size_t>::const_iterator 
            index_it = (*it)->indices_of_assigned_points().begin();

        while (index_it != (*it)->indices_of_assigned_points().end()) {
            PointData& p = *(points.begin() + (*index_it));

            // assign shape diagnostic color info
            boost::get<3>(p) = rgb;

            // insert Point_3 data for triangulation and attach PointData info
            TriTraits::Vertex_handle vertex = tr.insert(boost::get<0>(p));
            vertex->info() = p;

            index_it++; // next assigned point
        }
        
        std::cout << "Found triangulation with: \n\t" << 
            tr.number_of_vertices() << "\tvertices\n\t" <<
            tr.number_of_edges() << "\tedges\n\t" <<
            tr.number_of_facets() << "\tfacets" << std::endl;

        // build a Graph out of the triangulation that we can do a Minimum-Spanning-Tree on
        // examples taken from https://www.boost.org/doc/libs/1_80_0/libs/graph/example/kruskal-example.cpp
        using Graph = boost::adjacency_list<
                        boost::vecS,            // OutEdgeList
                        boost::vecS,            // VertexList
                        boost::undirectedS,     // Directed
                        boost::no_property,     // VertexProperties
                        boost::property< boost::edge_weight_t, double >  // EdgeProperties
                        >;
        using Edge = boost::graph_traits<Graph>::edge_descriptor;
        using E = std::pair< size_t, size_t >; // <: TODO - should be iterator index of vertex in Triangulation_3 instead of size_t?

        Graph g(tr.number_of_vertices());
        boost::property_map< Graph, boost::edge_weight_t >::type weightmap = boost::get(boost::edge_weight, g);

        // iterate over finite edges in the triangle, and add these 
        for (
                Triangulation_3::Finite_edges_iterator eit = tr.finite_edges_begin();
                eit != tr.finite_edges_end(); 
                eit++
        ) 
        {
            Triangulation_3::Segment s = tr.segment(*eit);

            Point_3 vtx = s.point(0);
            Point_3 n_vtx = s.point(1);

            // locate the (*eit), get vertex handles?
            // from https://www.appsloveworld.com/cplus/100/204/how-to-get-the-source-and-target-points-from-edge-iterator-in-cgal
            Triangulation_3::Vertex_handle vh1 = eit->first->vertex((eit->second + 1) % 3);
            Triangulation_3::Vertex_handle vh2 = eit->first->vertex((eit->second + 2) % 3);

            double weight = std::sqrt(CGAL::squared_distance(vtx, n_vtx));

            Edge e;
            bool inserted;
            boost::tie(e, inserted)
                = boost::add_edge(
                    boost::get<6>(vh1->info()),
                    boost::get<6>(vh2->info()),
                    g
                );
            weightmap[e] = weight;
        }

        // build Euclidean-Minimum-Spanning-Tree (EMST) as list of simplex edges between vertices
        //boost::property_map<Graph, boost::edge_weight_t>::type weight = boost::get(boost::edge_weight, g);
        std::vector<Edge> spanning_tree;

        boost::kruskal_minimum_spanning_tree(g, std::back_inserter(spanning_tree));

        // we can use something like a hash table to go from source -> target
        // for each of the edges, making traversal easier. 
        // from there, we can keep track or eventually find a source "key" which
        // does not correspond to any target "key" within the table
        std::unordered_map< size_t, std::vector<size_t> > map = {};

        // iterate minimum spanning tree to build unordered_map (hashtable)
        std::cout << "Found minimum spanning tree of " << spanning_tree.size() << " edges for #vertices " << tr.number_of_vertices() << std::endl;
        for (std::vector< Edge >::iterator ei = spanning_tree.begin();
            ei != spanning_tree.end(); ++ei)
        {

            size_t source = boost::source(*ei, g); 
            size_t target = boost::target(*ei, g);
            // << " with weight of " << weightmap[*ei] << std::endl;

            if ( map.find(source) == map.end() ) {
                map.insert(
                    {
                        source, 
                        std::vector({target})
                    }
                );
            } else {
                std::vector<size_t> target_vec = map[source];
                target_vec.push_back(target);
                map[source] = target_vec;
            }
        }

        // iterate over map to find an "origin" node
        size_t origin = 0; 

        for (const auto& it : map) {
            bool exit_flag = false;
            std::vector<size_t> check_targets = it.second;
            for (size_t target : check_targets) {
                if (map.find(target) == map.end()) {
                    origin = target;
                    exit_flag = true;
                    break;
                }
            }
            if (exit_flag) {
                break;
            }
        }

        std::cout << "Found origin of tree with value: " << origin << std::endl;