Fix of issue 8213 by consider length of a vector being zero if enough small, define diff_of_product in a specific file

Surface_mesh_simplification/include/CGAL/Cartesian/CrossProduct.h

@@ -0,0 +1,60 @@
+// Copyright (c) 2025  GeometryFactory (France). All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//  Author(s)     : Mael Rouxel-Labbé
+//
+#ifndef CGAL_CARTESIAN_CROSSPRODUCT_H
+#define CGAL_CARTESIAN_CROSSPRODUCT_H
+
+#include <CGAL/license/Surface_mesh_simplification.h>
+
+#include <optional>
+
+namespace CGAL {
+
+// 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)
+  {
+    // return a*b - c*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;
+  }
+
+} // namespace CGAL
+
+#endif // CGAL_CARTESIAN_CROSSPRODUCT_H //
+// EOF //
+
+

Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h

@@ -17,6 +17,7 @@
 #include <CGAL/Null_matrix.h>
 #include <CGAL/number_utils.h>
 #include <CGAL/Vector_3.h>
+#include <CGAL/Cartesian/CrossProduct.h>
 
 #include <optional>
 
@@ -219,6 +220,19 @@ std::optional< MatrixC33<R> > inverse_matrix(const MatrixC33<R>& m)
 
   if(! CGAL_NTS is_zero(det))
   {
+
+    // RT c00 = diff_of_products(m.r1().y(),m.r2().z(),m.r2().y(),m.r1().z()) / det;
+    // RT c01 = diff_of_products(m.r2().y(),m.r0().z(),m.r0().y(),m.r2().z()) / det;
+    // RT c02 = diff_of_products(m.r0().y(),m.r1().z(),m.r1().y(),m.r0().z()) / det;
+
+    // RT c10 = diff_of_products(m.r2().x(),m.r1().z(),m.r1().x(),m.r2().z()) / det;
+    // RT c11 = diff_of_products(m.r0().x(),m.r2().z(),m.r2().x(),m.r0().z()) / det;
+    // RT c12 = diff_of_products(m.r1().x(),m.r0().z(),m.r0().x(),m.r1().z()) / det;
+
+    // RT c20 = diff_of_products(m.r1().x(),m.r2().y(),m.r2().x(),m.r1().y()) / det;
+    // RT c21 = diff_of_products(m.r2().x(),m.r0().y(),m.r0().x(),m.r2().y()) / det;
+    // RT c22 = diff_of_products(m.r0().x(),m.r1().y(),m.r1().x(),m.r0().y()) / det;
+
     RT c00 = (m.r1().y()*m.r2().z() - m.r1().z()*m.r2().y()) / det;
     RT c01 = (m.r2().y()*m.r0().z() - m.r0().y()*m.r2().z()) / det;
     RT c02 = (m.r0().y()*m.r1().z() - m.r1().y()*m.r0().z()) / det;

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

@@ -18,6 +18,7 @@
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_profile.h>
 
 #include <CGAL/Cartesian/MatrixC33.h>
+#include <CGAL/Cartesian/CrossProduct.h>
 
 #include <limits>
 #include <vector>
@@ -101,6 +102,8 @@ private :
   void extract_triangle_data();
   void extract_boundary_data();
 
+  double maxBb;
+
   void add_boundary_preservation_constraints(const Boundary_data_vector& aBdry);
   void add_volume_preservation_constraints(const Triangle_data_vector& triangles);
   void add_boundary_and_volume_optimization_constraints(const Boundary_data_vector& aBdry,
@@ -119,42 +122,6 @@ 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)
   {
@@ -353,12 +320,18 @@ extract_triangle_data()
 {
   mTriangle_data.reserve(mProfile.triangles().size());
 
+  //TODO for obscur reason, computing this abs_max increase running time by 10%
+  double abs_max;
   for(const Triangle& tri : mProfile.triangles())
   {
     const Point_reference p0 = get_point(tri.v0);
     const Point_reference p1 = get_point(tri.v1);
     const Point_reference p2 = get_point(tri.v2);
 
+    abs_max=(std::max)({abs_max,std::abs(p0.x()),std::abs(p0.y()),std::abs(p0.z()),
+                             std::abs(p1.x()),std::abs(p1.y()),std::abs(p1.z()),
+                             std::abs(p2.x()),std::abs(p2.y()),std::abs(p2.z())});
+
     Vector v01 = p1 - p0;
     Vector v02 = p2 - p0;
 
@@ -370,6 +343,7 @@ extract_triangle_data()
 
     mTriangle_data.push_back(Triangle_data(lNormalV,lNormalL));
   }
+  maxBb= abs_max;
 }
 
 template<class TM, class K>
@@ -687,7 +661,11 @@ add_constraint_if_alpha_compatible(const Vector& Ai,
     FT l = CGAL_NTS sqrt(slai);
     CGAL_SMS_LT_TRACE(3, "      l: " << n_to_string(l));
 
-    if(!CGAL_NTS is_zero(l))
+    // Due to double number type, l may have a small value instead of zero (example sum of the faces normals of a tetrahedra for volumic constraint)
+    // if bin is greater than maxBb, we consider that l is zero
+    CGAL_SMS_LT_TRACE(3, "      error consider: " << (std::abs(bi) / (2*maxBb)));
+    if(l > (std::abs(bi) / (2*maxBb)))
+    // if(!CGAL_NTS is_zero(l))
     {
       Vector Ain = Ai / l;
       FT bin = bi / l;
@@ -862,6 +840,7 @@ add_constraint_from_gradient(const Matrix& H,
     }
       break;
   }
+
 }
 
 } // namespace internal

Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp

@@ -3,6 +3,8 @@
 
 #include <CGAL/Surface_mesh.h>
 
+CGAL::Bbox_3 bbox_g;
+double max_bbox_g;
 #define CGAL_CHECK_EXPENSIVE
 
 #define CGAL_SURFACE_SIMPLIFICATION_ENABLE_TRACE 5
@@ -15,7 +17,7 @@ void Surface_simplification_external_trace(const std::string& s)
 
 #include <CGAL/Surface_mesh_simplification/edge_collapse.h>
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounded_normal_change_filter.h>
-#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounding_box_filter.h>
+// #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounding_box_filter.h>
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Face_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>
@@ -48,19 +50,19 @@ int main(int argc, char** argv)
     return EXIT_FAILURE;
   }
 
