Correct usage of CGAL::Euler::split_edge()

[nimafo]

Hi,
I want to split my polyhedral surface edges at a given point.
I use P.split_edge(edge_to_split, splitPoint); but it returns error on the second argument (too many arguments in function call).
How can I correctly add a split a given halfedge at a given point?

This is the full code:

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <iostream>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3 Point_3;
typedef CGAL::Polyhedron_3<K> Polyhedron;

int main() {
    Polyhedron P;

    // Add vertices by creating faces
    Polyhedron::Halfedge_handle he0 = P.make_tetrahedron(Point_3(0, 0, 0), Point_3(1, 0, 0), Point_3(1, 1, 0), Point_3(0, 1, 0));

    // Split an edge at a specified point
    Point_3 splitPoint(0.5, 0.5, 0);

    // Find the edge to split (for simplicity, we choose an edge)
    Polyhedron::Halfedge_handle edge_to_split = he0;

    // Split the edge
    Polyhedron::Halfedge_handle new_halfedge = P.split_edge(edge_to_split, splitPoint);
    std::cout << typeid(new_halfedge).name() << std::endl;

    // Output the modified mesh
    std::cout << "Number of vertices: " << P.size_of_vertices() << std::endl;
    std::cout << "Number of half edges: " << P.size_of_halfedges() / 2 << std::endl;
    std::cout << "Number of faces: " << P.size_of_facets() << std::endl;

    return 0;
}

[sloriot]

Should be something like:

Polyhedron::Halfedge_handle new_halfedge = P.split_edge(edge_to_split);
new_halfedge->target()->point() = splitPoint;

[nimafo]

Thanks @sloriot,
But new_halfedge does not have target() :/

[sloriot]

vertex() sorry

[nimafo]

Thank you for your response @sloriot,
However, I still seem to be implementing it wrong,
What I am looking for is this:

Yet after intersecting, I get 16 half-edges instead of 12, and the resulting halfedges are not what I expect them to be.

What I am trying to do is intersect a line with a polyhedral surface (a quad here), and make two separate polyhedral surfaces as a result (if the line does intersect the surface).
If there are simpler ways to do this, please let me know.

This is my complete code so far:

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_incremental_builder_3.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/intersections.h> //temp
#include <cassert>
#include <vector> //temp
#include <set>

template <class HDS>
class Build_quad : public CGAL::Modifier_base<HDS>
{
public:
    Build_quad() {}

    void operator()(HDS &hds)
    {
        CGAL::Polyhedron_incremental_builder_3<HDS> B(hds, true);
        B.begin_surface(4, 1, 8);
        typedef typename HDS::Vertex Vertex;
        typedef typename Vertex::Point Point;
        B.add_vertex(Point(0, 0, 0));
        B.add_vertex(Point(1, 0, 0));
        B.add_vertex(Point(1, 1, 0));
        B.add_vertex(Point(0, 1, 0));
        B.begin_facet();
        B.add_vertex_to_facet(0);
        B.add_vertex_to_facet(1);
        B.add_vertex_to_facet(2);
        B.add_vertex_to_facet(3);
        B.end_facet();
        B.end_surface();
    }
};

typedef CGAL::Simple_cartesian<double> Kernel;
typedef Kernel::Point_3 Point_3;
typedef Kernel::Segment_3 Segment_3;
typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
typedef Polyhedron::HalfedgeDS HalfedgeDS;
typedef Kernel::Line_3 Line_3;
typedef Kernel::Vector_3 Vector_3;

// function declarations
void print_halfedgeInformation(const Polyhedron &P);

