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;