diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/intersection_3_1_impl.h b/Intersections_3/include/CGAL/Intersections_3/internal/intersection_3_1_impl.h index 8206961896d..a9066d1373c 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/intersection_3_1_impl.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/intersection_3_1_impl.h @@ -223,31 +223,55 @@ intersection_point(const typename K::Plane_3 &plane1, const FT &m00 = plane1.a(); const FT &m01 = plane1.b(); const FT &m02 = plane1.c(); - const FT &b0 = plane1.d(); + const FT &b0 = - plane1.d(); const FT &m10 = plane2.a(); const FT &m11 = plane2.b(); const FT &m12 = plane2.c(); - const FT &b1 = plane2.d(); + const FT &b1 = - plane2.d(); const FT &m20 = plane3.a(); const FT &m21 = plane3.b(); const FT &m22 = plane3.c(); - const FT &b2 = plane3.d(); - - const FT den = - determinant(m00, m01, m02, - m10, m11, m12, - m20, m21, m22); + const FT &b2 = - plane3.d(); + /* + const FT den = determinant(m00, m01, m02, + m10, m11, m12, + m20, m21, m22); if(den == FT(0)){ return boost::none; } const FT num1 = determinant(b0, m01, m02, - b1, m11, m12, - b2, m21, m22); + b1, m11, m12, + b2, m21, m22); const FT num2 = determinant(m00, b0, m02, - m10, b1, m12, - m20, b2, m22); + m10, b1, m12, + m20, b2, m22); const FT num3 = determinant(m00, m01, b0, - m10, m11, b1, - m20, m21, b2); + m10, m11, b1, + m20, m21, b2); + */ + + // Minors common to two determinants + const FT minor_0 = m00*m11 - m10*m01; + const FT minor_1 = m00*m21 - m20*m01; + const FT minor_2 = m10*m21 - m20*m11; + + const FT den = minor_0*m22 - minor_1*m12 + minor_2*m02; // determinant of M + + if(den == FT(0)){ + return boost::none; + } + + const FT num3 = minor_0*b2 - minor_1*b1 + minor_2*b0; // determinant of M with M[x:2] swapped with [b0,b1,b2] + + // Minors common to two determinants + const FT minor_3 = b0*m12 - b1*m02; + const FT minor_4 = b0*m22 - b2*m02; + const FT minor_5 = b1*m22 - b2*m12; + + // num1 has opposite signs because b0 and M[:1] have been swapped + const FT num1 = - minor_3*m21 + minor_4*m11 - minor_5*m01; // determinant of M with M[x:0] swapped with [b0,b1,b2] + const FT num2 = minor_3*m20 - minor_4*m10 + minor_5*m00; // determinant of M with M[x:1] swapped with [b0,b1,b2] + return boost::make_optional(typename K::Point_3(num1/den, num2/den, num3/den)); }