int main()
{
    const double epsilon = 1e-6;
    Polyhedron P;
    Build_quad<HalfedgeDS> quad;
    P.delegate(quad);
    // VALIDATION !
    print_halfedgeInformation(P);
    // --------------------------------------------------------------------------------------------------------------------------------------------------------------

    // Define the line given two points
    Point_3 p(0, 0.5, 0);
    Point_3 q(1, 0.5, 0);
    Line_3 line(p, q);
    // ------------------------------------------------------------------------ Find intersection points ------------------------------------------------------------
    int counter = 0;
    std::set<Point_3> Intersection_PointSet;
    for (Polyhedron::Halfedge_iterator he = P.halfedges_begin(); he != P.halfedges_end(); ++he)
    {
        if (he->is_border())
        {
            continue;
        }
        std::cout << " ------------------------ " << counter << " ------------------------ " << std::endl;
        counter++;
        Polyhedron::Vertex_handle v1 = he->vertex();
        Polyhedron::Vertex_handle v2 = he->next()->vertex();
        std::cout << "Half-edge: (" << v1->point() << ") -> (" << v2->point() << ")" << std::endl;
        Segment_3 segment(v1->point(), v2->point());
        if (v1->point() != v2->point())
        {
            auto result = CGAL::do_intersect(line, segment); // Checks whether line and segment intersect.
            if (result)
            {
                std::cout << "Does intersect!" << std::endl;
                // Find the intersection point
                auto intersection_result = CGAL::intersection(line, segment);
                if (const Point_3 *p = boost::get<Point_3>(&*intersection_result))
                {
                    std::cout << "Intersection point: " << *p << std::endl;
                    Intersection_PointSet.insert(*p);
                }
                else
                {
                    std::cout << "Intersection result is not a point!" << std::endl;
                }
            }
            else
            {
                std::cout << " Intersection NOT found !!" << std::endl;
            }
        }
        else
        {
            std::cout << "two vertices of the segment are the same! " << std::endl;
        }
    }

    // ------------------------------------------------------------------------ Intersect points with half edges  ------------------------------------------------------------

    std::cout << std::endl
              << std::endl;
    std::cout << "Number of intersection pointsSet are: " << Intersection_PointSet.size() << std::endl
              << std::endl
              << std::endl;

    // Intersect points in the "Intersection_PointSet" with all the half-edges in a copy of P

    Polyhedron P_(P);
    for (Polyhedron::Halfedge_iterator he = P_.halfedges_begin(); he != P_.halfedges_end(); ++he)
    {
        for (const Point_3 &element : Intersection_PointSet)
        {
            // Check to see if the halfedge has the point on it:
            Point_3 endpoint1 = he->vertex()->point();
            Point_3 endpoint2 = he->opposite()->vertex()->point();
            Segment_3 segment(endpoint1, endpoint2);

            if (segment.has_on(element))
            {
                std::cout << " Element " << element << " is on segment." << endpoint1 << " | " << endpoint2 << std::endl;
                // Intersect Point(element) with halfedge(he) and make new half edge? where does the half edge go?
                Polyhedron::Halfedge_handle new_halfedge = P.split_edge(he);
                // new_halfedge->vertex()->point() = element;
                // std::cout << "Vertex of the splitted half-edge: " << new_halfedge->vertex()->point() << std::endl;
            }
            else
            {
                // std::cout << "segment does not have element on it! " << std::endl;
                continue;
            }
        }
    }
    // VALIDATION !
    print_halfedgeInformation(P);
    return 0;
}

// VALIDATION !
// Function implementation
void print_halfedgeInformation(const Polyhedron &P)
{
    std::cout << "---------------------------------------------------" << std::endl;
    std::cout << "Is it empty?  " << P.is_empty() << std::endl;
    std::cout << "Is it closed?  " << P.is_closed() << std::endl;
    std::cout << "Is it pure bivalent?  " << P.is_pure_bivalent() << std::endl;
    std::cout << "Is it pure quad?  " << P.is_pure_quad() << std::endl;
    std::cout << "Is it pure triangle?  " << P.is_pure_triangle() << std::endl;
    std::cout << "Is it valid?  " << P.is_valid() << std::endl;
    std::cout << "Size of half edges?  " << P.size_of_halfedges() << std::endl;
    std::cout << "---------------------------------------------------" << std::endl;
    std::cout << "---------------------------------------------------" << std::endl;
    int counter = 0;
    for (Polyhedron::Halfedge_const_iterator he = P.halfedges_begin(); he != P.halfedges_end(); ++he)
    {
        std::cout << "---------------- Halfedge number: " << counter << "----------------" << std::endl;
        std::cout << he->opposite()->vertex()->point() << " TO " << he->vertex()->point() << std::endl;
        counter++;
    }
    std::cout << "---------------------------------------------------" << std::endl;
}

And this is the output of the program:

---------------------------------------------------
Is it empty?  0
Is it closed?  0
Is it pure bivalent?  1
Is it pure quad?  1
Is it pure triangle?  0
Is it valid?  1
Size of half edges?  8
---------------------------------------------------
---------------------------------------------------
---------------- Halfedge number: 0----------------
0 0 0 TO 1 0 0
---------------- Halfedge number: 1----------------
1 0 0 TO 0 0 0
---------------- Halfedge number: 2----------------
1 0 0 TO 1 1 0
---------------- Halfedge number: 3----------------
1 1 0 TO 1 0 0
---------------- Halfedge number: 4----------------
1 1 0 TO 0 1 0
---------------- Halfedge number: 5----------------
0 1 0 TO 1 1 0
---------------- Halfedge number: 6----------------
0 1 0 TO 0 0 0
---------------- Halfedge number: 7----------------
0 0 0 TO 0 1 0
---------------------------------------------------
 ------------------------ 0 ------------------------ 
