From 8e5c5b84f30115e31983034852bd44ea5b2496d4 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 29 May 2024 09:52:23 +0100 Subject: [PATCH 01/19] SMS: Fix placement --- .../CMakeLists.txt | 1 + .../data/issue_8213.off | 16 +++++++ .../issue_8213.cpp | 44 +++++++++++++++++++ 3 files changed, 61 insertions(+) create mode 100644 Surface_mesh_simplification/test/Surface_mesh_simplification/data/issue_8213.off create mode 100644 Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp diff --git a/Surface_mesh_simplification/test/Surface_mesh_simplification/CMakeLists.txt b/Surface_mesh_simplification/test/Surface_mesh_simplification/CMakeLists.txt index 9c1e010e5ef..dabc2c49c08 100644 --- a/Surface_mesh_simplification/test/Surface_mesh_simplification/CMakeLists.txt +++ b/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") diff --git a/Surface_mesh_simplification/test/Surface_mesh_simplification/data/issue_8213.off b/Surface_mesh_simplification/test/Surface_mesh_simplification/data/issue_8213.off new file mode 100644 index 00000000000..4669828bf08 --- /dev/null +++ b/Surface_mesh_simplification/test/Surface_mesh_simplification/data/issue_8213.off @@ -0,0 +1,16 @@ +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 + + + diff --git a/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp b/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp new file mode 100644 index 00000000000..be6b26d0fb7 --- /dev/null +++ b/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp @@ -0,0 +1,44 @@ +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace SMS = CGAL::Surface_mesh_simplification; + +typedef CGAL::Simple_cartesian Kernel; +typedef CGAL::Surface_mesh Surface_mesh; + +int main() { + + std::string filename("data/issue_8213.off"); + + Surface_mesh surface_mesh; + + std::ifstream in(filename); + in >> surface_mesh; + + const size_t target_number_of_faces = 3; + SMS::Face_count_stop_predicate stop(target_number_of_faces); + + std::cout << "Input mesh number of faces: " << surface_mesh.number_of_faces() << ", target number of faces: " << target_number_of_faces << std::endl; + + SMS::edge_collapse( + surface_mesh, + stop, + CGAL::parameters:: + filter(SMS::Bounded_normal_change_filter<>()) + .get_cost(SMS::LindstromTurk_cost()) + .get_placement(SMS::LindstromTurk_placement()) + ); + + std::cout.precision(17); + std::cout << surface_mesh << std::endl; + return EXIT_SUCCESS; +} From 99613eb99c3bd0da0884916b1aef5e6af34badc8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 30 May 2024 14:45:40 +0200 Subject: [PATCH 02/19] Generalize an expression to make it easier to use other kernels --- .../Policies/Edge_collapse/internal/Lindstrom_Turk_core.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h index d97688dde8b..21f4018c0e3 100644 --- a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h +++ b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h @@ -190,7 +190,7 @@ private : 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) + auto minor = [](auto ai, auto bi, auto aj, auto 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 From b80c1d8f482a38db3e0a7ee0d40668c9e4552455 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 30 May 2024 14:46:19 +0200 Subject: [PATCH 03/19] Comments & debug code --- .../internal/Lindstrom_Turk_core.h | 7 ++-- .../internal/Edge_collapse.h | 32 +++++++++++++++++-- 2 files changed, 33 insertions(+), 6 deletions(-) diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h index 21f4018c0e3..44691800463 100644 --- a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h +++ b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h @@ -394,11 +394,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 +422,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 lOptional_Ai = inverse_matrix(mConstraints_A); if(lOptional_Ai) { diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/internal/Edge_collapse.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/internal/Edge_collapse.h index 5c23e83e94f..5d43f1d73cf 100644 --- a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/internal/Edge_collapse.h +++ b/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 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 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()) { From 24e1c96f62310164706bcec06e8520ebaf1a7b30 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 30 May 2024 14:46:30 +0200 Subject: [PATCH 04/19] Fix trace function --- .../CGAL/Surface_mesh_simplification/internal/Common.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/internal/Common.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/internal/Common.h index 4ed45a02d1e..a44f13d7770 100644 --- a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/internal/Common.h +++ b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/internal/Common.h @@ -115,8 +115,9 @@ inline std::string optional_to_string(const std::optional& 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(); \ - /*Surface_simplification_external_trace(s)*/ std::cerr << s << std::endl; \ + std::ostringstream ss; ss << m; \ + std::string s = ss.str(); \ + Surface_simplification_external_trace(s); \ } #define CGAL_SMS_DEBUG_CODE(code) code From c40da7044175e1d52053ce3f13c407078cce46b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 30 May 2024 14:46:49 +0200 Subject: [PATCH 05/19] Enhance example --- .../issue_8213.cpp | 83 +++++++++++++------ 1 file changed, 57 insertions(+), 26 deletions(-) diff --git a/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp b/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp index be6b26d0fb7..ca4da904544 100644 --- a/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp +++ b/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp @@ -1,44 +1,75 @@ #include +#include + #include + +#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 #include #include #include #include +#include + #include #include #include namespace SMS = CGAL::Surface_mesh_simplification; +namespace PMP = CGAL::Polygon_mesh_processing; -typedef CGAL::Simple_cartesian Kernel; -typedef CGAL::Surface_mesh Surface_mesh; - -int main() { - - std::string filename("data/issue_8213.off"); - - Surface_mesh surface_mesh; - - std::ifstream in(filename); - in >> surface_mesh; - - const size_t target_number_of_faces = 3; - SMS::Face_count_stop_predicate stop(target_number_of_faces); - - std::cout << "Input mesh number of faces: " << surface_mesh.number_of_faces() << ", target number of faces: " << target_number_of_faces << std::endl; - - SMS::edge_collapse( - surface_mesh, - stop, - CGAL::parameters:: - filter(SMS::Bounded_normal_change_filter<>()) - .get_cost(SMS::LindstromTurk_cost()) - .get_placement(SMS::LindstromTurk_placement()) - ); +using Kernel = CGAL::Simple_cartesian; +// using Kernel = CGAL::Exact_predicates_exact_constructions_kernel; +using Surface_mesh = CGAL::Surface_mesh; +int main(int argc, char** argv) +{ std::cout.precision(17); - std::cout << surface_mesh << std::endl; + 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 bbox = PMP::bbox(sm); + + 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 stop(1); + SMS::edge_collapse(sm, stop, + CGAL::parameters::filter(SMS::Bounded_normal_change_filter<>()) + .get_cost(SMS::LindstromTurk_cost()) + .get_placement(SMS::LindstromTurk_placement())); + + 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(bbox, sm.point(v).bbox())) + { + std::cerr << "Error: " << sm.point(v) << " is outside the initial bbox" << std::endl; + return EXIT_FAILURE; + } + } + return EXIT_SUCCESS; } From 042cb277fb5514226bb48cf06f335723ef4134d6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 30 May 2024 14:58:00 +0200 Subject: [PATCH 06/19] Use the more robust formulation on another eligible cross product --- .../Policies/Edge_collapse/internal/Lindstrom_Turk_core.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h index 44691800463..8c04b537e7f 100644 --- a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h +++ b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h @@ -362,7 +362,7 @@ extract_triangle_data() Vector v01 = p1 - p0; Vector v02 = p2 - p0; - Vector lNormalV = cross_product(v01,v02); + Vector lNormalV = X_product(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 From b54a2eab1a056451a3ae6545d2824163623822e1 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Thu, 30 May 2024 14:59:22 +0100 Subject: [PATCH 07/19] Add a filter --- .../Edge_collapse/Bounding_box_filter.h | 77 +++++++++++++++++++ .../issue_8213.cpp | 15 ++-- 2 files changed, 87 insertions(+), 5 deletions(-) create mode 100644 Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounding_box_filter.h diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounding_box_filter.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounding_box_filter.h new file mode 100644 index 00000000000..c8b534f0ee1 --- /dev/null +++ b/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 +#include +#include + +#include + +namespace CGAL { +namespace Surface_mesh_simplification { + + +template +class Bounding_box_filter +{ +public: + Bounding_box_filter(const BaseFilter& base_filter = BaseFilter()) + : m_base_filter(base_filter) + {} + + + template + std::optional + operator()(const Profile& profile, std::optional 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::value_type Point; + typedef typename boost::property_traits::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(); + } + } + 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 diff --git a/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp b/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp index ca4da904544..b086fe22e73 100644 --- a/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp +++ b/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp @@ -15,6 +15,7 @@ void Surface_simplification_external_trace(const std::string& s) #include #include +#include #include #include #include @@ -52,12 +53,16 @@ int main(int argc, char** argv) 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 stop(1); - SMS::edge_collapse(sm, stop, - CGAL::parameters::filter(SMS::Bounded_normal_change_filter<>()) - .get_cost(SMS::LindstromTurk_cost()) - .get_placement(SMS::LindstromTurk_placement())); + SMS::Face_count_stop_predicate stop(1); + SMS::edge_collapse( + surface_mesh, + stop, + CGAL::parameters:: + filter(SMS::Bounded_normal_change_filter>(SMS::Bounding_box_filter<>())) + .get_cost(SMS::LindstromTurk_cost()) + .get_placement(SMS::LindstromTurk_placement()) + ); CGAL::IO::write_OFF(std::cout, sm, CGAL::parameters::stream_precision(17)); for(auto v : vertices(sm)) From 95b4eba11e8f83a6e7748ee34456dadf08b0bb26 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o=20Valque?= Date: Thu, 6 Feb 2025 09:35:02 +0100 Subject: [PATCH 08/19] Fix of issue 8213 by consider length of a vector being zero if enough small, define diff_of_product in a specific file --- .../include/CGAL/Cartesian/CrossProduct.h | 60 +++++++++++++++++++ .../include/CGAL/Cartesian/MatrixC33.h | 14 +++++ .../internal/Lindstrom_Turk_core.h | 53 +++++----------- .../issue_8213.cpp | 16 ++--- 4 files changed, 99 insertions(+), 44 deletions(-) create mode 100644 Surface_mesh_simplification/include/CGAL/Cartesian/CrossProduct.h diff --git a/Surface_mesh_simplification/include/CGAL/Cartesian/CrossProduct.h b/Surface_mesh_simplification/include/CGAL/Cartesian/CrossProduct.h new file mode 100644 index 00000000000..ebaed862213 --- /dev/null +++ b/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 + +#include + +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 + 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 // + + diff --git a/Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h b/Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h index 3c24abbc3c2..70399434bb2 100644 --- a/Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h +++ b/Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h @@ -17,6 +17,7 @@ #include #include #include +#include #include @@ -219,6 +220,19 @@ std::optional< MatrixC33 > inverse_matrix(const MatrixC33& 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; diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h index 8c04b537e7f..247fa23381e 100644 --- a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h +++ b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h @@ -18,6 +18,7 @@ #include #include +#include #include #include @@ -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 - 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 @@ -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 diff --git a/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp b/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp index b086fe22e73..d8bd43e4a74 100644 --- a/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp +++ b/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp @@ -3,6 +3,8 @@ #include +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 #include -#include +// #include #include #include #include @@ -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 stop(1); SMS::edge_collapse( - surface_mesh, + sm, stop, CGAL::parameters:: - filter(SMS::Bounded_normal_change_filter>(SMS::Bounding_box_filter<>())) - .get_cost(SMS::LindstromTurk_cost()) + // filter(SMS::Bounded_normal_change_filter<>()) //>(SMS::Bounding_box_filter<>())) + get_cost(SMS::LindstromTurk_cost()) .get_placement(SMS::LindstromTurk_placement()) ); 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; From cc5297554c7e2bee65cb44acd1866be58d1938bd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o=20Valque?= Date: Tue, 1 Apr 2025 09:48:25 +0200 Subject: [PATCH 09/19] Correct code with suggestion of sebastien --- .../Edge_collapse/internal/Lindstrom_Turk_core.h | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h index 247fa23381e..d9ea2691883 100644 --- a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h +++ b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h @@ -320,17 +320,16 @@ 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())}); + //TODO for obscur reason, computing this maxBb increase running time by 10% + maxBb=(std::max)({maxBb,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; @@ -343,7 +342,6 @@ extract_triangle_data() mTriangle_data.push_back(Triangle_data(lNormalV,lNormalL)); } - maxBb= abs_max; } template @@ -662,7 +660,7 @@ add_constraint_if_alpha_compatible(const Vector& Ai, CGAL_SMS_LT_TRACE(3, " l: " << n_to_string(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 + // if bi 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)) From c63c2a24e2afa83277b4107489267b9d9ff8f966 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Thu, 22 May 2025 14:51:31 +0100 Subject: [PATCH 10/19] Activate Leo's code --- .../edge_collapse_bounded_normal_change.cpp | 5 ++-- .../include/CGAL/Cartesian/MatrixC33.h | 23 ++++++++++--------- .../internal/Lindstrom_Turk_core.h | 16 ++++++------- 3 files changed, 23 insertions(+), 21 deletions(-) diff --git a/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_bounded_normal_change.cpp b/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_bounded_normal_change.cpp index 04561d2239c..e4d0de6ed33 100644 --- a/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_bounded_normal_change.cpp +++ b/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_bounded_normal_change.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include @@ -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; diff --git a/Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h b/Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h index 70399434bb2..b2b3d59f51b 100644 --- a/Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h +++ b/Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h @@ -220,19 +220,19 @@ std::optional< MatrixC33 > inverse_matrix(const MatrixC33& 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 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 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; @@ -244,6 +244,7 @@ std::optional< MatrixC33 > inverse_matrix(const MatrixC33& 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, diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h index d9ea2691883..df70aabd1c1 100644 --- a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h +++ b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h @@ -319,7 +319,7 @@ LindstromTurkCore:: 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); @@ -327,9 +327,9 @@ extract_triangle_data() const Point_reference p2 = get_point(tri.v2); //TODO for obscur reason, computing this maxBb increase running time by 10% - maxBb=(std::max)({maxBb,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())}); + 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; @@ -342,6 +342,7 @@ extract_triangle_data() mTriangle_data.push_back(Triangle_data(lNormalV,lNormalL)); } + maxBb *= 2.0; // to avoid numerical problems } template @@ -659,11 +660,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)); - // 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) + // Due to double number type, l may have a small value instead of zero (example sum of the face normals of a tetrahedron for volumic constraint) // if bi 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)) + CGAL_SMS_LT_TRACE(3, " error consider: " << (CGAL::abs(bi) / maxBb)); + if(l > (std::abs(bi) / maxBb)) { Vector Ain = Ai / l; FT bin = bi / l; From 0f86fa6d23f0ace24f75437f84b18e62edfb43ce Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Thu, 22 May 2025 15:20:14 +0100 Subject: [PATCH 11/19] clean up in the cross product alternatives --- .../include/CGAL/Cartesian/CrossProduct.h | 1 + .../internal/Lindstrom_Turk_core.h | 22 ++++++++++--------- 2 files changed, 13 insertions(+), 10 deletions(-) diff --git a/Surface_mesh_simplification/include/CGAL/Cartesian/CrossProduct.h b/Surface_mesh_simplification/include/CGAL/Cartesian/CrossProduct.h index ebaed862213..bde522d5930 100644 --- a/Surface_mesh_simplification/include/CGAL/Cartesian/CrossProduct.h +++ b/Surface_mesh_simplification/include/CGAL/Cartesian/CrossProduct.h @@ -52,6 +52,7 @@ namespace CGAL { return a*b - c*d; } + } // namespace CGAL #endif // CGAL_CARTESIAN_CROSSPRODUCT_H // diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h index df70aabd1c1..d11e1eaf931 100644 --- a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h +++ b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h @@ -152,7 +152,8 @@ private : } #endif - static Vector SL_cross_product(const Vector& a, const Vector& b) + // balanced solution based on abusing the fact that here we expect u and v to have similar coordinates + static Vector robust_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(); @@ -184,9 +185,10 @@ private : exv.exact(); return to_approx(exv); } -#endif - static Vector X_product(const Vector& u, const Vector& v) + + + static Vector experimental_cross_product(const Vector& u, const Vector& v) { #if 0 // this can create large errors and spiky meshes for kernels with inexact constructions @@ -201,18 +203,18 @@ private : 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 } +#endif + + static Vector point_cross_product(const Point& a, const Point& b) { - return X_product(a-ORIGIN, b-ORIGIN); + return robust_cross_product(a-ORIGIN, b-ORIGIN); } // This is the (uX)(Xu) product described in the Lindstrom-Turk paper @@ -334,7 +336,7 @@ extract_triangle_data() Vector v01 = p1 - p0; Vector v02 = p2 - p0; - Vector lNormalV = X_product(v01,v02); + Vector lNormalV = robust_cross_product(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 @@ -660,10 +662,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)); - // Due to double number type, l may have a small value instead of zero (example sum of the face normals of a tetrahedron for volumic constraint) + // 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 > (std::abs(bi) / maxBb)) + if(l > (CGAL::abs(bi) / maxBb)) { Vector Ain = Ai / l; FT bin = bi / l; From 6acfaaedac8774103d10775b8045ed0b4ab0a1a5 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Thu, 22 May 2025 15:37:41 +0100 Subject: [PATCH 12/19] fix minor --- .../Edge_collapse/internal/Lindstrom_Turk_core.h | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h index d11e1eaf931..e45d198b1ed 100644 --- a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h +++ b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h @@ -155,10 +155,14 @@ private : // balanced solution based on abusing the fact that here we expect u and v to have similar coordinates static Vector robust_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(); + const FT& ax=a.x(); + const FT& ay=a.y(); + const FT& az=a.z(); + const FT& bx=b.x(); + const FT& by=b.y(); + const FT& bz=b.z(); - auto minor = [](auto ai, auto bi, auto aj, auto bj) + auto minor = [](const FT& ai, const FT& bi, const FT& aj, const FT& 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 From 1d9c84f9ef54426e27c0ed2584a86a98b6dde11d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Thu, 3 Jul 2025 09:41:43 +0200 Subject: [PATCH 13/19] static -> inline --- .../include/CGAL/Cartesian/CrossProduct.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Surface_mesh_simplification/include/CGAL/Cartesian/CrossProduct.h b/Surface_mesh_simplification/include/CGAL/Cartesian/CrossProduct.h index bde522d5930..a0882d63b61 100644 --- a/Surface_mesh_simplification/include/CGAL/Cartesian/CrossProduct.h +++ b/Surface_mesh_simplification/include/CGAL/Cartesian/CrossProduct.h @@ -19,7 +19,7 @@ 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) + 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); @@ -27,7 +27,7 @@ namespace CGAL { return f + e; } - static double diff_of_products_cht(const double a, const double b, const double c, const double d) + 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; @@ -38,7 +38,7 @@ namespace CGAL { return r + e; } - static double diff_of_products(const double a, const double b, const double c, const double d) + inline 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 @@ -47,7 +47,7 @@ namespace CGAL { } template - static OFT diff_of_products(const OFT& a, const OFT& b, const OFT& c, const OFT& d) + inline OFT diff_of_products(const OFT& a, const OFT& b, const OFT& c, const OFT& d) { return a*b - c*d; } From 388632e0fa74520b24a49850b739c853266c802d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Thu, 3 Jul 2025 09:46:32 +0200 Subject: [PATCH 14/19] fix warning --- .../CGAL/Surface_mesh_simplification/internal/Edge_collapse.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/internal/Edge_collapse.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/internal/Edge_collapse.h index c2259c55504..2335bb5a843 100644 --- a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/internal/Edge_collapse.h +++ b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/internal/Edge_collapse.h @@ -610,7 +610,7 @@ loop() auto mPQ_clone = *mPQ; std::optional opt_th; - while(opt_th = mPQ_clone.extract_top()) + 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(); From 02883a41969f649bb28d633667d177313dda63c3 Mon Sep 17 00:00:00 2001 From: lvalque Date: Tue, 15 Jul 2025 18:35:10 +0200 Subject: [PATCH 15/19] rename CrossProduct.h to robust_cross_product.h and move it in the appropriate directory --- .../include/CGAL/Cartesian/CrossProduct.h | 61 -------- .../include/CGAL/Cartesian/MatrixC33.h | 20 +-- .../internal/Lindstrom_Turk_core.h | 70 +--------- .../CGAL/internal/robust_cross_product.h | 131 ++++++++++++++++++ 4 files changed, 144 insertions(+), 138 deletions(-) delete mode 100644 Surface_mesh_simplification/include/CGAL/Cartesian/CrossProduct.h create mode 100644 Surface_mesh_simplification/include/CGAL/internal/robust_cross_product.h diff --git a/Surface_mesh_simplification/include/CGAL/Cartesian/CrossProduct.h b/Surface_mesh_simplification/include/CGAL/Cartesian/CrossProduct.h deleted file mode 100644 index a0882d63b61..00000000000 --- a/Surface_mesh_simplification/include/CGAL/Cartesian/CrossProduct.h +++ /dev/null @@ -1,61 +0,0 @@ -// 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 - -#include - -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 - 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) - { - // 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 - inline 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 // - - diff --git a/Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h b/Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h index b2b3d59f51b..096ad525c9b 100644 --- a/Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h +++ b/Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h @@ -17,7 +17,7 @@ #include #include #include -#include +#include #include @@ -221,17 +221,17 @@ std::optional< MatrixC33 > inverse_matrix(const MatrixC33& 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 c00 = ::CGAL::Surface_mesh_simplification::internal::diff_of_products(m.r1().y(),m.r2().z(),m.r2().y(),m.r1().z()) / det; + RT c01 = ::CGAL::Surface_mesh_simplification::internal::diff_of_products(m.r2().y(),m.r0().z(),m.r0().y(),m.r2().z()) / det; + RT c02 = ::CGAL::Surface_mesh_simplification::internal::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 c10 = ::CGAL::Surface_mesh_simplification::internal::diff_of_products(m.r2().x(),m.r1().z(),m.r1().x(),m.r2().z()) / det; + RT c11 = ::CGAL::Surface_mesh_simplification::internal::diff_of_products(m.r0().x(),m.r2().z(),m.r2().x(),m.r0().z()) / det; + RT c12 = ::CGAL::Surface_mesh_simplification::internal::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 c20 = ::CGAL::Surface_mesh_simplification::internal::diff_of_products(m.r1().x(),m.r2().y(),m.r2().x(),m.r1().y()) / det; + RT c21 = ::CGAL::Surface_mesh_simplification::internal::diff_of_products(m.r2().x(),m.r0().y(),m.r0().x(),m.r2().y()) / det; + RT c22 = ::CGAL::Surface_mesh_simplification::internal::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; diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h index e45d198b1ed..83589a60d51 100644 --- a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h +++ b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/Lindstrom_Turk_core.h @@ -18,7 +18,7 @@ #include #include -#include +#include #include #include @@ -152,73 +152,10 @@ private : } #endif - // balanced solution based on abusing the fact that here we expect u and v to have similar coordinates - static Vector robust_cross_product(const Vector& a, const Vector& b) - { - const FT& ax=a.x(); - const FT& ay=a.y(); - const FT& az=a.z(); - const FT& bx=b.x(); - const FT& by=b.y(); - const FT& bz=b.z(); - - auto minor = [](const FT& ai, const FT& bi, const FT& aj, const FT& 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 to_exact; - CGAL::Cartesian_converter to_approx; - auto exv = cross_product(to_exact(a), to_exact(b)); - exv.exact(); - return to_approx(exv); - } - - - - static Vector experimental_cross_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 0 - // obviously too slow - return exact_cross_product(u, v); -#endif - } - -#endif - static Vector point_cross_product(const Point& a, const Point& b) { - return robust_cross_product(a-ORIGIN, b-ORIGIN); + return robust_cross_product(a-ORIGIN, b-ORIGIN); } // This is the (uX)(Xu) product described in the Lindstrom-Turk paper @@ -332,7 +269,6 @@ extract_triangle_data() const Point_reference p1 = get_point(tri.v1); const Point_reference p2 = get_point(tri.v2); - //TODO for obscur reason, computing this maxBb increase running time by 10% 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())}); @@ -340,7 +276,7 @@ extract_triangle_data() Vector v01 = p1 - p0; Vector v02 = p2 - p0; - Vector lNormalV = robust_cross_product(v01,v02); + Vector lNormalV = robust_cross_product(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 diff --git a/Surface_mesh_simplification/include/CGAL/internal/robust_cross_product.h b/Surface_mesh_simplification/include/CGAL/internal/robust_cross_product.h new file mode 100644 index 00000000000..b4f9c0e7b37 --- /dev/null +++ b/Surface_mesh_simplification/include/CGAL/internal/robust_cross_product.h @@ -0,0 +1,131 @@ +// 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 + +#include +#include + +namespace CGAL::Surface_mesh_simplification::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) + { + // 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 + 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 + typename Geom_traits::Vector_3 robust_cross_product(const typename Geom_traits::Vector_3& a, const typename Geom_traits::Vector_3& b) + { + using FT = typename Geom_traits::FT; + using Vector = typename Geom_traits::Vector_3; + + const FT& ax=a.x(); + const FT& ay=a.y(); + const FT& az=a.z(); + const FT& bx=b.x(); + const FT& by=b.y(); + const FT& bz=b.z(); + + auto minor = [](const FT& ai, const FT& bi, const FT& aj, const FT& 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 + Vector exact_cross_product(const Vector& a, const Vector& b) + { + CGAL::Cartesian_converter to_exact; + CGAL::Cartesian_converter to_approx; + auto exv = cross_product(to_exact(a), to_exact(b)); + exv.exact(); + return to_approx(exv); + } + + + + Vector experimental_cross_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 0 + // obviously too slow + return exact_cross_product(u, v); +#elif 1 + // balanced solution based on abusing the fact that here we expect u and v to have similar coordinates + return robust_cross_product(u, v); +#endif + } + +#endif + + +} // namespace CGAL + +#endif // CGAL_SMS_INTERNAL_ROBUST_CROSS_PRODUCT_H + + From f51867ac3674e0b7970d31cf1b5d181da2cef0d4 Mon Sep 17 00:00:00 2001 From: lvalque Date: Wed, 16 Jul 2025 12:09:09 +0200 Subject: [PATCH 16/19] delete Boundinx_box_filter.h --- .../Edge_collapse/Bounding_box_filter.h | 77 ------------------- 1 file changed, 77 deletions(-) delete mode 100644 Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounding_box_filter.h diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounding_box_filter.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounding_box_filter.h deleted file mode 100644 index c8b534f0ee1..00000000000 --- a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounding_box_filter.h +++ /dev/null @@ -1,77 +0,0 @@ -// 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 -#include -#include - -#include - -namespace CGAL { -namespace Surface_mesh_simplification { - - -template -class Bounding_box_filter -{ -public: - Bounding_box_filter(const BaseFilter& base_filter = BaseFilter()) - : m_base_filter(base_filter) - {} - - - template - std::optional - operator()(const Profile& profile, std::optional 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::value_type Point; - typedef typename boost::property_traits::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(); - } - } - 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 From 40180827842f8802da9ed3a940d94658f9d58da6 Mon Sep 17 00:00:00 2001 From: lvalque Date: Wed, 16 Jul 2025 12:10:08 +0200 Subject: [PATCH 17/19] Variants as if maccro instead of comment lines --- .../include/CGAL/internal/robust_cross_product.h | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/Surface_mesh_simplification/include/CGAL/internal/robust_cross_product.h b/Surface_mesh_simplification/include/CGAL/internal/robust_cross_product.h index b4f9c0e7b37..9dbb4f460ef 100644 --- a/Surface_mesh_simplification/include/CGAL/internal/robust_cross_product.h +++ b/Surface_mesh_simplification/include/CGAL/internal/robust_cross_product.h @@ -41,10 +41,15 @@ namespace CGAL::Surface_mesh_simplification::internal{ inline double diff_of_products(const double a, const double b, const double c, const double d) { - // return a*b - c*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); - // return diff_of_products_cht(a, b, c, d); + #elif 0 + return diff_of_products_cht(a, b, c, d); + #endif } template @@ -86,7 +91,8 @@ namespace CGAL::Surface_mesh_simplification::internal{ } #if 0 - Vector exact_cross_product(const Vector& a, const Vector& b) + template + typename Geom_traits::Vector_3 exact_cross_product(const typename Geom_traits::Vector_3& a, const typename Geom_traits::Vector_3& b) { CGAL::Cartesian_converter to_exact; CGAL::Cartesian_converter to_approx; @@ -97,7 +103,8 @@ namespace CGAL::Surface_mesh_simplification::internal{ - Vector experimental_cross_product(const Vector& u, const Vector& v) + template + typename Geom_traits::Vector_3 experimental_cross_product(const typename Geom_traits::Vector_3& a, const typename Geom_traits::Vector_3& b) { #if 0 // this can create large errors and spiky meshes for kernels with inexact constructions From 13f5968eaac711649a5151662860dc1f4661a11d Mon Sep 17 00:00:00 2001 From: Mael Date: Thu, 17 Jul 2025 11:37:58 +0200 Subject: [PATCH 18/19] Clean up --- .../CGAL/internal/robust_cross_product.h | 2 +- .../data/issue_8213.off | 21 ++++++++----------- .../issue_8213.cpp | 16 ++++---------- 3 files changed, 14 insertions(+), 25 deletions(-) diff --git a/Surface_mesh_simplification/include/CGAL/internal/robust_cross_product.h b/Surface_mesh_simplification/include/CGAL/internal/robust_cross_product.h index 9dbb4f460ef..ea3b05fbcd2 100644 --- a/Surface_mesh_simplification/include/CGAL/internal/robust_cross_product.h +++ b/Surface_mesh_simplification/include/CGAL/internal/robust_cross_product.h @@ -18,7 +18,7 @@ namespace CGAL::Surface_mesh_simplification::internal{ -// a*b - c*d + // 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) { diff --git a/Surface_mesh_simplification/test/Surface_mesh_simplification/data/issue_8213.off b/Surface_mesh_simplification/test/Surface_mesh_simplification/data/issue_8213.off index 4669828bf08..14e2cf1fbe3 100644 --- a/Surface_mesh_simplification/test/Surface_mesh_simplification/data/issue_8213.off +++ b/Surface_mesh_simplification/test/Surface_mesh_simplification/data/issue_8213.off @@ -1,16 +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 - - +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 diff --git a/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp b/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp index d8bd43e4a74..24f89279cf8 100644 --- a/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp +++ b/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp @@ -3,8 +3,6 @@ #include -CGAL::Bbox_3 bbox_g; -double max_bbox_g; #define CGAL_CHECK_EXPENSIVE #define CGAL_SURFACE_SIMPLIFICATION_ENABLE_TRACE 5 @@ -17,7 +15,6 @@ void Surface_simplification_external_trace(const std::string& s) #include #include -// #include #include #include #include @@ -50,20 +47,15 @@ int main(int argc, char** argv) return EXIT_FAILURE; } - auto bb=PMP::bbox(sm); + 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 stop(1); - SMS::edge_collapse( - sm, - stop, - CGAL::parameters:: - // filter(SMS::Bounded_normal_change_filter<>()) //>(SMS::Bounding_box_filter<>())) - get_cost(SMS::LindstromTurk_cost()) - .get_placement(SMS::LindstromTurk_placement()) + SMS::edge_collapse(sm, stop, + CGAL::parameters::get_cost(SMS::LindstromTurk_cost()) + .get_placement(SMS::LindstromTurk_placement()) ); CGAL::IO::write_OFF(std::cout, sm, CGAL::parameters::stream_precision(17)); From 263dc8b49bd94aa2cb260df2c400975e6f3666f4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 17 Jul 2025 12:16:18 +0200 Subject: [PATCH 19/19] Misc fixes --- .../include/CGAL/Cartesian/MatrixC33.h | 59 ++--- .../CGAL/internal/robust_cross_product.h | 215 +++++++++--------- .../issue_8213.cpp | 7 +- 3 files changed, 139 insertions(+), 142 deletions(-) diff --git a/Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h b/Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h index 096ad525c9b..a41c8e768ca 100644 --- a/Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h +++ b/Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h @@ -13,11 +13,12 @@ #include +#include + #include #include #include #include -#include #include @@ -184,17 +185,19 @@ MatrixC33 cofactors_matrix(const MatrixC33& m) { typedef typename R::RT RT; - RT c00 = determinant(m.r1().y(),m.r1().z(),m.r2().y(),m.r2().z()); - 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()); + using ::CGAL::Surface_mesh_simplification::internal::diff_of_products; - RT c10 = -determinant(m.r0().y(),m.r0().z(),m.r2().y(),m.r2().z()); - RT c11 = determinant(m.r0().x(),m.r0().z(),m.r2().x(),m.r2().z()); - RT c12 = -determinant(m.r0().x(),m.r0().y(),m.r2().x(),m.r2().y()); + RT c00 = diff_of_products(m.r1().y(), m.r2().z(), m.r2().y(), m.r1().z()); + RT c01 = -diff_of_products(m.r1().x(), m.r2().z(), m.r2().x(), m.r1().z()); + 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 c21 = -determinant(m.r0().x(),m.r0().z(),m.r1().x(),m.r1().z()); - RT c22 = determinant(m.r0().x(),m.r0().y(),m.r1().x(),m.r1().y()); + RT c10 = -diff_of_products(m.r0().y(), m.r2().z(), m.r2().y(), m.r0().z()); + RT c11 = diff_of_products(m.r0().x(), m.r2().z(), m.r2().x(), m.r0().z()); + 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(c00,c01,c02, c10,c11,c12, @@ -212,7 +215,9 @@ std::optional< MatrixC33 > inverse_matrix(const MatrixC33& m) { typedef typename R::RT RT; typedef MatrixC33 Matrix; - typedef std::optional result_type; + typedef std::optional result_type; + + using ::CGAL::Surface_mesh_simplification::internal::diff_of_products; result_type rInverse; @@ -220,31 +225,17 @@ std::optional< MatrixC33 > inverse_matrix(const MatrixC33& m) if(! CGAL_NTS is_zero(det)) { -#if 1 - RT c00 = ::CGAL::Surface_mesh_simplification::internal::diff_of_products(m.r1().y(),m.r2().z(),m.r2().y(),m.r1().z()) / det; - RT c01 = ::CGAL::Surface_mesh_simplification::internal::diff_of_products(m.r2().y(),m.r0().z(),m.r0().y(),m.r2().z()) / det; - RT c02 = ::CGAL::Surface_mesh_simplification::internal::diff_of_products(m.r0().y(),m.r1().z(),m.r1().y(),m.r0().z()) / 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 = ::CGAL::Surface_mesh_simplification::internal::diff_of_products(m.r2().x(),m.r1().z(),m.r1().x(),m.r2().z()) / det; - RT c11 = ::CGAL::Surface_mesh_simplification::internal::diff_of_products(m.r0().x(),m.r2().z(),m.r2().x(),m.r0().z()) / det; - RT c12 = ::CGAL::Surface_mesh_simplification::internal::diff_of_products(m.r1().x(),m.r0().z(),m.r0().x(),m.r1().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 = ::CGAL::Surface_mesh_simplification::internal::diff_of_products(m.r1().x(),m.r2().y(),m.r2().x(),m.r1().y()) / det; - RT c21 = ::CGAL::Surface_mesh_simplification::internal::diff_of_products(m.r2().x(),m.r0().y(),m.r0().x(),m.r2().y()) / det; - RT c22 = ::CGAL::Surface_mesh_simplification::internal::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; - - RT c10 = (m.r1().z()*m.r2().x() - m.r1().x()*m.r2().z()) / det; - RT c11 = (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 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 + 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; rInverse = result_type(Matrix(c00,c01,c02, c10,c11,c12, diff --git a/Surface_mesh_simplification/include/CGAL/internal/robust_cross_product.h b/Surface_mesh_simplification/include/CGAL/internal/robust_cross_product.h index ea3b05fbcd2..e8e2caf095f 100644 --- a/Surface_mesh_simplification/include/CGAL/internal/robust_cross_product.h +++ b/Surface_mesh_simplification/include/CGAL/internal/robust_cross_product.h @@ -14,124 +14,129 @@ #include #include + #include -namespace CGAL::Surface_mesh_simplification::internal{ +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; - } +// 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 - 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 - typename Geom_traits::Vector_3 robust_cross_product(const typename Geom_traits::Vector_3& a, const typename Geom_traits::Vector_3& b) - { - using FT = typename Geom_traits::FT; - using Vector = typename Geom_traits::Vector_3; - - const FT& ax=a.x(); - const FT& ay=a.y(); - const FT& az=a.z(); - const FT& bx=b.x(); - const FT& by=b.y(); - const FT& bz=b.z(); - - auto minor = [](const FT& ai, const FT& bi, const FT& aj, const FT& 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); - } +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 - template - typename Geom_traits::Vector_3 exact_cross_product(const typename Geom_traits::Vector_3& a, const typename Geom_traits::Vector_3& b) - { - CGAL::Cartesian_converter to_exact; - CGAL::Cartesian_converter to_approx; - auto exv = cross_product(to_exact(a), to_exact(b)); - exv.exact(); - return to_approx(exv); - } - - - - template - typename Geom_traits::Vector_3 experimental_cross_product(const typename Geom_traits::Vector_3& a, const typename Geom_traits::Vector_3& b) - { -#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); + // this can create large errors with inexact constructions + return a*b - c*d; + // the next two are equivalent in results and speed #elif 1 - // balanced solution based on abusing the fact that here we expect u and v to have similar coordinates - return robust_cross_product(u, v); + return diff_of_products_kahan(a, b, c, d); +#elif 0 + return diff_of_products_cht(a, b, c, d); #endif - } +} +template +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 +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 +typename GeomTraits::Vector_3 exact_cross_product(const typename GeomTraits::Vector_3& a, + const typename GeomTraits::Vector_3& b) +{ + CGAL::Cartesian_converter to_exact; + CGAL::Cartesian_converter to_approx; + auto exv = cross_product(to_exact(a), to_exact(b)); + exv.exact(); + return to_approx(exv); +} #endif +template +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(u, v); +#endif +} } // namespace CGAL +} // namespace Surface_mesh_simplification +} // namespace internal #endif // CGAL_SMS_INTERNAL_ROBUST_CROSS_PRODUCT_H diff --git a/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp b/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp index 24f89279cf8..2e9811ca6ee 100644 --- a/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp +++ b/Surface_mesh_simplification/test/Surface_mesh_simplification/issue_8213.cpp @@ -30,6 +30,7 @@ namespace PMP = CGAL::Polygon_mesh_processing; using Kernel = CGAL::Simple_cartesian; // using Kernel = CGAL::Exact_predicates_exact_constructions_kernel; + using Surface_mesh = CGAL::Surface_mesh; int main(int argc, char** argv) @@ -53,10 +54,10 @@ int main(int argc, char** argv) std::cout << "Input mesh has " << num_faces(sm) << " faces" << std::endl; SMS::Face_count_stop_predicate stop(1); - SMS::edge_collapse(sm, stop, + SMS::edge_collapse(sm, stop, CGAL::parameters::get_cost(SMS::LindstromTurk_cost()) - .get_placement(SMS::LindstromTurk_placement()) - ); + .get_placement(SMS::LindstromTurk_placement())); + CGAL::IO::write_OFF(std::cout, sm, CGAL::parameters::stream_precision(17)); for(auto v : vertices(sm))