Function line_triangle_intersection

src/geode/geometry/intersection.cpp:344–413

    }

    IntersectionResult< Point3D > line_triangle_intersection(
        const InfiniteLine3D& line, const Triangle3D& triangle )
    {
        // http://www.geometrictools.com/LibMathematics/Intersection/Intersection.html
        // Compute the offset origin, edges, and normal.
        const auto& vertices = triangle.vertices();
        const Vector3D edge1{ vertices[0], vertices[1] };
        const Vector3D edge2{ vertices[0], vertices[2] };
        const auto normal = edge1.cross( edge2 );

        // Solve Q + t*D = b1*E1 + b2*E2 (Q = diff, D = segment direction,
        // E1 = edge1, E2 = edge2, N = Cross(E1,E2)) by
        //   |Dot(D,N)|*b1 = sign(Dot(D,N))*Dot(D,Cross(Q,E2))
        //   |Dot(D,N)|*b2 = sign(Dot(D,N))*Dot(D,Cross(E1,Q))
        //   |Dot(D,N)|*t = -sign(Dot(D,N))*Dot(Q,N)
        auto d_dot_n = line.direction().dot( normal );
        signed_index_t sign;
        if( d_dot_n > 0. )
        {
            sign = 1;
        }
        else if( d_dot_n < -0. )
        {
            sign = -1;
            d_dot_n = -d_dot_n;
        }
        else
        {
            // Segment and triangle are parallel
            return { INTERSECTION_TYPE::parallel };
        }

        const Vector3D diff{ vertices[0], line.origin() };
        const auto d_dot_q_cross_e2 =
            sign * line.direction().dot( diff.cross( edge2 ) );
        if( d_dot_q_cross_e2 >= -GLOBAL_EPSILON )
        {
            const auto d_dot_e1_cross_q =
                sign * line.direction().dot( edge1.cross( diff ) );
            if( d_dot_e1_cross_q >= -GLOBAL_EPSILON
                && d_dot_q_cross_e2 + d_dot_e1_cross_q <= d_dot_n )
            {
                // InfiniteLine intersects triangle.
                const auto q_dot_n = -sign * diff.dot( normal );
                const auto inv = 1. / d_dot_n;
                const auto seg_parameter = q_dot_n * inv;

                auto result = line.origin() + line.direction() * seg_parameter;
                CorrectnessInfo< Point3D >::Correctness first_correctness{
                    point_line_distance( result, line ) <= GLOBAL_EPSILON,
                    point_line_projection( result, line )
                };
                const auto tri_lambdas =
                    safe_triangle_barycentric_coordinates( result, triangle );
                const auto correctness_point =
                    vertices[0].get() * tri_lambdas[0]
                    + vertices[1].get() * tri_lambdas[1]
                    + vertices[2].get() * tri_lambdas[2];