Half-edge: (1 0 0) -> (1 1 0)
Does intersect!
Intersection point: 1 0.5 0
 ------------------------ 1 ------------------------ 
Half-edge: (1 1 0) -> (0 1 0)
 Intersection NOT found !!!!!!!!!!!!!!!!!
 ------------------------ 2 ------------------------ 
Half-edge: (0 1 0) -> (0 0 0)
Does intersect!
Intersection point: 0 0.5 0
 ------------------------ 3 ------------------------ 
Half-edge: (0 0 0) -> (1 0 0)
 Intersection NOT found !!!!!!!!!!!!!!!!!


Number of intersection pointsSet are: 2


 Element 1 0.5 0 is on segment.1 1 0 | 1 0 0
 Element 1 0.5 0 is on segment.1 0 0 | 1 1 0
 Element 0 0.5 0 is on segment.0 0 0 | 0 1 0
 Element 0 0.5 0 is on segment.0 1 0 | 0 0 0
---------------------------------------------------
Is it empty?  0
Is it closed?  0
Is it pure bivalent?  1
Is it pure quad?  1
Is it pure triangle?  0
Is it valid?  0
Size of half edges?  16
---------------------------------------------------
---------------------------------------------------
---------------- Halfedge number: 0----------------
0 0 0 TO 1 0 0
---------------- Halfedge number: 1----------------
1 0 0 TO 0 0 0
---------------- Halfedge number: 2----------------
1 0 0 TO 1 1 0
---------------- Halfedge number: 3----------------
1 1 0 TO 1 0 0
---------------- Halfedge number: 4----------------
1 1 0 TO 0 1 0
---------------- Halfedge number: 5----------------
0 1 0 TO 1 1 0
---------------- Halfedge number: 6----------------
0 1 0 TO 0 0 0
---------------- Halfedge number: 7----------------
0 0 0 TO 0 1 0
---------------- Halfedge number: 8----------------
1 0 0 TO 1 0 0
---------------- Halfedge number: 9----------------
1 0 0 TO 1 0 0
---------------- Halfedge number: 10----------------
1 1 0 TO 1 1 0
---------------- Halfedge number: 11----------------
1 1 0 TO 1 1 0
---------------- Halfedge number: 12----------------
0 1 0 TO 0 1 0
---------------- Halfedge number: 13----------------
0 1 0 TO 0 1 0
---------------- Halfedge number: 14----------------
0 0 0 TO 0 0 0
---------------- Halfedge number: 15----------------
0 0 0 TO 0 0 0
---------------------------------------------------

[sloriot]

vertex() is the target of the edge so you are not splitting the intersected halfedge.Also you cannot iterate over halfedge and split halfedge in the loop, this does invalidate the iterator.

[nimafo]

@sloriot Thanks a lot for your quick response,

I see,
So back to my original question, how can I split an edge in a polyhedral surface like the above figure? is split_edge() the right way to do so?
Can you please give me an example of the proper usage?

PS: This is the whole pipeline I want to have eventually:

[sloriot]

you're almost there. You need to pick h->vertex() and h->opposite()->vertex() for h endpoints. Then collect the split points and the halfedge in a separate container and do the split iterating over that container and not the mesh halfedge container.

[nimafo]

@sloriot Thanks again for your quick response,

After trying a few times I could not really understand and implement what you said here. It is most probably because of my lack of sufficient experience with C++ and CGAL.

Can you please explain that in a clearer way?
Thanks in advance!

[sloriot]

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_incremental_builder_3.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/intersections.h> //temp
#include <cassert>
#include <vector> //temp
#include <set>

template <class HDS>
class Build_quad : public CGAL::Modifier_base<HDS>
{
public:
    Build_quad() {}

    void operator()(HDS &hds)
    {
        CGAL::Polyhedron_incremental_builder_3<HDS> B(hds, true);
        B.begin_surface(4, 1, 8);
        typedef typename HDS::Vertex Vertex;
        typedef typename Vertex::Point Point;
        B.add_vertex(Point(0, 0, 0));
        B.add_vertex(Point(1, 0, 0));
        B.add_vertex(Point(1, 1, 0));
        B.add_vertex(Point(0, 1, 0));
        B.begin_facet();
        B.add_vertex_to_facet(0);
        B.add_vertex_to_facet(1);
        B.add_vertex_to_facet(2);
        B.add_vertex_to_facet(3);
        B.end_facet();
        B.end_surface();
    }
};

