SMS: Make Lindstrom Turk placement more robust

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/CrossProduct.h

@@ -0,0 +1,61 @@
+// 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))
   {
+#if 1
+    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;
+#else
     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;
@@ -230,6 +244,7 @@ std::optional< MatrixC33<R> > inverse_matrix(const MatrixC33<R>& m)
     RT c20 = (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 c22 = (m.r0().x()*m.r1().y() - m.r0().y()*m.r1().x()) / det;
+#endif
 
     rInverse = result_type(Matrix(c00,c01,c02,
                                   c10,c11,c12,

Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounding_box_filter.h

@@ -0,0 +1,77 @@
+// Copyright (c) 2024  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)     : Andreas Fabri
+//
+#ifndef CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_BOUNDING_BOX_FILTER_H
+#define CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_BOUNDING_BOX_FILTER_H
+
+#include <CGAL/license/Surface_mesh_simplification.h>
+#include<CGAL/Surface_mesh_simplification/internal/Common.h>
+#include <CGAL/property_map.h>
+
+#include <optional>
+
+namespace CGAL {
+namespace Surface_mesh_simplification {
+
+
+template<class BaseFilter = internal::Dummy_filter>
+class Bounding_box_filter
+{
+public:
+  Bounding_box_filter(const BaseFilter& base_filter = BaseFilter())
+    : m_base_filter(base_filter)
+  {}
+
+
+  template <typename Profile>
+  std::optional<typename Profile::Point>
+  operator()(const Profile& profile, std::optional<typename Profile::Point> op) const
+  {
+    typedef typename Profile::VertexPointMap                              Vertex_point_map;
+
+    typedef typename Profile::Geom_traits                                 Geom_traits;
+    typedef typename Geom_traits::Vector_3                                Vector;
+
+    typedef typename boost::property_traits<Vertex_point_map>::value_type Point;
+    typedef typename boost::property_traits<Vertex_point_map>::reference  Point_reference;
+
+    const Geom_traits& gt = profile.geom_traits();
+    const Vertex_point_map& vpm = profile.vertex_point_map();
+
+    op = m_base_filter(profile, op);
+    if(op)
+    {
+      const Point& placement = *op;
+      Bbox_3 bb;
+      for(auto vd : profile.link()){
+        Point_reference p = get(vpm, vd);
+        bb += p.bbox();
+      }
+      double wx = bb.xmax() - bb.xmin();
+      double wy = bb.ymax() - bb.ymin();
+      double wz = bb.zmax() - bb.zmin();
+      bb = Bbox_3(bb.xmin()-  wx/10.0, bb.ymin() - wy/10.0, bb.zmin()- wz/10.0, bb.xmax() + wx/10.0, bb.ymax() + wy/10.0, bb.zmax()+ wz/10.0);
+      if(!do_overlap(bb, placement.bbox())){
+             return std::optional<typename Profile::Point>();
+      }
+    }
+    return op;
+  }
+
+
+
+private:
+  const BaseFilter m_base_filter;
+};
+
+} // namespace Surface_mesh_simplification
+} // namespace CGAL
+
+#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_BOUNDING_BOX_FILTER_H