mirror of https://github.com/CGAL/cgal
rewrite predicate
This commit is contained in:
Andreas Fabri
parent
2d0dfd2822
commit
f747db2774
|
|
@ -1,4 +1,4 @@
|
||||||
// Copyright (c) 20025 GeometryFactory (France).
|
// Copyright (c) 2025 GeometryFactory (France).
|
||||||
// All rights reserved.
|
// All rights reserved.
|
||||||
//
|
//
|
||||||
// This file is part of CGAL (www.cgal.org)
|
// This file is part of CGAL (www.cgal.org)
|
||||||
|
|
|
||||||
|
|
@ -1,276 +1,109 @@
|
||||||
|
// Copyright (c) 2025 GeometryFactory (France).
|
||||||
|
// All rights reserved.
|
||||||
|
//
|
||||||
|
// This file is part of CGAL (www.cgal.org)
|
||||||
|
//
|
||||||
|
// $URL$
|
||||||
|
// $Id$
|
||||||
|
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
|
||||||
|
//
|
||||||
|
//
|
||||||
|
// Author(s) : Andreas Fabri
|
||||||
|
|
||||||
inline double determinant( double a00, double a01, double a02, double a03, double a04, double a10, double a11, double a12, double a13, double a14, double a20, double a21, double a22, double a23, double a24, double a30, double a31, double a32, double a33, double a34, double a40, double a41, double a42, double a43, double a44) {
|
typedef double RT;
|
||||||
double m01;
|
|
||||||
m01 = ((a10 * a01) - (a00 * a11));
|
RT
|
||||||
double m02;
|
determinant_with_1_in_row_0(
|
||||||
m02 = ((a20 * a01) - (a00 * a21));
|
RT a01, RT a02, RT a03, RT a04, RT a05,
|
||||||
double m03;
|
RT a11, RT a12, RT a13, RT a14, RT a15,
|
||||||
m03 = ((a30 * a01) - (a00 * a31));
|
RT a21, RT a22, RT a23, RT a24, RT a25,
|
||||||
double m04;
|
RT a31, RT a32, RT a33, RT a34, RT a35,
|
||||||
m04 = ((a40 * a01) - (a00 * a41));
|
RT a41, RT a42, RT a43, RT a44, RT a45,
|
||||||
double m12;
|
RT a51, RT a52, RT a53, RT a54, RT a55)
|
||||||
m12 = ((a20 * a11) - (a10 * a21));
|
{
|
||||||
double m13;
|
// First compute the det2x2
|
||||||
m13 = ((a30 * a11) - (a10 * a31));
|
const RT m01 = a11 - a01;
|
||||||
double m14;
|
const RT m02 = a21 - a01;
|
||||||
m14 = ((a40 * a11) - (a10 * a41));
|
const RT m03 = a31 - a01;
|
||||||
double m23;
|
const RT m04 = a41 - a01;
|
||||||
m23 = ((a30 * a21) - (a20 * a31));
|
const RT m05 = a51 - a01;
|
||||||
double m24;
|
const RT m12 = a21 - a11;
|
||||||
m24 = ((a40 * a21) - (a20 * a41));
|
const RT m13 = a31 - a11;
|
||||||
double m34;
|
const RT m14 = a41 - a11;
|
||||||
m34 = ((a40 * a31) - (a30 * a41));
|
const RT m15 = a51 - a11;
|
||||||
double m012;
|
const RT m23 = a31 - a21;
|
||||||
m012 = (((m12 * a02) - (m02 * a12)) + (m01 * a22));
|
const RT m24 = a41 - a21;
|
||||||
double m013;
|
const RT m25 = a51 - a21;
|
||||||
m013 = (((m13 * a02) - (m03 * a12)) + (m01 * a32));
|
const RT m34 = a41 - a31;
|
||||||
double m014;
|
const RT m35 = a51 - a31;
|
||||||
m014 = (((m14 * a02) - (m04 * a12)) + (m01 * a42));
|
const RT m45 = a51 - a41;
|
||||||
double m023;
|
// Now compute the minors of rank 3
|
||||||
m023 = (((m23 * a02) - (m03 * a22)) + (m02 * a32));
|
const RT m012 = m01*a22 - m02*a12 + m12*a02;
|
||||||
double m024;
|
const RT m013 = m01*a32 - m03*a12 + m13*a02;
|
||||||
m024 = (((m24 * a02) - (m04 * a22)) + (m02 * a42));
|
const RT m014 = m01*a42 - m04*a12 + m14*a02;
|
||||||
double m034;
|
const RT m015 = m01*a52 - m05*a12 + m15*a02;
|
||||||
m034 = (((m34 * a02) - (m04 * a32)) + (m03 * a42));
|
const RT m023 = m02*a32 - m03*a22 + m23*a02;
|
||||||
double m123;
|
const RT m024 = m02*a42 - m04*a22 + m24*a02;
|
||||||
m123 = (((m23 * a12) - (m13 * a22)) + (m12 * a32));
|
const RT m025 = m02*a52 - m05*a22 + m25*a02;
|
||||||
double m124;
|
const RT m034 = m03*a42 - m04*a32 + m34*a02;
|
||||||
m124 = (((m24 * a12) - (m14 * a22)) + (m12 * a42));
|
const RT m035 = m03*a52 - m05*a32 + m35*a02;
|
||||||
double m134;
|
const RT m045 = m04*a52 - m05*a42 + m45*a02;
|
||||||
m134 = (((m34 * a12) - (m14 * a32)) + (m13 * a42));
|
const RT m123 = m12*a32 - m13*a22 + m23*a12;
|
||||||
double m234;
|
const RT m124 = m12*a42 - m14*a22 + m24*a12;
|
||||||
m234 = (((m34 * a22) - (m24 * a32)) + (m23 * a42));
|
const RT m125 = m12*a52 - m15*a22 + m25*a12;
|
||||||
double m0123;
|
const RT m134 = m13*a42 - m14*a32 + m34*a12;
|
||||||
m0123 = ((((m123 * a03) - (m023 * a13)) + (m013 * a23)) - (m012 * a33));
|
const RT m135 = m13*a52 - m15*a32 + m35*a12;
|
||||||
double m0124;
|
const RT m145 = m14*a52 - m15*a42 + m45*a12;
|
||||||
m0124 = ((((m124 * a03) - (m024 * a13)) + (m014 * a23)) - (m012 * a43));
|
const RT m234 = m23*a42 - m24*a32 + m34*a22;
|
||||||
double m0134;
|
const RT m235 = m23*a52 - m25*a32 + m35*a22;
|
||||||
m0134 = ((((m134 * a03) - (m034 * a13)) + (m014 * a33)) - (m013 * a43));
|
const RT m245 = m24*a52 - m25*a42 + m45*a22;
|
||||||
double m0234;
|
const RT m345 = m34*a52 - m35*a42 + m45*a32;
|
||||||
m0234 = ((((m234 * a03) - (m034 * a23)) + (m024 * a33)) - (m023 * a43));
|
// Now compute the minors of rank 4
|
||||||
double m1234;
|
const RT m0123 = m012*a33 - m013*a23 + m023*a13 - m123*a03;
|
||||||
m1234 = ((((m234 * a13) - (m134 * a23)) + (m124 * a33)) - (m123 * a43));
|
const RT m0124 = m012*a43 - m014*a23 + m024*a13 - m124*a03;
|
||||||
double m01234;
|
const RT m0125 = m012*a53 - m015*a23 + m025*a13 - m125*a03;
|
||||||
m01234 = (((((m1234 * a04) - (m0234 * a14)) + (m0134 * a24)) - (m0124 * a34)) + (m0123 * a44));
|
const RT m0134 = m013*a43 - m014*a33 + m034*a13 - m134*a03;
|
||||||
return m01234;
|
const RT m0135 = m013*a53 - m015*a33 + m035*a13 - m135*a03;
|
||||||
|
const RT m0145 = m014*a53 - m015*a43 + m045*a13 - m145*a03;
|
||||||
|
const RT m0234 = m023*a43 - m024*a33 + m034*a23 - m234*a03;
|
||||||
|
const RT m0235 = m023*a53 - m025*a33 + m035*a23 - m235*a03;
|
||||||
|
const RT m0245 = m024*a53 - m025*a43 + m045*a23 - m245*a03;
|
||||||
|
const RT m0345 = m034*a53 - m035*a43 + m045*a33 - m345*a03;
|
||||||
|
const RT m1234 = m123*a43 - m124*a33 + m134*a23 - m234*a13;
|
||||||
|
const RT m1235 = m123*a53 - m125*a33 + m135*a23 - m235*a13;
|
||||||
|
const RT m1245 = m124*a53 - m125*a43 + m145*a23 - m245*a13;
|
||||||
|
const RT m1345 = m134*a53 - m135*a43 + m145*a33 - m345*a13;
|
||||||
|
const RT m2345 = m234*a53 - m235*a43 + m245*a33 - m345*a23;
|
||||||
|
// Now compute the minors of rank 5
|
||||||
|
const RT m01234 = m0123*a44 - m0124*a34 + m0134*a24 - m0234*a14 + m1234*a04;
|
||||||
|
const RT m01235 = m0123*a54 - m0125*a34 + m0135*a24 - m0235*a14 + m1235*a04;
|
||||||
|
const RT m01245 = m0124*a54 - m0125*a44 + m0145*a24 - m0245*a14 + m1245*a04;
|
||||||
|
const RT m01345 = m0134*a54 - m0135*a44 + m0145*a34 - m0345*a14 + m1345*a04;
|
||||||
|
const RT m02345 = m0234*a54 - m0235*a44 + m0245*a34 - m0345*a24 + m2345*a04;
|
||||||
|
const RT m12345 = m1234*a54 - m1235*a44 + m1245*a34 - m1345*a24 + m2345*a14;
|
||||||
|
// Now compute the minors of rank 6
|
||||||
|
const RT m012345 = m01234*a55 - m01235*a45 + m01245*a35 - m01345*a25
|
||||||
|
+ m02345*a15 - m12345*a05;
|
||||||
|
return m012345;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
inline int orientationC5( double p0, double p1, double p2, double p3, double p4, double q0, double q1, double q2, double q3, double q4, double r0, double r1, double r2, double r3, double r4, double s0, double s1, double s2, double s3, double s4, double t0, double t1, double t2, double t3, double t4, double u0, double u1, double u2, double u3, double u4) {
|
|
||||||
double pq0;
|
inline int orientationC5( double p0, double p1, double p2, double p3, double p4,
|
||||||
pq0 = (q0 - p0);
|
double q0, double q1, double q2, double q3, double q4,
|
||||||
double pq1;
|
double r0, double r1, double r2, double r3, double r4,
|
||||||
pq1 = (q1 - p1);
|
double s0, double s1, double s2, double s3, double s4,
|
||||||
double pq2;
|
double t0, double t1, double t2, double t3, double t4,
|
||||||
pq2 = (q2 - p2);
|
double u0, double u1, double u2, double u3, double u4)
|
||||||
double pq3;
|
{
|
||||||
pq3 = (q3 - p3);
|
RT det = determinant_with_1_in_row_0(p0, p1, p2, p3, p4,
|
||||||
double pq4;
|
q0, q1, q2, q3, q4,
|
||||||
pq4 = (q4 - p4);
|
r0, r1, r2, r3, r4,
|
||||||
double pr0;
|
s0, s1, s2, s3, s4,
|
||||||
pr0 = (r0 - p0);
|
t0, t1, t2, t3, t4,
|
||||||
double pr1;
|
u0, u1, u2, u3, u4);
|
||||||
pr1 = (r1 - p1);
|
if (det > 0) return 1;
|
||||||
double pr2;
|
if (det < 0) return -1;
|
||||||
pr2 = (r2 - p2);
|
return 0;
|
||||||
double pr3;
|
|
||||||
pr3 = (r3 - p3);
|
|
||||||
double pr4;
|
|
||||||
pr4 = (r4 - p4);
|
|
||||||
double ps0;
|
|
||||||
ps0 = (s0 - p0);
|
|
||||||
double ps1;
|
|
||||||
ps1 = (s1 - p1);
|
|
||||||
double ps2;
|
|
||||||
ps2 = (s2 - p2);
|
|
||||||
double ps3;
|
|
||||||
ps3 = (s3 - p3);
|
|
||||||
double ps4;
|
|
||||||
ps4 = (s4 - p4);
|
|
||||||
double pt0;
|
|
||||||
pt0 = (t0 - p0);
|
|
||||||
double pt1;
|
|
||||||
pt1 = (t1 - p1);
|
|
||||||
double pt2;
|
|
||||||
pt2 = (t2 - p2);
|
|
||||||
double pt3;
|
|
||||||
pt3 = (t3 - p3);
|
|
||||||
double pt4;
|
|
||||||
pt4 = (t4 - p4);
|
|
||||||
double pu0;
|
|
||||||
pu0 = (t0 - p0);
|
|
||||||
double pu1;
|
|
||||||
pu1 = (u1 - p1);
|
|
||||||
double pu2;
|
|
||||||
pu2 = (u2 - p2);
|
|
||||||
double pu3;
|
|
||||||
pu3 = (u3 - p3);
|
|
||||||
double pu4;
|
|
||||||
pu4 = (u4 - p4);
|
|
||||||
double det;
|
|
||||||
det = determinant( pq0, pq1, pq2, pq3, pq4, pr0, pr1, pr2, pr3, pr4, ps0, ps1, ps2, ps3, ps4, pt0, pt1, pt2, pt3, pt4, pu0, pu1, pu2, pu3, pu4 );
|
|
||||||
int int_tmp_result;
|
|
||||||
double eps;
|
|
||||||
double max1 = fabs(pq0);
|
|
||||||
if( (max1 < fabs(pr0)) )
|
|
||||||
{
|
|
||||||
max1 = fabs(pr0);
|
|
||||||
}
|
|
||||||
if( (max1 < fabs(ps0)) )
|
|
||||||
{
|
|
||||||
max1 = fabs(ps0);
|
|
||||||
}
|
|
||||||
if( (max1 < fabs(pt0)) )
|
|
||||||
{
|
|
||||||
max1 = fabs(pt0);
|
|
||||||
}
|
|
||||||
if( (max1 < fabs(pu0)) )
|
|
||||||
{
|
|
||||||
max1 = fabs(pu0);
|
|
||||||
}
|
|
||||||
double max2 = fabs(pq1);
|
|
||||||
if( (max2 < fabs(pr1)) )
|
|
||||||
{
|
|
||||||
max2 = fabs(pr1);
|
|
||||||
}
|
|
||||||
if( (max2 < fabs(ps1)) )
|
|
||||||
{
|
|
||||||
max2 = fabs(ps1);
|
|
||||||
}
|
|
||||||
if( (max2 < fabs(pt1)) )
|
|
||||||
{
|
|
||||||
max2 = fabs(pt1);
|
|
||||||
}
|
|
||||||
if( (max2 < fabs(pu1)) )
|
|
||||||
{
|
|
||||||
max2 = fabs(pu1);
|
|
||||||
}
|
|
||||||
double max3 = fabs(pq2);
|
|
||||||
if( (max3 < fabs(pr2)) )
|
|
||||||
{
|
|
||||||
max3 = fabs(pr2);
|
|
||||||
}
|
|
||||||
if( (max3 < fabs(ps2)) )
|
|
||||||
{
|
|
||||||
max3 = fabs(ps2);
|
|
||||||
}
|
|
||||||
if( (max3 < fabs(pt2)) )
|
|
||||||
{
|
|
||||||
max3 = fabs(pt2);
|
|
||||||
}
|
|
||||||
if( (max3 < fabs(pu2)) )
|
|
||||||
{
|
|
||||||
max3 = fabs(pu2);
|
|
||||||
}
|
|
||||||
double max4 = fabs(pq3);
|
|
||||||
if( (max4 < fabs(pr3)) )
|
|
||||||
{
|
|
||||||
max4 = fabs(pr3);
|
|
||||||
}
|
|
||||||
if( (max4 < fabs(ps3)) )
|
|
||||||
{
|
|
||||||
max4 = fabs(ps3);
|
|
||||||
}
|
|
||||||
if( (max4 < fabs(pt3)) )
|
|
||||||
{
|
|
||||||
max4 = fabs(pt3);
|
|
||||||
}
|
|
||||||
if( (max4 < fabs(pu3)) )
|
|
||||||
{
|
|
||||||
max4 = fabs(pu3);
|
|
||||||
}
|
|
||||||
double max5 = fabs(pq4);
|
|
||||||
if( (max5 < fabs(pr4)) )
|
|
||||||
{
|
|
||||||
max5 = fabs(pr4);
|
|
||||||
}
|
|
||||||
if( (max5 < fabs(ps4)) )
|
|
||||||
{
|
|
||||||
max5 = fabs(ps4);
|
|
||||||
}
|
|
||||||
if( (max5 < fabs(pt4)) )
|
|
||||||
{
|
|
||||||
max5 = fabs(pt4);
|
|
||||||
}
|
|
||||||
if( (max5 < fabs(pu4)) )
|
|
||||||
{
|
|
||||||
max5 = fabs(pu4);
|
|
||||||
}
|
|
||||||
double lower_bound_1;
|
|
||||||
double upper_bound_1;
|
|
||||||
lower_bound_1 = max5;
|
|
||||||
upper_bound_1 = max5;
|
|
||||||
if( (max1 < lower_bound_1) )
|
|
||||||
{
|
|
||||||
lower_bound_1 = max1;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
if( (max1 > upper_bound_1) )
|
|
||||||
{
|
|
||||||
upper_bound_1 = max1;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
if( (max2 < lower_bound_1) )
|
|
||||||
{
|
|
||||||
lower_bound_1 = max2;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
if( (max2 > upper_bound_1) )
|
|
||||||
{
|
|
||||||
upper_bound_1 = max2;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
if( (max3 < lower_bound_1) )
|
|
||||||
{
|
|
||||||
lower_bound_1 = max3;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
if( (max3 > upper_bound_1) )
|
|
||||||
{
|
|
||||||
upper_bound_1 = max3;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
if( (max4 < lower_bound_1) )
|
|
||||||
{
|
|
||||||
lower_bound_1 = max4;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
if( (max4 > upper_bound_1) )
|
|
||||||
{
|
|
||||||
upper_bound_1 = max4;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
if( (lower_bound_1 < 9.99657131447050328602e-60) )
|
|
||||||
{
|
|
||||||
return FPG_UNCERTAIN_VALUE;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
if( (upper_bound_1 > 3.21387608851797912384e+60) )
|
|
||||||
{
|
|
||||||
return FPG_UNCERTAIN_VALUE;
|
|
||||||
}
|
|
||||||
eps = (2.22889232457534740153e-13 * ((((max1 * max2) * max3) * max4) * max5));
|
|
||||||
if( (det > eps) )
|
|
||||||
{
|
|
||||||
int_tmp_result = 1;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
if( (det < -eps) )
|
|
||||||
{
|
|
||||||
int_tmp_result = -1;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
return FPG_UNCERTAIN_VALUE;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
return int_tmp_result;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue