How to use the output of Polygon_mesh_slicer to control points in Surface_mesh_deformation CGAL?

I'm working on a mesh processing project that requires identifying points that are intersected by a plane and then deforming the mesh using those points. I am using CGAL::Polygon_mesh_slicer<> to find the points of interest and want to use those points as control points for the deformation in CGAL::Surface_mesh_deformation<>.

I have been scouring the documentation and forums to find something to get me started. In the documentation it shows that the iterator past to the Polygon_mesh_slicer to get the polylines must be of the type std::vector< traits::Point_3> reference manual, while the input for vertices in the Surface_mesh_deformation must be a boost::graph_traits< Triangle_mesh >::vertex_descriptor from the reference manual
I don't know how to get relate the information from Point_3 to vertex_descriptor.

This is a function that I threw together to illustrate the problem. I have run it using a Polyhedron mesh with a typedef of:

CGAL::Polyhedron_3< CGAL::Exact_predicates_inexact_constructions_kernel >

and a plane using:

CGAL::Exact_predicates_inexact_constructions_kernel::Plane_3

I would really appreciate any help thank you!

 template  <typename cgalMesh, typename cgalPlane, typename K>
 void functionForSOforum(cgalMesh& Mesh, cgalPlane& plane, K kernel)
 {
     typedef typename K::Point_3                                                     Point_3;
     typedef typename CGAL::Polyhedron_3<K>                                       Polyhedron;
     typedef typename K::Plane_3                                                    Plane;
     typedef typename K::Triangle_3                                              Triangle;

     typedef typename boost::graph_traits<cgalMesh>::vertex_descriptor     vertex_descriptor;
     typedef typename boost::graph_traits<cgalMesh>::vertex_iterator         vertex_iterator;
     typedef typename boost::graph_traits<cgalMesh>::halfedge_descriptor halfedge_descriptor;
     typedef typename boost::graph_traits<cgalMesh>::halfedge_iterator     halfedge_iterator;

     // Define the maps
     typedef std::map<vertex_descriptor, std::size_t>                   Vertex_id_map;
     typedef std::map<halfedge_descriptor, std::size_t>                  Hedge_id_map;
     typedef boost::associative_property_map<Vertex_id_map>            Vertex_id_pmap;
     typedef boost::associative_property_map<Hedge_id_map>              Hedge_id_pmap;

     typedef std::vector<Point_3>                                       Polyline_type;
     typedef std::list< Polyline_type>                                      Polylines;

     typedef CGAL::Polygon_mesh_slicer<Polyhedron, K>   Slicer;

     typedef CGAL::Surface_mesh_deformation<Polyhedron, Vertex_id_pmap, Hedge_id_pmap> Surface_mesh_deformation;

     // Get the least square plane describing the mesh
     linear_least_squares_fitting_3(Mesh.points_begin(),Mesh.points_end(), plane, CGAL::Dimension_tag<0>());

     // Find the intersection of the plane and the mesh
     Slicer slicer(Mesh);

      Polylines polylines;
      slicer(plane, std::back_inserter(polylines));

      // Set up the Region of interest as the whole mesh
      // *** from the CGAL docs ***
      Vertex_id_map vertex_index_map;
      vertex_iterator vb, ve;
      std::size_t counter = 0;
      for(boost::tie(vb, ve) = vertices(Mesh); vb != ve; ++vb, ++counter)
        vertex_index_map[*vb]=counter;

      Hedge_id_map hedge_index_map;
      counter = 0;
      halfedge_iterator eb, ee;
      for(boost::tie(eb, ee) = halfedges(Mesh); eb != ee; ++eb, ++counter)
        hedge_index_map[*eb]=counter;

      Surface_mesh_deformation deform_mesh( Mesh,
                                            Vertex_id_pmap(vertex_index_map),
                                            Hedge_id_pmap(hedge_index_map) );
      // Insert the whole mesh as region of interest
        boost::tie(vb, ve) = vertices(Mesh);
        deform_mesh.insert_roi_vertices(vb, ve);

      // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      // ~~~~~     Want to put the points from the slicer here    ~~~~~~~~~

        int randomPoint = 200;

        // Insert a control vertex using the
         deform_mesh.insert_control_vertex(*std::next(vb, randomPoint));

 }

Solution

The first step is to refine your input mesh with the intersection of the mesh with your plane. This can be done using the AABB-tree package. More specifically, consider the edges of your mesh as primitives by using the class AABB_halfedge_graph_segment_primitive, like in this example (note that you can replace Polyhedron_3 by Surface_mesh). Then calling the all_intersections() method from the AABB-tree will give you all the intersection points and edges to be refined.

To split an edge you can use the function split_edge(). Note that you have to keep the vertex descriptors returned by the function for the deformation process.

Then a call to triangulate_faces() will restore the triangulated property of your mesh. (You can probably use split_face instead but this require some house keeping and sort the intersection points per edge first).

Once this is done, you can use the vertex descriptors previously saved in the deformation algorithm.

Note that the intersection edges won't be present in the mesh. If you need them to be then extra steps are needed (or directly use the function corefine with a planar mesh.