SMS: Make Lindstrom Turk placement more robust (#8237)

## Summary of Changes

Reduced data set to illustrate Issue #8213. 
The collapse of edge `v1-v2` results in a completely wrong placement of
the vertex.


![image](https://github.com/CGAL/cgal/assets/3263539/97ff96a1-7f6b-4182-b7b8-5267cb53f8ef)

The next step is to make the placement robust.

## Release Management

* Affected package(s): Surface_mesh_simplification.
* Issue(s) solved (if any): fix #8213
* License and copyright ownership:   unchanged

Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_bounded_normal_change.cpp

@@ -1,5 +1,6 @@
 #include <CGAL/Simple_cartesian.h>
 #include <CGAL/Surface_mesh.h>
+#include <CGAL/IO/polygon_mesh_io.h>
 #include <CGAL/Timer.h>
 
 #include <CGAL/Surface_mesh_simplification/edge_collapse.h>
@@ -36,8 +37,8 @@ int main(int argc, char** argv)
 {
   Surface_mesh surface_mesh;
   const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/fold.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;

Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h

@@ -13,6 +13,8 @@
 
 #include <CGAL/license/Surface_mesh_simplification.h>
 
+#include <CGAL/internal/robust_cross_product.h>
+
 #include <CGAL/determinant.h>
 #include <CGAL/Null_matrix.h>
 #include <CGAL/number_utils.h>
@@ -183,17 +185,19 @@ MatrixC33<R> cofactors_matrix(const MatrixC33<R>& m)
 {
   typedef typename R::RT RT;
 
-  RT c00 =  determinant(m.r1().y(),m.r1().z(),m.r2().y(),m.r2().z());
+  using ::CGAL::Surface_mesh_simplification::internal::diff_of_products;
-  RT c01 = -determinant(m.r1().x(),m.r1().z(),m.r2().x(),m.r2().z());
-  RT c02 =  determinant(m.r1().x(),m.r1().y(),m.r2().x(),m.r2().y());
 
-  RT c10 = -determinant(m.r0().y(),m.r0().z(),m.r2().y(),m.r2().z());
+  RT c00 =  diff_of_products(m.r1().y(), m.r2().z(), m.r2().y(), m.r1().z());
-  RT c11 =  determinant(m.r0().x(),m.r0().z(),m.r2().x(),m.r2().z());
+  RT c01 = -diff_of_products(m.r1().x(), m.r2().z(), m.r2().x(), m.r1().z());
-  RT c12 = -determinant(m.r0().x(),m.r0().y(),m.r2().x(),m.r2().y());
+  RT c02 =  diff_of_products(m.r1().x(), m.r2().y(), m.r2().x(), m.r1().y());
 
-  RT c20 =  determinant(m.r0().y(),m.r0().z(),m.r1().y(),m.r1().z());
+  RT c10 = -diff_of_products(m.r0().y(), m.r2().z(), m.r2().y(), m.r0().z());
-  RT c21 = -determinant(m.r0().x(),m.r0().z(),m.r1().x(),m.r1().z());
+  RT c11 =  diff_of_products(m.r0().x(), m.r2().z(), m.r2().x(), m.r0().z());
-  RT c22 =  determinant(m.r0().x(),m.r0().y(),m.r1().x(),m.r1().y());
+  RT c12 = -diff_of_products(m.r0().x(), m.r2().y(), m.r2().x(), m.r0().y());
+
+  RT c20 =  diff_of_products(m.r0().y(), m.r1().z(), m.r1().y(), m.r0().z());
+  RT c21 = -diff_of_products(m.r0().x(), m.r1().z(), m.r1().x(), m.r0().z());
+  RT c22 =  diff_of_products(m.r0().x(), m.r1().y(), m.r1().x(), m.r0().y());
 
   return MatrixC33<R>(c00,c01,c02,
                       c10,c11,c12,
@@ -213,23 +217,25 @@ std::optional< MatrixC33<R> > inverse_matrix(const MatrixC33<R>& m)
   typedef MatrixC33<R>                                          Matrix;
   typedef std::optional<Matrix>                                 result_type;
 
+  using ::CGAL::Surface_mesh_simplification::internal::diff_of_products;
+
   result_type rInverse;
 
   RT det = m.determinant();
 
   if(! CGAL_NTS is_zero(det))
   {
-    RT c00 = (m.r1().y()*m.r2().z() - m.r1().z()*m.r2().y()) / det;
+    RT c00 = diff_of_products(m.r1().y(),m.r2().z(),m.r2().y(),m.r1().z()) / det;
-    RT c01 = (m.r2().y()*m.r0().z() - m.r0().y()*m.r2().z()) / det;
+    RT c01 = diff_of_products(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;
+    RT c02 = diff_of_products(m.r0().y(),m.r1().z(),m.r1().y(),m.r0().z()) / det;
 
-    RT c10 = (m.r1().z()*m.r2().x() - m.r1().x()*m.r2().z()) / det;
+    RT c10 = diff_of_products(m.r2().x(),m.r1().z(),m.r1().x(),m.r2().z()) / det;
-    RT c11 = (m.r0().x()*m.r2().z() - m.r2().x()*m.r0().z()) / det;
+    RT c11 = diff_of_products(m.r0().x(),m.r2().z(),m.r2().x(),m.r0().z()) / det;
-    RT c12 = (m.r1().x()*m.r0().z() - m.r0().x()*m.r1().z()) / det;
+    RT c12 = diff_of_products(m.r1().x(),m.r0().z(),m.r0().x(),m.r1().z()) / det;
 
-    RT c20 = (m.r1().x()*m.r2().y() - m.r2().x()*m.r1().y()) / det;
+    RT c20 = diff_of_products(m.r1().x(),m.r2().y(),m.r2().x(),m.r1().y()) / det;
-    RT c21 = (m.r2().x()*m.r0().y() - m.r0().x()*m.r2().y()) / det;
+    RT c21 = diff_of_products(m.r2().x(),m.r0().y(),m.r0().x(),m.r2().y()) / det;
-    RT c22 = (m.r0().x()*m.r1().y() - m.r0().y()*m.r1().x()) / det;
+    RT c22 = diff_of_products(m.r0().x(),m.r1().y(),m.r1().x(),m.r0().y()) / det;
 
     rInverse = result_type(Matrix(c00,c01,c02,
                                   c10,c11,c12,

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/internal/robust_cross_product.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)
   {
@@ -185,67 +152,10 @@ private :
   }
 #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 compute_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 = compute_minor(ay, by, az, bz);
-    FT y = compute_minor(az, bz, ax, bx);
-    FT z = compute_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 X_product(a-ORIGIN, b-ORIGIN);
+    return robust_cross_product<Geom_traits>(a-ORIGIN, b-ORIGIN);
   }
 
   // This is the (uX)(Xu) product described in the Lindstrom-Turk paper
@@ -352,17 +262,21 @@ LindstromTurkCore<TM,K>::
 extract_triangle_data()
 {
   mTriangle_data.reserve(mProfile.triangles().size());
-
+  maxBb = 0.0;
   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);
 
+    maxBb=(std::max)({maxBb,CGAL::abs(p0.x()),CGAL::abs(p0.y()),CGAL::abs(p0.z()),
+                            CGAL::abs(p1.x()),CGAL::abs(p1.y()),CGAL::abs(p1.z()),
+                            CGAL::abs(p2.x()),CGAL::abs(p2.y()),CGAL::abs(p2.z())});
+
     Vector v01 = p1 - p0;
     Vector v02 = p2 - p0;
 
-    Vector lNormalV = cross_product(v01,v02);
+    Vector lNormalV = robust_cross_product<Geom_traits>(v01,v02);
     FT lNormalL = point_cross_product(p0,p1) * (p2 - ORIGIN);
 
     CGAL_SMS_LT_TRACE(1, "  Extracting triangle v" << tri.v0 << "->v" << tri.v1 << "->v" << tri.v2
@@ -370,6 +284,7 @@ extract_triangle_data()
 
     mTriangle_data.push_back(Triangle_data(lNormalV,lNormalL));
   }
+  maxBb *= 2.0; // to avoid numerical problems
 }
 
 template<class TM, class K>
@@ -394,11 +309,12 @@ compute_placement()
   //  A1 * v = b1
   //  A2 * v = b2
   //
-  // Which in matrix form is :  A * v = b
+  // Which in matrix form is:  A * v = b
   //
   // (with 'A' a 3x3 matrix and 'b' a vector)
   //
-  // The member variable mConstrinas contains A and b. Indidivual constraints (Ai,bi) can be added to it.
+  // The member variables mConstraints_A and mConstraints_b contain A and b.
+  // Indidivual constraints (Ai,bi) can be added to it.
   // Once 3 such constraints have been added 'v' is directly solved a:
   //
   //  v = b*inverse(A)
@@ -421,7 +337,7 @@ compute_placement()
   // In that case there is simply no good vertex placement
   if(mConstraints_n == 3)
   {
-    // If the matrix is singular it's inverse cannot be computed so an 'absent' value is returned.
+    // If the matrix is singular its inverse cannot be computed so an 'absent' value is returned.
     std::optional<Matrix> lOptional_Ai = inverse_matrix(mConstraints_A);
     if(lOptional_Ai)
     {
@@ -686,7 +602,10 @@ 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 face normals of a tetrahedron for volume constraint)
+    // if bi is greater than maxBb, we consider that l is zero
+    CGAL_SMS_LT_TRACE(3, "      error consider: " << (CGAL::abs(bi) / maxBb));
+    if(l > (CGAL::abs(bi) / maxBb))
     {
       Vector Ain = Ai / l;
       FT bin = bi / l;
@@ -861,6 +780,7 @@ add_constraint_from_gradient(const Matrix& H,
     }
       break;
   }
+
 }
 
 } // namespace internal

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

@@ -114,8 +114,9 @@ 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::ostringstream ss; ss << m; std::string s = ss.str(); \
+      std::ostringstream ss; ss << m; \
-    /*Surface_simplification_external_trace(s)*/ std::cerr << s << std::endl; \
+      std::string s = ss.str(); \
+      Surface_simplification_external_trace(s); \
     }
 
   #define CGAL_SMS_DEBUG_CODE(code) code

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

@@ -296,6 +296,8 @@ private:
 
   void insert_in_PQ(const halfedge_descriptor h, Edge_data& data)
   {
+    CGAL_SMS_TRACE(5, "Insert " << edge_to_string(h) << " in PQ");
+
     CGAL_assertion(is_primary_edge(h));
     CGAL_expensive_assertion(!data.is_in_PQ());
     CGAL_expensive_assertion(!mPQ->contains(h));
@@ -594,12 +596,33 @@ loop()
 
   std::optional<halfedge_descriptor> opt_h;
 
+// #define CGAL_SURF_SIMPL_INTERMEDIATE_STEPS_PRINTING
 #ifdef CGAL_SURF_SIMPL_INTERMEDIATE_STEPS_PRINTING
   int i_rm = 0;
 #endif
 
-  while((opt_h = pop_from_PQ()))
+  for(;;)
   {
+#ifdef CGAL_SURFACE_SIMPLIFICATION_ENABLE_TRACE
+    if(5 <= CGAL_SURFACE_SIMPLIFICATION_ENABLE_TRACE)
+    {
+      CGAL_SMS_TRACE_IMPL("== Current queue ==");
+
+      auto mPQ_clone = *mPQ;
+      std::optional<halfedge_descriptor> opt_th;
+      while((opt_th = mPQ_clone.extract_top()))
+      {
+        CGAL_SMS_TRACE_IMPL("\t" + edge_to_string(*opt_th));
+        Cost_type tcost = get_data(*opt_th).cost();
+        if(tcost)
+          CGAL_SMS_TRACE_IMPL("\t" + std::to_string(CGAL::to_double(*tcost)));
+      }
+    }
+#endif
+
+    if(!(opt_h = pop_from_PQ()))
+      break;
+
     CGAL_SMS_TRACE(1, "Popped " << edge_to_string(*opt_h));
     CGAL_assertion(!is_constrained(*opt_h));
 
@@ -639,7 +662,7 @@ loop()
 
             m_visitor.OnNonCollapsable(profile);
 
-            CGAL_SMS_TRACE(1, edge_to_string(*opt_h) << " NOT Collapsible" );
+            CGAL_SMS_TRACE(1, edge_to_string(*opt_h) << " NOT Collapsible (filter)" );
           }
 
 #ifdef CGAL_SURF_SIMPL_INTERMEDIATE_STEPS_PRINTING
@@ -660,7 +683,7 @@ loop()
 
         m_visitor.OnNonCollapsable(profile);
 
-        CGAL_SMS_TRACE(1, edge_to_string(*opt_h) << " NOT Collapsible" );
+        CGAL_SMS_TRACE(1, edge_to_string(*opt_h) << " NOT Collapsible (topology)" );
       }
     }
     else
@@ -823,7 +846,10 @@ is_collapse_topologically_valid(const Profile& profile)
   {
     /// ensure two constrained edges cannot get merged
     if(is_edge_adjacent_to_a_constrained_edge(profile, m_ecm))
+    {
+      CGAL_SMS_TRACE(3,"  edge to collapse is adjacent to a constrained edge.");
       return false;
+    }
 
     if(profile.is_v0_v1_a_border())
     {

Surface_mesh_simplification/include/CGAL/internal/robust_cross_product.h

@@ -0,0 +1,143 @@
+// 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_SMS_INTERNAL_ROBUST_CROSS_PRODUCT_H
+#define CGAL_SMS_INTERNAL_ROBUST_CROSS_PRODUCT_H
+
+#include <CGAL/license/Surface_mesh_simplification.h>
+
+#include <CGAL/Surface_mesh_simplification/internal/Common.h>
+
+#include <optional>
+
+namespace CGAL {
+namespace Surface_mesh_simplification {
+namespace internal {
+
+// 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
+inline 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;
+}
+
+inline 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;
+}
+
+inline double diff_of_products(const double a, const double b, const double c, const double d)
+{
+#if 0
+  // this can create large errors with inexact constructions
+  return a*b - c*d;
+  // the next two are equivalent in results and speed
+#elif 1
+  return diff_of_products_kahan(a, b, c, d);
+#elif 0
+  return diff_of_products_cht(a, b, c, d);
+#endif
+}
+
+template <typename OFT>
+inline OFT diff_of_products(const OFT& a, const OFT& b, const OFT& c, const OFT& d)
+{
+  return a*b - c*d;
+}
+
+// balanced solution based on abusing the fact that here we expect u and v to have similar coordinates
+template<class GeomTraits>
+typename GeomTraits::Vector_3 similar_coordinates_cross_product(const typename GeomTraits::Vector_3& u,
+                                                                const typename GeomTraits::Vector_3& v)
+{
+  using FT = typename GeomTraits::FT;
+  using Vector = typename GeomTraits::Vector_3;
+
+  const FT& ux = u.x();
+  const FT& uy = u.y();
+  const FT& uz = u.z();
+  const FT& vx = v.x();
+  const FT& vy = v.y();
+  const FT& vz = v.z();
+
+  auto minor = [](const FT& ui, const FT& vi, const FT& uj, const FT& vj)
+  {
+    // 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 ui * (vj - uj) + uj * (ui - vi);
+  };
+
+  // ay*
+  FT x = minor(uy, vy, uz, vz);
+  FT y = minor(uz, vz, ux, vx);
+  FT z = minor(ux, vx, uy, vy);
+
+  return Vector(x, y, z);
+}
+
+#if 0
+template<class GeomTraits>
+typename GeomTraits::Vector_3 exact_cross_product(const typename GeomTraits::Vector_3& a,
+                                                  const typename GeomTraits::Vector_3& b)
+{
+  CGAL::Cartesian_converter<GeomTraits, CGAL::Exact_predicates_exact_constructions_kernel> to_exact;
+  CGAL::Cartesian_converter<CGAL::Exact_predicates_exact_constructions_kernel, GeomTraits> to_approx;
+  auto exv = cross_product(to_exact(a), to_exact(b));
+  exv.exact();
+  return to_approx(exv);
+}
+#endif
+
+template<class GeomTraits>
+typename GeomTraits::Vector_3 robust_cross_product(const typename GeomTraits::Vector_3& u,
+                                                   const typename GeomTraits::Vector_3& 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 0
+  // obviously too slow
+  return exact_cross_product(u, v);
+#elif 1
+  // balanced solution based on abusing the fact that in this package, we usually have that
+  // u and v to have similar coordinates
+  return similar_coordinates_cross_product<GeomTraits>(u, v);
+#endif
+}
+
+} // namespace CGAL
+} // namespace Surface_mesh_simplification
+} // namespace internal
+
+#endif // CGAL_SMS_INTERNAL_ROBUST_CROSS_PRODUCT_H
+
+

Surface_mesh_simplification/test/Surface_mesh_simplification/CMakeLists.txt

@@ -6,6 +6,7 @@ project(Surface_mesh_simplification_Tests)
 
 find_package(CGAL REQUIRED)
 
+create_single_source_cgal_program("issue_8213.cpp")
 create_single_source_cgal_program("edge_collapse_topology.cpp")
 create_single_source_cgal_program("test_edge_collapse_bounded_distance.cpp")
 create_single_source_cgal_program("test_edge_collapse_Envelope.cpp")

Surface_mesh_simplification/test/Surface_mesh_simplification/data/issue_8213.off

@@ -0,0 +1,13 @@
+OFF
+5 4 0
+
+0.60153250345565623 3.2925554343503594 -0.93390733763467004
+0.50125687092531912 3.266008536541555 -0.80580753798383942
+0.57499779785916183 3.2558452065056969 -0.97860403852322797
+0.56586410588624558 3.2541065339825863 -0.99341202997519495
+0.56756366821062887 3.2478315549358072 -0.99100621040927039
+
+3    0 1 2
+3    1 3 2
+3    1 0 3
+3    0 4 3

Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp

@@ -0,0 +1,75 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+#include <CGAL/Surface_mesh.h>
+
+#define CGAL_CHECK_EXPENSIVE
+
+#define CGAL_SURFACE_SIMPLIFICATION_ENABLE_TRACE 5
+#define CGAL_SURFACE_SIMPLIFICATION_ENABLE_LT_TRACE 4
+
+void Surface_simplification_external_trace(const std::string& s)
+{
+  std::cout << s << std::endl;
+}
+
+#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/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>
+
+#include <CGAL/Polygon_mesh_processing/bbox.h>
+
+#include <cstdlib>
+#include <iostream>
+#include <fstream>
+
+namespace SMS = CGAL::Surface_mesh_simplification;
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+using Kernel = CGAL::Simple_cartesian<double>;
+// using Kernel = CGAL::Exact_predicates_exact_constructions_kernel;
+
+using Surface_mesh = CGAL::Surface_mesh<Kernel::Point_3>;
+
+int main(int argc, char** argv)
+{
+  std::cout.precision(17);
+  std::cerr.precision(17);
+
+  const char* filename = (argc > 1) ? argv[1] : "data/issue_8213.off";
+
+  Surface_mesh sm;
+
+  if(!CGAL::IO::read_polygon_mesh(filename, sm))
+  {
+    std::cerr << "Error: failed to read input data" << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  CGAL::Bbox_3 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(sm, stop,
+                      CGAL::parameters::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));
+
+  for(auto v : vertices(sm))
+  {
+    // 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(bb, sm.point(v).bbox()))
+    {
+      std::cerr << "Error: " << sm.point(v) << " is outside the initial bbox" << std::endl;
+      return EXIT_FAILURE;
+    }
+  }
+
+  return EXIT_SUCCESS;
+}