typedef CGAL::Simple_cartesian<double> Kernel;
typedef Kernel::Point_3 Point_3;
typedef Kernel::Segment_3 Segment_3;
typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
typedef Polyhedron::HalfedgeDS HalfedgeDS;
typedef Kernel::Line_3 Line_3;
typedef Kernel::Vector_3 Vector_3;

// function declarations
void print_halfedgeInformation(const Polyhedron &P);

int main()
{
    const double epsilon = 1e-6;
    Polyhedron P;
    Build_quad<HalfedgeDS> quad;
    P.delegate(quad);
    // VALIDATION !
    print_halfedgeInformation(P);
    // --------------------------------------------------------------------------------------------------------------------------------------------------------------

    // Define the line given two points
    Point_3 p(0, 0.5, 0);
    Point_3 q(1, 0.5, 0);
    Line_3 line(p, q);
    // ------------------------------------------------------------------------ Find intersection points ------------------------------------------------------------
    int counter = 0;
    std::vector<std::pair<Point_3, Polyhedron::Halfedge_handle>> Intersection_PointSet;
    for (Polyhedron::Edge_iterator eit = P.edges_begin(); eit != P.edges_end(); ++eit)
    {
        std::cout << " ------------------------ " << counter << " ------------------------ " << std::endl;
        counter++;
        Polyhedron::Halfedge_handle he = eit;
        Polyhedron::Vertex_handle v1 = he->opposite()->vertex();
        Polyhedron::Vertex_handle v2 = he->vertex();
        std::cout << "Half-edge: (" << v1->point() << ") -> (" << v2->point() << ")" << std::endl;
        Segment_3 segment(v1->point(), v2->point());
        if (v1->point() != v2->point())
        {
            auto result = CGAL::do_intersect(line, segment); // Checks whether line and segment intersect.
            if (result)
            {
                std::cout << "Does intersect!" << std::endl;
                // Find the intersection point
                auto intersection_result = CGAL::intersection(line, segment);
                if (const Point_3 *p = boost::get<Point_3>(&*intersection_result))
                {
                    std::cout << "Intersection point: " << *p << std::endl;
                    Intersection_PointSet.emplace_back(*p, he);
                }
                else
                {
                    std::cout << "Intersection result is not a point!" << std::endl;
                }
            }
            else
            {
                std::cout << " Intersection NOT found !!" << std::endl;
            }
        }
        else
        {
            std::cout << "two vertices of the segment are the same! " << std::endl;
        }
    }

    // ------------------------------------------------------------------------ Intersect points with half edges  ------------------------------------------------------------

    std::cout << std::endl
              << std::endl;
    std::cout << "Number of intersection pointsSet are: " << Intersection_PointSet.size() << std::endl
              << std::endl
              << std::endl;

    // Intersect points in the "Intersection_PointSet" with all the half-edges in a copy of P
    for (const auto& p : Intersection_PointSet)
      P.split_edge(p.second)->vertex()->point()=p.first;

    // VALIDATION !
    print_halfedgeInformation(P);
    return 0;
}

// VALIDATION !
// Function implementation
void print_halfedgeInformation(const Polyhedron &P)
{
    std::cout << "---------------------------------------------------" << std::endl;
    std::cout << "Is it empty?  " << P.is_empty() << std::endl;
    std::cout << "Is it closed?  " << P.is_closed() << std::endl;
    std::cout << "Is it pure bivalent?  " << P.is_pure_bivalent() << std::endl;
    std::cout << "Is it pure quad?  " << P.is_pure_quad() << std::endl;
    std::cout << "Is it pure triangle?  " << P.is_pure_triangle() << std::endl;
    std::cout << "Is it valid?  " << P.is_valid() << std::endl;
    std::cout << "Size of half edges?  " << P.size_of_halfedges() << std::endl;
    std::cout << "---------------------------------------------------" << std::endl;
    std::cout << "---------------------------------------------------" << std::endl;
    int counter = 0;
    for (Polyhedron::Halfedge_const_iterator he = P.halfedges_begin(); he != P.halfedges_end(); ++he)
    {
        std::cout << "---------------- Halfedge number: " << counter << "----------------" << std::endl;
        std::cout << he->opposite()->vertex()->point() << " TO " << he->vertex()->point() << std::endl;
        counter++;
    }
    std::cout << "---------------------------------------------------" << std::endl;
}

Here is the updated code. Note that I got rid of P_. If you really need to keep P unchanged, then make a copy of P before looking for intersections.

[nimafo]

It is exactly what I was looking for. Thanks a lot, @sloriot!