Reformulate a cross product to increase precision

Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_surface_mesh.cpp

@@ -16,10 +16,12 @@ namespace SMS = CGAL::Surface_mesh_simplification;
 
 int main(int argc, char** argv)
 {
+  std::cout.precision(17);
+  std::cerr.precision(17);
+
   Surface_mesh surface_mesh;
   const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/cube-meshed.off");
-  std::ifstream is(filename);
+  if(!CGAL::IO::read_polygon_mesh(filename, surface_mesh))
-  if(!is || !(is >> surface_mesh))
   {
     std::cerr << "Failed to read input mesh: " << filename << std::endl;
     return EXIT_FAILURE;
@@ -33,6 +35,8 @@ int main(int argc, char** argv)
 
   std::chrono::steady_clock::time_point start_time = std::chrono::steady_clock::now();
 
+  std::cout << num_vertices(surface_mesh) << " vertices, " << num_edges(surface_mesh) << " edges (BEFORE)" << std::endl;
+
   // In this example, the simplification stops when the number of undirected edges
   // drops below 10% of the initial count
   double stop_ratio = (argc > 2) ? std::stod(argv[2]) : 0.1;
@@ -45,7 +49,7 @@ int main(int argc, char** argv)
   std::cout << "\nFinished!\n" << r << " edges removed.\n" << surface_mesh.number_of_edges() << " final edges.\n";
   std::cout << "Time elapsed: " << std::chrono::duration_cast<std::chrono::milliseconds>(end_time - start_time).count() << "ms" << std::endl;
 
-  CGAL::IO::write_polygon_mesh((argc > 3) ? argv[3] : "out.off", surface_mesh, CGAL::parameters::stream_precision(17));
+  CGAL::IO::write_polygon_mesh((argc > 3) ? argv[3] : "out_SC.off", surface_mesh, CGAL::parameters::stream_precision(17));
 
   return EXIT_SUCCESS;
 }

Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h

@@ -131,6 +131,13 @@ public:
     return Vector_3(v*m.r0(), v*m.r1(), v*m.r2());
   }
 
+  friend std::ostream& operator<<(std::ostream & os, const MatrixC33& m)
+{
+  return os << m.r0() << std::endl
+            << m.r1() << std::endl
+            << m.r2() << std::endl;
+}
+
   RT determinant() const
   {
     return CGAL::determinant(r0().x(), r0().y(), r0().z(),

Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h

@@ -13,6 +13,8 @@
 
 #include <CGAL/license/Surface_mesh_simplification.h>
 
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
 #include <CGAL/Surface_mesh_simplification/internal/Common.h>
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/LindstromTurk_params.h>
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_profile.h>
@@ -119,9 +121,103 @@ private :
   const Geom_traits& geom_traits() const { return mProfile.geom_traits(); }
   const TM& surface() const { return mProfile.surface(); }
 
+#if 0
+  // a*b - c*d
+  // The next two functions are from https://stackoverflow.com/questions/63665010/accurate-floating-point-computation-of-the-sum-and-difference-of-two-products
+  static double diff_of_products_kahan(const double a, const double b, const double c, const double d)
+  {
+    double w = d * c;
+    double e = std::fma(c, -d, w);
+    double f = std::fma(a, b, -w);
+    return f + e;
+  }
+
+  static double diff_of_products_cht(const double a, const double b, const double c, const double d)
+  {
+    double p1 = a * b;
+    double p2 = c * d;
+    double e1 = std::fma (a, b, -p1);
+    double e2 = std::fma (c, -d, p2);
+    double r = p1 - p2;
+    double e = e1 + e2;
+    return r + e;
+  }
+
+  static double diff_of_products(const double a, const double b, const double c, const double d)
+  {
+    // the next two are equivalent in results and speed
+    return diff_of_products_kahan(a, b, c, d);
+    // return diff_of_products_cht(a, b, c, d);
+  }
+
+  template <typename OFT>
+  static OFT diff_of_products(const OFT& a, const OFT& b, const OFT& c, const OFT& d)
+  {
+    return a*b - c*d;
+  }
+#endif
+
+  static Vector SL_cross_product(const Vector& a, const Vector& b)
+  {
+    const FT ax=a.x(), ay=a.y(), az=a.z();
+    const FT bx=b.x(), by=b.y(), bz=b.z();
+
+    auto minor = [](double ai, double bi, double aj, double bj)
+    {
+      // The main idea is that we expect ai and bi (and aj and bj) to have roughly the same magnitude
+      // since this function is used to compute the cross product of two vectors that are defined
+      // as (ORIGIN, pa) and (ORIGIN, pb) and pa and pb are part of the same triangle.
+      //
+      // We can abuse this fact to trade 2 extra subtractions to lower the error.
+      return ai * (bj - aj) + aj * (ai - bi);
+    };
+
+    // ay*
+    FT x = minor(ay, by, az, bz);
+    FT y = minor(az, bz, ax, bx);
+    FT z = minor(ax, bx, ay, by);
+
+    return Vector(x, y, z);
+  }
+
+#if 0
+  static Vector exact_cross_product(const Vector& a, const Vector& b)
+  {
+    CGAL::Cartesian_converter<Geom_traits, CGAL::Exact_predicates_exact_constructions_kernel> to_exact;
+    CGAL::Cartesian_converter<CGAL::Exact_predicates_exact_constructions_kernel, Geom_traits> to_approx;
+    auto exv = cross_product(to_exact(a), to_exact(b));
+    exv.exact();
+    return to_approx(exv);
+  }
+#endif
+
+  static Vector X_product(const Vector& u, const Vector& v)
+  {
+#if 0
+    // this can create large errors and spiky meshes for kernels with inexact constructions
+    return CGAL::cross_product(u,v);
+#elif 0
+    // improves the problem mentioned above a bit, but not enough
+    return { std::fma(u.y(), v.z(), -u.z()*v.y()),
+             std::fma(u.z(), v.x(), -u.x()*v.z()),
+             std::fma(u.x(), v.y(), -u.y()*v.x()) };
+#elif 0
+    // this is the best without resorting to exact, but it inflicts a 20% slowdown
+    return { diff_of_products(u.y(), v.z(), u.z(), v.y()),
+             diff_of_products(u.z(), v.x(), u.x(), v.z()),
+             diff_of_products(u.x(), v.y(), u.y(), v.x()) };
+#elif 1
+    // balanced solution based on abusing the fact that here we expect u and v to have similar coordinates
+    return SL_cross_product(u, v);
+#elif 0
+    // obviously too slow
+    return exact_cross_product(u, v);
+#endif
+  }
+
   static Vector point_cross_product(const Point& a, const Point& b)
   {
-    return cross_product(a-ORIGIN, b-ORIGIN);
+    return X_product(a-ORIGIN, b-ORIGIN);
   }
 
   // This is the (uX)(Xu) product described in the Lindstrom-Turk paper

Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/internal/Common.h

@@ -114,10 +114,7 @@ inline std::string optional_to_string(const std::optional<T>& o) {
 
   namespace internal { namespace { bool cgal_enable_sms_trace = true; } }
   #define CGAL_SMS_TRACE_IMPL(m) \
-    if(::internal::cgal_enable_sms_trace) { \
+    std::cerr << m << std::endl;
-    std::ostringstream ss; ss << m; std::string s = ss.str(); \
-    /*Surface_simplification_external_trace(s)*/ std::cerr << s << std::endl; \
-    }
 
   #define CGAL_SMS_DEBUG_CODE(code) code
 #else