Triangulation_3: Store subdeterminant for in_sphere predicate

Cartesian_kernel/include/CGAL/Cartesian/function_objects.h

@@ -27,6 +27,8 @@
 #include <CGAL/Distance_3/Point_3_Point_3.h>
 #include <CGAL/Distance_3/internal/squared_distance_utils_3.h>
 
+#include <array>
+
 namespace CGAL {
 
 namespace CartesianKernelFunctors {
@@ -4380,6 +4382,18 @@ namespace CartesianKernelFunctors {
                                        s.x(), s.y(), s.z(),
                                        test.x(), test.y(), test.z());
     }
+  Oriented_side
+    operator()( const Point_3& p, const Point_3& q, const Point_3& r,
+                const Point_3& s, const Point_3& test,
+                const std::array<double,4>& det) const
+    {
+      return side_of_oriented_sphereC3(p.x(), p.y(), p.z(),
+                                       q.x(), q.y(), q.z(),
+                                       r.x(), r.y(), r.z(),
+                                       s.x(), s.y(), s.z(),
+                                       test.x(), test.y(), test.z(),
+                                       det);
+    }
   };
 
 

Cartesian_kernel/include/CGAL/predicates/kernel_ftC3.h

@@ -21,6 +21,7 @@
 #include <CGAL/predicates/sign_of_determinant.h>
 #include <CGAL/predicates/kernel_ftC2.h>
 #include <CGAL/constructions/kernel_ftC3.h>
+#include <array>
 
 namespace CGAL {
 
@@ -369,6 +370,19 @@ side_of_oriented_sphereC3(const FT &px, const FT &py, const FT &pz,
   // Note that the determinant above is det(P,R,Q,S) (Q and R are swapped)!
 }
 
+inline
+Oriented_side
+side_of_oriented_sphereC3(double px, double py, double pz,
+                          double tx, double ty, double tz,
+                          const std::array<double,4>& det)
+{
+  double ptx = tx - px;
+  double pty = ty - py;
+  double ptz = tz - pz;
+  return sign(ptx * det[0] - pty * det[1] + ptz * det[2] - (ptx*ptx + pty*pty + ptz*ptz) * det[3]);
+}
+
+
 template <class FT >
 CGAL_KERNEL_MEDIUM_INLINE
 typename Same_uncertainty_nt<Bounded_side, FT>::type

Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Side_of_oriented_sphere_3.h

@@ -148,6 +148,122 @@ class Side_of_oriented_sphere_3
       return Base::operator()(p, q, r, s, t);
   }
 
