diff --git a/Packages/Interval_arithmetic/changes.txt b/Packages/Interval_arithmetic/changes.txt index e0485971829..8976acf7311 100644 --- a/Packages/Interval_arithmetic/changes.txt +++ b/Packages/Interval_arithmetic/changes.txt @@ -1,3 +1,6 @@ +Version 4.122 on 5 October 2001 +- Added Coplanar_side_of_bounded_circle_3 to Static_filters. + Version 4.121 on 4 October 2001 - Added Coplanar_orientation_3 to Static_filters. diff --git a/Packages/Interval_arithmetic/include/CGAL/Static_filters.h b/Packages/Interval_arithmetic/include/CGAL/Static_filters.h index a7e21294157..3c6b3fb0928 100644 --- a/Packages/Interval_arithmetic/include/CGAL/Static_filters.h +++ b/Packages/Interval_arithmetic/include/CGAL/Static_filters.h @@ -38,6 +38,7 @@ #include #include #include +#include // This traits class gathers optimized predicates written by hand, using // a few steps of filtering. It should work if the initial traits has @@ -78,6 +79,8 @@ public : typedef SF_Side_of_oriented_sphere_3 Side_of_oriented_sphere_3; typedef SF_Coplanar_orientation_3 Coplanar_orientation_3; + typedef SF_Side_of_bounded_circle_3 + Coplanar_side_of_bounded_circle_3; const Orientation_2 & orientation_2_object() const @@ -99,6 +102,10 @@ public : coplanar_orientation_3_object() const { return _coplanar_orientation_3; } + const Coplanar_side_of_bounded_circle_3 & + coplanar_side_of_bounded_circle_3_object() const + { return _coplanar_side_of_bounded_circle_3; } + // These should not be const, but unfortunately Triangulation_?.geom_traits() // only give a const& access (should this be changed ?). // In the mean time, I put the data members mutable. @@ -134,6 +141,7 @@ public : _orientation_3.update(max3x, max3y, max3z); _side_of_oriented_sphere_3.update(max3x, max3y, max3z); _coplanar_orientation_3.update(max3x, max3y, max3z); + _coplanar_side_of_bounded_circle_3.update(max3x, max3y, max3z); } } @@ -144,11 +152,12 @@ private: // A data member for each predicate. // Their state is related to the state of *this. - mutable Orientation_2 _orientation_2; - mutable Orientation_3 _orientation_3; - mutable Side_of_oriented_circle_2 _side_of_oriented_circle_2; - mutable Side_of_oriented_sphere_3 _side_of_oriented_sphere_3; - mutable Coplanar_orientation_3 _coplanar_orientation_3; + mutable Orientation_2 _orientation_2; + mutable Orientation_3 _orientation_3; + mutable Side_of_oriented_circle_2 _side_of_oriented_circle_2; + mutable Side_of_oriented_sphere_3 _side_of_oriented_sphere_3; + mutable Coplanar_orientation_3 _coplanar_orientation_3; + mutable Coplanar_side_of_bounded_circle_3 _coplanar_side_of_bounded_circle_3; }; CGAL_END_NAMESPACE diff --git a/Packages/Interval_arithmetic/include/CGAL/Static_filters/Coplanar_side_of_bounded_circle_3.h b/Packages/Interval_arithmetic/include/CGAL/Static_filters/Coplanar_side_of_bounded_circle_3.h new file mode 100644 index 00000000000..18d6e7d0c4f --- /dev/null +++ b/Packages/Interval_arithmetic/include/CGAL/Static_filters/Coplanar_side_of_bounded_circle_3.h @@ -0,0 +1,177 @@ +// ============================================================================ +// +// Copyright (c) 2001 The CGAL Consortium +// +// This software and related documentation is part of an INTERNAL release +// of the Computational Geometry Algorithms Library (CGAL). It is not +// intended for general use. +// +// ---------------------------------------------------------------------------- +// +// release : +// release_date : +// +// file : include/CGAL/Static_filters/Coplanar_side_of_bounded_circle_3.h +// revision : $Revision$ +// revision_date : $Date$ +// package : Interval Arithmetic +// author(s) : Sylvain Pion +// coordinator : INRIA Sophia-Antipolis () +// +// ============================================================================ + +#ifndef CGAL_STATIC_FILTERS_COPLANAR_SIDE_OF_BOUNDED_CIRCLE_3_H +#define CGAL_STATIC_FILTERS_COPLANAR_SIDE_OF_BOUNDED_CIRCLE_3_H + +#include +#include +#include +#include +// #include // Only used to precompute constants + +CGAL_BEGIN_NAMESPACE + +template +class SF_Side_of_bounded_circle_3 +{ + double _static_epsilon; + + // Computes the epsilon for In_circle_3. + static double cir_3() + { + // NOTE : This produces a buggy degree warning. + // Maybe there's a bug in the formula. + typedef CGAL::Static_filter_error F; + F t1 = F(1)-F(1); // First translation + F sq = t1*t1+t1*t1+t1*t1; // squares + F n1 = t1*t1 - t1*t1; // normal vector + F sq_n1 = n1*n1 + n1*n1 + n1*n1; + F det = det4x4_by_formula(t1, t1, t1, sq, + t1, t1, t1, sq, + t1, t1, t1, sq, + n1, n1, n1, sq_n1); // Full det + double err = det.error(); + std::cerr << "*** epsilon for In_circle_3 = " << err << std::endl; + return err; + } + + static const double epsilon; + +protected: + + template < class R > + friend class Static_filters; + + // These operations are reserved to Static_filters<>, because the context of + // a predicate is linked to the one of the Static_filter<> it is a member of. + SF_Side_of_bounded_circle_3(const SF_Side_of_bounded_circle_3 &) {} + + SF_Side_of_bounded_circle_3& operator=(const SF_Side_of_bounded_circle_3 &) + {} + + SF_Side_of_bounded_circle_3() + { + _static_epsilon = HUGE_VAL; + } + +public: + typedef Bounded_side result_type; + + void update(double dx, double dy, double dz) + { + double d = std::max(std::max(dx, dy), dz); + _static_epsilon = epsilon*d*d*d*d*d*d; + } + + Bounded_side operator()(const Point &p, const Point &q, + const Point &r, const Point &t) const + { + return opti_coplanar_in_circleC3( + to_double(p.x()), to_double(p.y()), to_double(p.z()), + to_double(q.x()), to_double(q.y()), to_double(q.z()), + to_double(r.x()), to_double(r.y()), to_double(r.z()), + to_double(t.x()), to_double(t.y()), to_double(t.z())); + } + + 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) const + { + CGAL_PROFILER(calls, "In_circle_3 calls") + + double ptx = px - tx; + double pty = py - ty; + double ptz = pz - tz; + double pt2 = CGAL_NTS square(ptx) + CGAL_NTS square(pty) + + CGAL_NTS square(ptz); + double qtx = qx - tx; + double qty = qy - ty; + double qtz = qz - tz; + double qt2 = CGAL_NTS square(qtx) + CGAL_NTS square(qty) + + CGAL_NTS square(qtz); + double rtx = rx - tx; + double rty = ry - ty; + double rtz = rz - tz; + double rt2 = CGAL_NTS square(rtx) + CGAL_NTS square(rty) + + CGAL_NTS square(rtz); + double pqx = qx - px; + double pqy = qy - py; + double pqz = qz - pz; + double prx = rx - px; + double pry = ry - py; + double prz = rz - pz; + double vx = pqy*prz - pqz*pry; + double vy = pqz*prx - pqx*prz; + double vz = pqx*pry - pqy*prx; + double v2 = CGAL_NTS square(vx) + CGAL_NTS square(vy) + + CGAL_NTS square(vz); + + double det = det4x4_by_formula(ptx,pty,ptz,pt2, + rtx,rty,rtz,rt2, + qtx,qty,qtz,qt2, + vx,vy,vz,v2); + + // Try a fully static bound first. + if (det > _static_epsilon) return ON_BOUNDED_SIDE; + if (det < -_static_epsilon) return ON_UNBOUNDED_SIDE; + + CGAL_PROFILER(st_fail, "In_circle_3 static failures") + + // Compute the semi-static bound. + double maxx = fabs(px); + if (maxx < fabs(qx)) maxx = fabs(qx); + if (maxx < fabs(rx)) maxx = fabs(rx); + if (maxx < fabs(tx)) maxx = fabs(tx); + double maxy = fabs(py); + if (maxy < fabs(qy)) maxy = fabs(qy); + if (maxy < fabs(ry)) maxy = fabs(ry); + if (maxy < fabs(ty)) maxy = fabs(ty); + double maxz = fabs(pz); + if (maxz < fabs(qz)) maxz = fabs(qz); + if (maxz < fabs(rz)) maxz = fabs(rz); + if (maxz < fabs(tz)) maxz = fabs(tz); + + double d = std::max(maxx, std::max(maxy, maxz)); + double eps = epsilon*d*d*d*d*d*d; + + if (det > eps) return ON_BOUNDED_SIDE; + if (det < -eps) return ON_UNBOUNDED_SIDE; + + CGAL_PROFILER(fail, "In_circle_3 semi-static failures") + + typedef Simple_cartesian > K; + typedef K::Point_3 P; + + return coplanar_side_of_bounded_circle(P(px,py,pz), P(qx,qy,qz), + P(rx,ry,rz), P(tx,ty,tz)); + } +}; + +template +const double SF_Side_of_bounded_circle_3::epsilon = 3.27418e-11; + +CGAL_END_NAMESPACE + +#endif // CGAL_STATIC_FILTERS_COPLANAR_SIDE_OF_BOUNDED_CIRCLE_3_H diff --git a/Packages/Interval_arithmetic/test/Interval_arithmetic/Static_filters.C b/Packages/Interval_arithmetic/test/Interval_arithmetic/Static_filters.C index 12e3be8ccde..d1b15d1191e 100644 --- a/Packages/Interval_arithmetic/test/Interval_arithmetic/Static_filters.C +++ b/Packages/Interval_arithmetic/test/Interval_arithmetic/Static_filters.C @@ -1,29 +1,24 @@ #define CGAL_PROFILE -#include +#include #include +#include +#include typedef CGAL::Simple_cartesian K; +typedef CGAL::Random_points_in_square_2 Rand_2; -CGAL::Random R; - -K::Point_3 prand() -{ - return K::Point_3(R.get_double(), R.get_double(), R.get_double()); -} - -K::Point_2 prand2() -{ - return K::Point_2(R.get_double(), R.get_double()); -} - -CGAL::Orientation ooo; +CGAL::Orientation oo; CGAL::Oriented_side os; +CGAL::Bounded_side bs; int main() { + CGAL::Random_points_in_square_2 Rand_2(1000.0); + CGAL::Random_points_in_cube_3 Rand_3(1000.0); CGAL::Static_filters k; + const CGAL::Static_filters::Orientation_3 & my_o3 = k.orientation_3_object(); const CGAL::Static_filters::Orientation_2 & my_o2 = @@ -34,32 +29,31 @@ int main() k.side_of_oriented_sphere_3_object(); const CGAL::Static_filters::Coplanar_orientation_3 & my_co3 = k.coplanar_orientation_3_object(); + const CGAL::Static_filters::Coplanar_side_of_bounded_circle_3 & my_cobc3 = + k.coplanar_side_of_bounded_circle_3_object(); + + // my_cobc3.cir_3(); for (int i=0; i<100; ++i) { - K::Point_3 p(prand()), q(prand()), r(prand()), s(prand()), t(prand()); - K::Point_2 p2(prand2()), q2(prand2()), r2(prand2()), s2(prand2()); - k.register_object(p); - k.register_object(q); - k.register_object(r); - k.register_object(s); - k.register_object(t); - k.register_object(p2); - k.register_object(q2); - k.register_object(r2); - k.register_object(s2); - // CGAL::Protect_FPU_rounding Z; + K::Point_3 p = *Rand_3++; k.register_object(p); + K::Point_3 q = *Rand_3++; k.register_object(q); + K::Point_3 r = *Rand_3++; k.register_object(r); + K::Point_3 s = *Rand_3++; k.register_object(s); + K::Point_3 t = *Rand_3++; k.register_object(t); + + K::Point_2 p2 = *Rand_2++; k.register_object(p2); + K::Point_2 q2 = *Rand_2++; k.register_object(q2); + K::Point_2 r2 = *Rand_2++; k.register_object(r2); + K::Point_2 s2 = *Rand_2++; k.register_object(s2); + for (int j=0; j<1000; ++j) { - // assert( my_o(p, q, r, s) == ore(pe, qe, re, se) ); - ooo = my_o3(p, q, r, s); // 2.62 s -> 2.22 s - ooo = my_o2(p2, q2, r2); // 2.62 s -> 2.22 s + oo = my_o3(p, q, r, s); + oo = my_o2(p2, q2, r2); os = my_c2(p2, q2, r2, s2); os = my_s3(p, q, r, s, t); - ooo = my_co3(p, q, r); - ooo = my_co3(p, q, r, s); - - // ooo = ore(pe, qe, re, se); // 15.07 s (!prot:13.5) static: 4.35 s (!p:3 -> 1.45) - // ooo = orf(pf, qf, rf, sf); // 0.68 s -> 0.85 s - // ooo = inexact_o(p, q, r, s); // 1.66 s -> 0.20 s + oo = my_co3(p, q, r); + oo = my_co3(p, q, r, s); + bs = my_cobc3(p, q, r, s); } }