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(!is || !(is >> surface_mesh))
+  if(!CGAL::IO::read_polygon_mesh(filename, 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;

Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h

@@ -131,6 +131,13 @@ class MatrixC33
     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

@@ -119,9 +119,133 @@ 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
+
+#ifdef __AVX__
+  static Vector SL_cross_product_avx(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();
+
+    __m256d a = _mm256_set_pd(ay, az, ax, 1.0);
+    __m256d b = _mm256_set_pd(bz, bx, by, 1.0);
+    __m256d c = _mm256_set_pd(az, ax, ay, 1.0);
+    __m256d d = _mm256_set_pd(by, bz, bx, 1.0);
+
+    __m256d s1 = _mm256_sub_pd(b, c);
+    __m256d s2 = _mm256_sub_pd(a, d);
+
+    b = _mm256_mul_pd(a, s1);
+    d = _mm256_mul_pd(c, s2);
+    a = _mm256_add_pd(b, d);
+
+    double res[4];
+    _mm256_storeu_pd(res, a);
+
+//            a  * (b  - c ) + c  * ( a - d);
+//    FT x =  ay * (bz - az) + az * (ay - by);
+//    FT y =  az * (bx - ax) + ax * (az - bz);
+//    FT z =  ax * (by - ay) + ay * (ax - bx);
+
+    return Vector(res[3], res[2], res[1]);
+  }
+#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/test/Surface_mesh_simplification/CMakeLists.txt

@@ -12,6 +12,7 @@ create_single_source_cgal_program("test_edge_collapse_Envelope.cpp")
 create_single_source_cgal_program("test_edge_collapse_Polyhedron_3.cpp")
 create_single_source_cgal_program("test_edge_profile_link.cpp")
 create_single_source_cgal_program("test_edge_deprecated_stop_predicates.cpp")
+create_single_source_cgal_program("test_edge_collapse_stability.cpp")
 
 find_package(Eigen3 3.1.0 QUIET) #(3.1.0 or greater)
 include(CGAL_Eigen3_support)

Surface_mesh_simplification/test/Surface_mesh_simplification/data/far_xy.off

@@ -0,0 +1,90 @@
+OFF
+36 50 0
+
+104605.80000019073 1217329.6000003815 134.27659606933594
+104605.75 1217329.6000003815 134.27859497070312
+104605.55000019073 1217329.7999992371 134.29460144042969
+104605.80000019073 1217329.7999992371 134.28489685058594
+104605.75 1217329.7999992371 134.28680419921875
+104605.75 1217329.75 134.28469848632812
+104605.55000019073 1217329.75 134.2926025390625
+104605.55000019073 1217329.7000007629 134.29060363769531
+104605.80000019073 1217329.75 134.28280639648438
+104605.84999990463 1217329.7999992371 134.28300476074219
+104605.69999980927 1217329.7999992371 134.2886962890625
+104605.55000019073 1217329.6499996185 134.28860473632812
+104605.55000019073 1217329.6000003815 134.28660583496094
+104605.84999990463 1217329.6000003815 134.27470397949219
+104605.59999990463 1217329.7999992371 134.2926025390625
+104605.75 1217329.7000007629 134.2825927734375
+104605.65000009537 1217329.75 134.28860473632812
+104605.80000019073 1217329.7000007629 134.28070068359375
+104605.84999990463 1217329.75 134.28089904785156
+104605.69999980927 1217329.7000007629 134.28460693359375
+104605.90000009537 1217329.7999992371 134.28109741210938
+104605.90000009537 1217329.6499996185 134.27490234375
+104605.90000009537 1217329.7000007629 134.27690124511719
+104605.65000009537 1217329.5499992371 134.28059387207031
+104605.59999990463 1217329.5 134.28050231933594
+104605.90000009537 1217329.5499992371 134.27070617675781
+104605.69999980927 1217329.6000003815 134.28059387207031
+104605.65000009537 1217329.6499996185 134.28460693359375
+104605.59999990463 1217329.7000007629 134.28860473632812
+104605.65000009537 1217329.4500007629 134.27650451660156
+104605.59999990463 1217329.6000003815 134.28460693359375
+104605.55000019073 1217329.5 134.28250122070312
+104605.55000019073 1217329.5499992371 134.28450012207031
+104605.90000009537 1217329.75 134.27900695800781
+104605.84999990463 1217329.7000007629 134.27879333496094
+104605.84999990463 1217329.6499996185 134.27679443359375
+3  30 23 26
+3  3 4 5
+3  26 23 1
+3  22 33 34
+3  8 9 3
+3  17 0 35
+3  17 1 0
+3  4 10 5
+3  15 1 17
+3  8 3 5
+3  15 26 1
+3  6 14 2
+3  0 25 13
+3  35 0 13
+3  6 16 14
+3  28 16 6
+3  7 28 6
+3  8 5 15
+3  16 10 14
+3  8 15 17
+3  16 5 10
+3  18 8 17
+3  16 19 5
+3  19 15 5
+3  18 20 9
+3  18 9 8
+3  11 28 7
+3  11 27 28
+3  21 35 13
+3  21 22 35
+3  22 34 35
+3  31 24 32
+3  24 23 32
+3  24 29 23
+3  29 1 23
+3  19 26 15
+3  27 26 19
+3  28 27 19
+3  28 19 16
+3  29 0 1
+3  30 26 27
+3  29 25 0
+3  18 33 20
+3  34 33 18
+3  34 18 17
+3  34 17 35
+3  32 23 30
+3  32 30 12
+3  12 30 11
+3  30 27 11
+

Surface_mesh_simplification/test/Surface_mesh_simplification/test_edge_collapse_stability.cpp

@@ -0,0 +1,71 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Surface_mesh.h>
+
+// Simplification function
+#include <CGAL/Surface_mesh_simplification/edge_collapse.h>
+#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_count_stop_predicate.h>
+#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_cost.h>
+#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_placement.h>
+
+#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
+
+//bbox
+#include <CGAL/Polygon_mesh_processing/bbox.h>
+
+#include <iostream>
+#include <fstream>
+
+namespace SMS = CGAL::Surface_mesh_simplification;
+
+typedef CGAL::Simple_cartesian<double>                        Kernel;
+
+typedef Kernel::Point_3                                       Point_3;
+typedef CGAL::Surface_mesh<Point_3>                           Surface;
+
+typedef SMS::LindstromTurk_cost<Surface>                      Cost;
+typedef SMS::LindstromTurk_placement<Surface>                 Placement;
+
+int main(int argc, char** argv)
+{
+  std::cout.precision(17);
+  std::cerr.precision(17);
+
+  std::string filename = (argc > 1) ? argv[1] : "data/far_xy.off";
+
+  Surface mesh;
+  if(!CGAL::IO::read_polygon_mesh(filename, mesh))
+  {
+    std::cerr << "Failed to read input mesh: " << filename << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  std::cout << "input has " << num_vertices(mesh) << " vertices." << std::endl;
+
+  CGAL::Iso_cuboid_3<Kernel> bbox(CGAL::Polygon_mesh_processing::bbox(mesh));
+
+  // scale a bit the bounding box bbox, because the kernel is SC
+  bbox = { bbox.min() - 0.01 * (bbox.max() - bbox.min()),
+           bbox.max() + 0.01 * (bbox.max() - bbox.min()) };
+
+  std::cout << "Bbox: " << bbox << std::endl;
+
+  SMS::Edge_count_stop_predicate<Surface> stop(num_halfedges(mesh)/10);
+  Placement placement_ref;
+
+  SMS::edge_collapse(mesh, stop,
+                     CGAL::parameters::get_cost(Cost())
+                                      .get_placement(placement_ref));
+
+  CGAL::IO::write_polygon_mesh("out.off", mesh, CGAL::parameters::stream_precision(17));
+
+  for(auto v : vertices(mesh))
+  {
+    if(bbox.has_on_unbounded_side(mesh.point(v)))
+    {
+      std::cerr << "Error: " << mesh.point(v) << " is outside" << std::endl;
+    }
+  }
+
+  std::cout << "output has " << vertices(mesh).size() << " vertices." << std::endl;
+  return EXIT_SUCCESS;
+}