+
+
+  Oriented_side
+  operator()(const Point_3 &p, const Point_3 &q, const Point_3 &r,
+             const Point_3 &s, const Point_3 &t, const std::array<double,4>& subdet) const
+  {
+      CGAL_BRANCH_PROFILER_3("semi-static failures/attempts/calls to   : Side_of_oriented_sphere_3", tmp);
+      double px, py, pz, qx, qy, qz, rx, ry, rz, sx, sy, sz, tx, ty, tz;
+
+      if (fit_in_double(p.x(), px) && fit_in_double(p.y(), py) &&
+          fit_in_double(p.z(), pz) &&
+          fit_in_double(q.x(), qx) && fit_in_double(q.y(), qy) &&
+          fit_in_double(q.z(), qz) &&
+          fit_in_double(r.x(), rx) && fit_in_double(r.y(), ry) &&
+          fit_in_double(r.z(), rz) &&
+          fit_in_double(s.x(), sx) && fit_in_double(s.y(), sy) &&
+          fit_in_double(s.z(), sz) &&
+          fit_in_double(t.x(), tx) && fit_in_double(t.y(), ty) &&
+          fit_in_double(t.z(), tz))
+      {
+          CGAL_BRANCH_PROFILER_BRANCH_1(tmp);
+
+          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 rtx = rx - tx;
+          double rty = ry - ty;
+          double rtz = rz - tz;
+
+          double stx = sx - tx;
+          double sty = sy - ty;
+          double stz = sz - tz;
+
+          // Compute the semi-static bound.
+          double maxx = CGAL::abs(ptx);
+          double maxy = CGAL::abs(pty);
+          double maxz = CGAL::abs(ptz);
+
+          double aqtx = CGAL::abs(qtx);
+          double artx = CGAL::abs(rtx);
+          double astx = CGAL::abs(stx);
+
+          double aqty = CGAL::abs(qty);
+          double arty = CGAL::abs(rty);
+          double asty = CGAL::abs(sty);
+
+          double aqtz = CGAL::abs(qtz);
+          double artz = CGAL::abs(rtz);
+          double astz = CGAL::abs(stz);
+
+#ifdef CGAL_USE_SSE2_MAX
+          CGAL::Max<double> mmax;
+          maxx = mmax(maxx, aqtx, artx, astx);
+          maxy = mmax(maxy, aqty, arty, asty);
+          maxz = mmax(maxz, aqtz, artz, astz);
+#else
+          if (maxx < aqtx) maxx = aqtx;
+          if (maxx < artx) maxx = artx;
+          if (maxx < astx) maxx = astx;
+
+          if (maxy < aqty) maxy = aqty;
+          if (maxy < arty) maxy = arty;
+          if (maxy < asty) maxy = asty;
+
+          if (maxz < aqtz) maxz = aqtz;
+          if (maxz < artz) maxz = artz;
+          if (maxz < astz) maxz = astz;
+#endif
+
+          double eps = 1.2466136531027298e-13 * maxx * maxy * maxz;
+
+#ifdef CGAL_USE_SSE2_MAX
+          /*
+          CGAL::Min<double> mmin;
+          double tmp = mmin(maxx, maxy, maxz);
+          maxz = mmax(maxx, maxy, maxz);
+          maxx = tmp;
+          */
+          sse2minmax(maxx,maxy,maxz);
+          // maxy can contain ANY element
+
+#else
+          // Sort maxx < maxy < maxz.
+          if (maxx > maxz)
+              std::swap(maxx, maxz);
+          if (maxy > maxz)
+              std::swap(maxy, maxz);
+          else if (maxy < maxx)
+              std::swap(maxx, maxy);
+#endif
+          double det =  - ptx * subdet[0] + pty * subdet[1] - ptz * subdet[2] -pt2 * subdet[3];
+
+          // Protect against underflow in the computation of eps.
+          if (maxx < 1e-58) /* sqrt^5(min_double/eps) */ {
+            if (maxx == 0)
+              return ON_ORIENTED_BOUNDARY;
+          }
+          // Protect against overflow in the computation of det.
+          else if (maxz < 1e61) /* sqrt^5(max_double/4 [hadamard]) */ {
+            eps *= (maxz * maxz);
+            if (det > eps)  return ON_POSITIVE_SIDE;
+            if (det < -eps) return ON_NEGATIVE_SIDE;
+          }
+
+          CGAL_BRANCH_PROFILER_BRANCH_2(tmp);
+      }
+      return Base::operator()(p, q, r, s, t);
+  }
+
+
   // Computes the epsilon for Side_of_oriented_sphere_3.
   static double compute_epsilon()
   {

Kernel_23/include/CGAL/determinant.h

@@ -20,6 +20,8 @@
 
 #include <CGAL/kernel_config.h>
 
+#include <array>
+
 namespace CGAL {
 
 template <class RT>
@@ -399,6 +401,39 @@ determinant(
   return m0123456;
 }
 
+inline void subdeterminants(
+ double ax,  double ay,  double az,
+ double bx,  double by,  double bz,
+ double cx,  double cy,  double cz,
+ double dx,  double dy,  double dz,
+ std::array<double,4>& det)
+{
+  double bax = bx-ax;
+  double bay = by-ay;
+  double baz = bz-az;
+  double cax = cx-ax;
+  double cay = cy-ay;
+  double caz = cz-az;
+  double dax = dx-ax;
+  double day = dy-ay;
+  double daz = dz-az;
+
+  double ba2 = bax*bax + bay*bay + baz*baz;
+  double ca2 = cax*cax + cay*cay + caz*caz;
+  double da2 = dax*dax + day*day + daz*daz;
+
+  double bc = baz*ca2 - caz*ba2;
+  double bd = baz*da2 - daz*ba2;
+  double cd = caz*da2 - daz*ca2;
+  det[0] = bay*cd - cay*bd + day*bc;
+  det[1] = bax*cd - cax*bd + dax*bc;
+
+  double sd0 = bax*cay - bay*cax;
+  double sd1 = bax*day - bay*dax;
+  double sd2 = cax*day - cay*dax;
+  det[2] = ba2 * sd2 - ca2 * sd1 + da2 * sd0;
+  det[3] = baz * sd2 - caz * sd1 + daz * sd0;
+}
 
 } //namespace CGAL
 

TDS_3/include/CGAL/Triangulation_data_structure_3.h

@@ -664,11 +664,11 @@ class Triangulation_data_structure_3
       f.first->set_neighbor(f.second, nc);
 
       vertex_pair_facet_map.set({u, v}, {local_facet_index,
-                                         static_cast<unsigned char>(nc->index(w))});
+                                         static_cast<unsigned char>(2)});
       vertex_pair_facet_map.set({v, w}, {local_facet_index,
-                                         static_cast<unsigned char>(nc->index(u))});
+                                         static_cast<unsigned char>(1)});
       vertex_pair_facet_map.set({w, u}, {local_facet_index,
-                                         static_cast<unsigned char>(nc->index(v))});
+                                         static_cast<unsigned char>(0)});
     }
 
     for(auto it = vertex_pair_facet_map.begin(); it != vertex_pair_facet_map.end(); ++it){

Triangulation_3/benchmark/Triangulation_3/CMakeLists.txt

@@ -15,6 +15,7 @@ create_single_source_cgal_program("ocean.cpp")
 create_single_source_cgal_program("incident_edges.cpp")
 create_single_source_cgal_program("simple_2.cpp")
 create_single_source_cgal_program("simple.cpp")
+create_single_source_cgal_program("DT_storing_subdeterminant.cpp")
 create_single_source_cgal_program("Triangulation_benchmark_3.cpp")
 create_single_source_cgal_program("segment_traverser_benchmark.cpp" )
 

Triangulation_3/benchmark/Triangulation_3/DT_storing_subdeterminant.cpp

@@ -0,0 +1,50 @@
+#define CGAL_DELAUNAY_3_USE_SUBDETERMINANTS 1
+// #define CGAL_PROFILE
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Delaunay_triangulation_3.h>
+#include <CGAL/Triangulation_cell_base_with_info_3.h>
+#include <CGAL/Timer.h>
+#include <vector>
+#include <array>
+#include <iostream>
+#include <string>
+#include <fstream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel         K;
+typedef CGAL::Triangulation_vertex_base_3<K>                        Vb;
+#ifdef CGAL_DELAUNAY_3_USE_SUBDETERMINANTS
+typedef CGAL::Triangulation_cell_base_with_info_3<std::array<double,4>, K>    Cbb;
+typedef CGAL::Delaunay_triangulation_cell_base_3<K,Cbb>             Cb;
+#else
+typedef CGAL::Delaunay_triangulation_cell_base_3<K>                 Cb;
+#endif
+typedef CGAL::Triangulation_data_structure_3<Vb, Cb>                Tds;
+typedef CGAL::Delaunay_triangulation_3<K, Tds>                      Delaunay;
+typedef Delaunay::Point                                             Point;
+typedef CGAL::Timer                                                 Timer;
+
+int main(int argc, char* argv[])
+{
+  const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("points_3/ocean_r.xyz");
+  std::ifstream in(filename.c_str());
+  std::vector<Point> points;
+  Point p;
+
+  while(in >> p){
+    points.push_back(p);
+  }
+
+  std::cout << points.size() << " points read\n";
+
+  Timer timer;
+  timer.start();
+  {
+    Delaunay dt;
+    dt.insert(points.begin(), points.end());
+  }
+  std::cerr << timer.time() << " sec" << std::endl;
+  std::cout << "done" << std::endl;
+  return 0;
+}
+