-  CGAL::Bbox_3 bbox = PMP::bbox(sm);
-
+  auto bb=PMP::bbox(sm);
+  std::cout << "Bbox:" << bb << std::endl;
   std::cout << "Input mesh has " << num_vertices(sm) << " vertices" << std::endl;
   std::cout << "Input mesh has " << num_faces(sm) << " faces" << std::endl;
 
 
   SMS::Face_count_stop_predicate<Surface_mesh> stop(1);
    SMS::edge_collapse(
-    surface_mesh,
+    sm,
     stop,
     CGAL::parameters::
-       filter(SMS::Bounded_normal_change_filter<SMS::Bounding_box_filter<>>(SMS::Bounding_box_filter<>()))
-      .get_cost(SMS::LindstromTurk_cost<Surface_mesh>())
+      //  filter(SMS::Bounded_normal_change_filter<>())  //<SMS::Bounding_box_filter<>>(SMS::Bounding_box_filter<>()))
+      get_cost(SMS::LindstromTurk_cost<Surface_mesh>())
       .get_placement(SMS::LindstromTurk_placement<Surface_mesh>())
     );
   CGAL::IO::write_OFF(std::cout, sm, CGAL::parameters::stream_precision(17));
@@ -69,7 +71,7 @@ int main(int argc, char** argv)
   {
     // To be within the bounding box isn't a guarantee, but here it is a sufficient test
     // to check if things went awry
-    if(!CGAL::do_overlap(bbox, sm.point(v).bbox()))
+    if(!CGAL::do_overlap(bb, sm.point(v).bbox()))
     {
       std::cerr << "Error: " << sm.point(v) << " is outside the initial bbox" << std::endl;
       return EXIT_FAILURE;