- Added Coplanar_side_of_bounded_circle_3 to Static_filters.

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.
 

Packages/Interval_arithmetic/include/CGAL/Static_filters.h

@@ -38,6 +38,7 @@
 #include <CGAL/Static_filters/Side_of_oriented_circle_2.h>
 #include <CGAL/Static_filters/Side_of_oriented_sphere_3.h>
 #include <CGAL/Static_filters/Coplanar_orientation_3.h>
+#include <CGAL/Static_filters/Coplanar_side_of_bounded_circle_3.h>
 
 // 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<Point_3>     Side_of_oriented_sphere_3;
   typedef SF_Coplanar_orientation_3<Point_3, Orientation_2>
                                                     Coplanar_orientation_3;
+  typedef SF_Side_of_bounded_circle_3<Point_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);
       }
   }
 
@@ -149,6 +157,7 @@ private:
   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

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 (<Mariette.Yvinec@sophia.inria.fr>)
+//
+// ============================================================================
+
+#ifndef CGAL_STATIC_FILTERS_COPLANAR_SIDE_OF_BOUNDED_CIRCLE_3_H
+#define CGAL_STATIC_FILTERS_COPLANAR_SIDE_OF_BOUNDED_CIRCLE_3_H
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Filtered_exact.h>
+#include <CGAL/MP_Float.h>
+#include <CGAL/Profile_counter.h>
+// #include <CGAL/Static_filter_error.h> // Only used to precompute constants
+
+CGAL_BEGIN_NAMESPACE
+
+template <class Point>
+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<Filtered_exact<double, MP_Float> > 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 <class Point>
+const double SF_Side_of_bounded_circle_3<Point>::epsilon = 3.27418e-11;
+
+CGAL_END_NAMESPACE
+
+#endif // CGAL_STATIC_FILTERS_COPLANAR_SIDE_OF_BOUNDED_CIRCLE_3_H

Packages/Interval_arithmetic/test/Interval_arithmetic/Static_filters.C

@@ -1,29 +1,24 @@
 
 #define CGAL_PROFILE
 
-#include <CGAL/Random.h>
+#include <CGAL/Simple_cartesian.h>
 #include <CGAL/Static_filters.h>
+#include <CGAL/point_generators_2.h>
+#include <CGAL/point_generators_3.h>
 
 typedef CGAL::Simple_cartesian<double> K;
+typedef CGAL::Random_points_in_square_2<K::Point_2> Rand_2;
 
-CGAL::Random R;
+CGAL::Orientation   oo;
-
-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::Oriented_side os;
+CGAL::Bounded_side  bs;
 
 int main()
 {
+  CGAL::Random_points_in_square_2<K::Point_2> Rand_2(1000.0);
+  CGAL::Random_points_in_cube_3<K::Point_3>   Rand_3(1000.0);
   CGAL::Static_filters<K> k;
+
   const CGAL::Static_filters<K>::Orientation_3 & my_o3 =
         k.orientation_3_object();
   const CGAL::Static_filters<K>::Orientation_2 & my_o2 =
@@ -34,32 +29,31 @@ int main()
         k.side_of_oriented_sphere_3_object();
   const CGAL::Static_filters<K>::Coplanar_orientation_3 & my_co3 =
         k.coplanar_orientation_3_object();
+  const CGAL::Static_filters<K>::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_3 p = *Rand_3++; k.register_object(p);
-    K::Point_2 p2(prand2()), q2(prand2()), r2(prand2()), s2(prand2());
+    K::Point_3 q = *Rand_3++; k.register_object(q);
-    k.register_object(p);
+    K::Point_3 r = *Rand_3++; k.register_object(r);
-    k.register_object(q);
+    K::Point_3 s = *Rand_3++; k.register_object(s);
-    k.register_object(r);
+    K::Point_3 t = *Rand_3++; k.register_object(t);
-    k.register_object(s);
+
-    k.register_object(t);
+    K::Point_2 p2 = *Rand_2++; k.register_object(p2);
-    k.register_object(p2);
+    K::Point_2 q2 = *Rand_2++; k.register_object(q2);
-    k.register_object(q2);
+    K::Point_2 r2 = *Rand_2++; k.register_object(r2);
-    k.register_object(r2);
+    K::Point_2 s2 = *Rand_2++; k.register_object(s2);
-    k.register_object(s2);
+
-    // CGAL::Protect_FPU_rounding<false> Z;
     for (int j=0; j<1000; ++j) {
-      // assert( my_o(p, q, r, s) == ore(pe, qe, re, se) );
+      oo = my_o3(p, q, r, s);
-      ooo = my_o3(p, q, r, s);      //  2.62 s -> 2.22 s
+      oo = my_o2(p2, q2, r2);
-      ooo = my_o2(p2, q2, r2);      //  2.62 s -> 2.22 s
       os = my_c2(p2, q2, r2, s2);
       os = my_s3(p, q, r, s, t);
-      ooo = my_co3(p, q, r);
+      oo = my_co3(p, q, r);
-      ooo = my_co3(p, q, r, s);
+      oo = my_co3(p, q, r, s);
-
+      bs = my_cobc3(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
     }
   }