1 // Copyright (c) 2001,2004 INRIA Sophia-Antipolis (France). 2 // All rights reserved. 3 // 4 // This file is part of CGAL (www.cgal.org) 5 // 6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Filtered_kernel/include/CGAL/internal/Static_filters/Coplanar_side_of_bounded_circle_3.h $ 7 // $Id: Coplanar_side_of_bounded_circle_3.h e9d41d7 2020-04-21T10:03:00+02:00 Maxime Gimeno 8 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial 9 // 10 // 11 // Author(s) : Sylvain Pion 12 13 #ifndef CGAL_INTERNAL_STATIC_FILTERS_COPLANAR_SIDE_OF_BOUNDED_CIRCLE_3_H 14 #define CGAL_INTERNAL_STATIC_FILTERS_COPLANAR_SIDE_OF_BOUNDED_CIRCLE_3_H 15 16 #include <CGAL/Simple_cartesian.h> 17 #include <CGAL/Filtered_kernel.h> 18 #include <CGAL/MP_Float.h> 19 #include <CGAL/Profile_counter.h> 20 // #include <CGAL/internal/Static_filters/Static_filter_error.h> // Only used to precompute constants 21 22 namespace CGAL { namespace internal { namespace Static_filters_predicates { 23 24 template <class Point> 25 class Side_of_bounded_circle_3 26 { 27 // Computes the epsilon for In_circle_3. cir_3()28 static double cir_3() 29 { 30 // NOTE : This produces a buggy degree warning. 31 // Maybe there's a bug in the formula. 32 typedef CGAL::Static_filter_error F; 33 F t1 = F(1)-F(1); // First translation 34 F sq = t1*t1+t1*t1+t1*t1; // squares 35 F n1 = t1*t1 - t1*t1; // normal vector 36 F sq_n1 = n1*n1 + n1*n1 + n1*n1; 37 F det = determinant(t1, t1, t1, sq, 38 t1, t1, t1, sq, 39 t1, t1, t1, sq, 40 n1, n1, n1, sq_n1); // Full det 41 double err = det.error(); 42 std::cerr << "*** epsilon for In_circle_3 = " << err << std::endl; 43 return err; 44 } 45 46 public: 47 typedef Bounded_side result_type; 48 operator()49 Bounded_side operator()(const Point &p, const Point &q, 50 const Point &r, const Point &t) const 51 { 52 return opti_coplanar_in_circleC3( 53 to_double(p.x()), to_double(p.y()), to_double(p.z()), 54 to_double(q.x()), to_double(q.y()), to_double(q.z()), 55 to_double(r.x()), to_double(r.y()), to_double(r.z()), 56 to_double(t.x()), to_double(t.y()), to_double(t.z())); 57 } 58 59 Bounded_side opti_coplanar_in_circleC3(double px,double py,double pz,double qx,double qy,double qz,double rx,double ry,double rz,double tx,double ty,double tz)60 opti_coplanar_in_circleC3(double px, double py, double pz, 61 double qx, double qy, double qz, 62 double rx, double ry, double rz, 63 double tx, double ty, double tz) const 64 { 65 CGAL_PROFILER("In_circle_3 calls"); 66 67 double ptx = px - tx; 68 double pty = py - ty; 69 double ptz = pz - tz; 70 double pt2 = CGAL_NTS square(ptx) + CGAL_NTS square(pty) + 71 CGAL_NTS square(ptz); 72 double qtx = qx - tx; 73 double qty = qy - ty; 74 double qtz = qz - tz; 75 double qt2 = CGAL_NTS square(qtx) + CGAL_NTS square(qty) + 76 CGAL_NTS square(qtz); 77 double rtx = rx - tx; 78 double rty = ry - ty; 79 double rtz = rz - tz; 80 double rt2 = CGAL_NTS square(rtx) + CGAL_NTS square(rty) + 81 CGAL_NTS square(rtz); 82 double pqx = qx - px; 83 double pqy = qy - py; 84 double pqz = qz - pz; 85 double prx = rx - px; 86 double pry = ry - py; 87 double prz = rz - pz; 88 double vx = pqy*prz - pqz*pry; 89 double vy = pqz*prx - pqx*prz; 90 double vz = pqx*pry - pqy*prx; 91 double v2 = CGAL_NTS square(vx) + CGAL_NTS square(vy) + 92 CGAL_NTS square(vz); 93 94 double det = determinant(ptx,pty,ptz,pt2, 95 rtx,rty,rtz,rt2, 96 qtx,qty,qtz,qt2, 97 vx,vy,vz,v2); 98 99 // Compute the semi-static bound. 100 double maxx = CGAL::abs(px); 101 double maxy = CGAL::abs(py); 102 double maxz = CGAL::abs(pz); 103 104 double aqx = CGAL::abs(qx); 105 double aqy = CGAL::abs(qy); 106 double aqz = CGAL::abs(qz); 107 108 double arx = CGAL::abs(rx); 109 double ary = CGAL::abs(ry); 110 double arz = CGAL::abs(rz); 111 112 113 double atx = CGAL::abs(tx); 114 double aty = CGAL::abs(ty); 115 double atz = CGAL::abs(tz); 116 117 if (maxx < aqx) maxx = aqx; 118 if (maxx < arx) maxx = arx; 119 if (maxx < atx) maxx = atx; 120 121 if (maxy < aqy) maxy = aqy; 122 if (maxy < ary) maxy = ary; 123 if (maxy < aty) maxy = aty; 124 125 if (maxz < aqz) maxz = aqz; 126 if (maxz < arz) maxz = arz; 127 if (maxz < atz) maxz = atz; 128 129 double d = (std::max)(maxx, (std::max)(maxy, maxz)); 130 double eps = 3.27418e-11 * d * d * d * d * d * d; 131 132 if (det > eps) return ON_BOUNDED_SIDE; 133 if (det < -eps) return ON_UNBOUNDED_SIDE; 134 135 CGAL_PROFILER("In_circle_3 semi-static failures"); 136 137 typedef Filtered_kernel<Simple_cartesian<double>, MP_Float> K; 138 typedef K::Point_3 P; 139 140 return coplanar_side_of_bounded_circle(P(px,py,pz), P(qx,qy,qz), 141 P(rx,ry,rz), P(tx,ty,tz)); 142 } 143 }; 144 145 } } } // namespace CGAL::internal::Static_filters_predicates 146 147 #endif // CGAL_INTERNAL_STATIC_FILTERS_COPLANAR_SIDE_OF_BOUNDED_CIRCLE_3_H 148