Wrong intersection point when computing the intersection of two 3D triangles

[ImperatorS79]

Issue Details

I am currently using CGAL::box_self_intersection_d to test if a surface made of triangles self-intersect. It might not be the best tool from CGAL to do the job, but it should work nonetheless. Normally, the only intersections should be points and segments of the triangles, but I obtain the following for one of the intersection:

Intersection point: 1.5400e-02 0.0000e+00 1.5400e-02
Triangle A: 8.0606e-03 1.7500e-01 8.0606e-03 ; 1.4600e-02 1.7500e-01 0.0000e+00 ; 1.5400e-02 1.7500e-01 1.5400e-02
Triangle B: 1.4600e-02 1.7500e-01 0.0000e+00 ; 2.4625e-02 1.6399e-01 0.0000e+00 ; 9.9918e-03 1.6476e-01 0.0000e+00

There is clearly a problem, because the point is absolutely not correct: the intersection should be 1.4600e-02 1.7500e-01 0.0000e+00.

In debug mode, the following precondition in the CGAL code is not satisfied (file Segment_3_Segment_3_intersection.h):

CGAL_precondition(!s1.is_degenerate() && !s2.is_degenerate());

The problem arise with Exact_predicates_inexact_constructions_kernel, with Exact_predicates_exact_constructions_kernel I have no problems (apart from the fact it is sloooow ^^').

Source Code

EDIT: Here is an example of source code which does not return the correct answer:

#include <CGAL/box_intersection_d.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <iostream>
#include <vector>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Box_intersection_d::Box_with_handle_d<double, 3, std::vector<K::Triangle_3>::iterator> Box3D;

int main()
{
    struct Report
    {
        std::vector<K::Triangle_3>* m_pTriangles;
        bool* m_pIntersect;

        Report(std::vector<K::Triangle_3>& triangles, bool& intersect):
        m_pTriangles(&triangles), m_pIntersect(&intersect)
        {

        }

        void operator()(const Box3D& a, const Box3D& b) const
        {
            if(!a.handle()->is_degenerate() && !b.handle()->is_degenerate())
            {
                K::Triangle_3 triA = *(a.handle());
                K::Triangle_3 triB = *(b.handle());

                if(CGAL::do_intersect(triA, triB))
                {
                    std::cout << "Triangles which intersect: " << triA << " vs " << triB << std::endl;
                    auto intersection = CGAL::intersection(triA, triB);

                    if(intersection.has_value())
                    {
                        bool ok = false;
                        if(const K::Point_3* pPoint = std::get_if<K::Point_3>(&*intersection))
                            std::cout << "Point: " << *pPoint << std::endl;
                        else if(const K::Segment_3* pSegment = std::get_if<K::Segment_3>(&*intersection))
                            std::cout << "Segment: " << *pSegment << std::endl;
                        else if(const K::Triangle_3* pTriangle = std::get_if<K::Triangle_3>(&*intersection))
                            std::cout << "Triangle: " << *pTriangle << std::endl;
                        else if(const std::vector<K::Point_3>* pPointVec = std::get_if<std::vector<K::Point_3>>(&*intersection))
                        {
                            std::cout << "Polygon: " << std::endl;
                            for(std::size_t i = 0 ; i < pPointVec->size() ; ++i)
                                std::cout << pPointVec->at(i) << std::endl; 
                        }
                        else
                            std::cout << "Not handled currently" << std::endl;
                    }
                    else
                        std::cout << "Intersect but no intersection ?" << std::endl;

                    *m_pIntersect = true;
                }    
            }
        }
    };

    std::vector<K::Triangle_3> triangles(2);

    K::Point_3 pA = K::Point_3(0.0000000000000000e+00, 1.0208333333314730e-01, 2.9199999999999998e-01);
    K::Point_3 pB = K::Point_3(0.0000000000000000e+00, 8.7499999999775646e-02, 2.9199999999999998e-01);
    K::Point_3 pC = K::Point_3(1.0617182181505985e-02, 8.7299626865148805e-02, 2.9199898847343825e-01);

    K::Point_3 pD = K::Point_3(0.0000000000000000e+00, 7.2916666666482610e-02, 2.9199999999999998e-01);
    K::Point_3 pE = K::Point_3(0.0000000000000000e+00, 8.7499999999775646e-02, 2.9199999999999998e-01);
    K::Point_3 pF = K::Point_3(1.0617182181505985e-02, 8.7299626865148805e-02, 2.9199898847343825e-01);

    triangles[0] = K::Triangle_3(pA, pB, pC);
    triangles[1] = K::Triangle_3(pD, pE, pF);

    std::vector<Box3D> boxes;
    for(auto i = triangles.begin(); i != triangles.end(); ++i)
        boxes.push_back(Box3D(i->bbox(), i));
    
    // Run the self intersection algorithm with all defaults
    bool intersect = false;
    CGAL::box_self_intersection_d(boxes.begin(), boxes.end(), Report(triangles, intersect));
    return static_cast<int>(intersect);
}

The intersection should be the segment made of point B and C, but I obtain the point (0, 0.0874995, -0) which has nothing to do with that segment.

Environment

• Operating system : Ubuntu 24.04 LTS 64 bits
• Compiler: gcc 13.3.0
• Release or debug mode: all
• Specific flags used (if any): NA
• CGAL version: 6.1
• Boost version: 1.83.0.1
• Other libraries versions if used (Eigen, TBB, etc.): I use a lot of them, but I am not sure it is relevant here.

[MaelRL]

If you use a kernel providing inexact constructions, it is not unexpected that a construction would return a wrong result.

But if you are interested in whether a surface self intersects, you do not need to actually construct anything, you can use the do_intersect functions, which are predicates and will thus be correct even with CGAL::Exact_predicates_inexact_constructions_kernel.

Also the following functions can be used for this test:

• PMP::does_self_intersect() if you have a triangle mesh.
• PMP::does_triangle_soup_self_intersect() if you have a triangle soup.

[ImperatorS79]

I actually need to construct the intersection in order to check for "dummy" intersection of the triangles.

Even with inexact constructions, I would not expect:

1. That an assertion would be triggered in debug, meaning something in the intersection algorithm is not working correctly.
2. That the intersection computed would be "so far" the actual intersection point (which is here simply one of the triangle vertex)

I will try the functions from PMP, but this does not change this issue nonetheless. This is a bug.

EDIT: The PMP::does_self_intersect() function depends of Intersecting Sequences of dD Iso-oriented Boxes, so this would have the same problem.

[MaelRL]

To investigate what's the issue with this specific assertion trigger, we will need the code & input data, but if you have written something that assumes exact constructions while using a kernel not providing this, anything can happen.

If the construction is inexact, it might very well be wildly wrong: the code might be taking a "wrong" branch, or doing a division by something close to 0 for which an epsilon error can result in a large error, etc. That is the whole point of using certified predicates in which you control the error and take the correct decisions. There is more information on https://www.cgal.org/exact.html and https://doc.cgal.org/latest/Manual/index.html.

The PMP functions are written using predicates, and combinatorial checks to avoid returning a false positive for neighboring elements, these functions do not actually compute any intersection.

[ImperatorS79]

I added source code demonstrating my issue. The fact that I obtain an intersection point which has nothing to do with the intersection segment is bothering me, even with inexact construction.

[ImperatorS79]

To investigate what's the issue with this specific assertion trigger, we will need the code & input data, but if you have written something that assumes exact constructions while using a kernel not providing this, anything can happen.

If the construction is inexact, it might very well be wildly wrong: the code might be taking a "wrong" branch, or doing a division by something close to 0 for which an epsilon error can result in a large error, etc. That is the whole point of using certified predicates in which you control the error and take the correct decisions. There is more information on https://www.cgal.org/exact.html and https://doc.cgal.org/latest/Manual/index.html.

The PMP functions are written using predicates, and combinatorial checks to avoid returning a false positive for neighboring elements, these functions do not actually compute any intersection.

An additional problem I have is that my surface mesh is non-manifold (i.e. and edge might have more than two faces). This prevents me from using does_self_intersect (And I just discovered that the face of the surface mesh must be provided with a specific orientation to be added in the surface mesh...), or isotropic_remeshing.

Is not there a data structure in CGAL which could handle such kind of surface mesh (with potentially more than two faces per edge ?) in order to check if the mesh self-intersect and to remesh the surface using isotropic_remeshing ? I don't see anything which would prevent an algorithm from checking if such a meshs self-intersect, or to remesh it.

[MaelRL]

An additional problem I have is that my surface mesh is non-manifold (i.e. and edge might have more than two faces). This prevents me from using does_self_intersect (And I just discovered that the face of the surface mesh must be provided with a specific orientation to be added in the surface mesh...), or isotropic_remeshing.

For self-intersection detection only, you can use a polygon soup representation, here are some functions that might be useful:

• CGAL::IO::read_polygon_soup()
• CGAL::Polygon_mesh_processing::polygon_mesh_to_polygon_soup()
• CGAL::Polygon_mesh_processing::does_triangle_soup_self_intersect()

Is not there a data structure in CGAL which could handle such kind of surface mesh (with potentially more than two faces per edge ?) in order to check if the mesh self-intersect and to remesh the surface using isotropic_remeshing ? I don't see anything which would prevent an algorithm from checking if such a meshs self-intersect, or to remesh it.

You can obtain a valid combinatorial mesh structure by reading a polygon soup and using orientation functions, which will duplicate your non manifold elements, see the following example: https://doc.cgal.org/latest/Polygon_mesh_processing/Polygon_mesh_processing_2orient_polygon_soup_example_8cpp-example.html

Note however that you'll have a valid mesh structure, but CGAL::Polygon_mesh_processing::does_self_intersect() will consider that there is a self-intersection at the non-manifold elements.

You can nevertheless call isotropic remeshing on this.