Merge pull request #9048 from LeoValque/SMS-slow_Garland-Heckbert

New line quadrics for Garland-Heckbert policies

Documentation/doc/biblio/geom.bib

@@ -152118,3 +152118,15 @@ keywords = {polygonal surface mesh, Surface reconstruction, kinetic framework, s
   pages={11},
   year={2015}
 }
+
+@article{liu2025linequadrics,
+  journal = {Computer Graphics Forum},
+  title = {{Controlling Quadric Error Simplification with Line Quadrics}},
+  author = {Liu, Hsueh-Ti Derek and Rahimzadeh, Mehdi and Zordan, Victor},
+  year = {2025},
+  publisher = {The Eurographics Association and John Wiley & Sons Ltd.},
+  ISSN = {1467-8659},
+  DOI = {10.1111/cgf.70184},
+  url = {https://doi.org/10.1111/cgf.70184},
+  year={2025}
+}

Installation/CHANGES.md

@@ -12,6 +12,15 @@ Release date: July 2026
   - `import_from_triangulation_3()` → `triangulation_3_to_lcc()`
 - The old function names are still available but marked as deprecated for backward compatibility.
 
+
+### [Surface Mesh Simplification](https://doc.cgal.org/6.2/Manual/packages.html#PkgSurfaceMeshSimplification)
+
+-   Added the class `CGAL::Surface_mesh_simplification::GarlandHeckbert_plane_and_line_policies`, which provides improved output for `CGAL::Surface_mesh_simplification::edge_collapse()`.
+    That class works the same as previous `GarlandHeckbert_policies`.
+    Its constructor accepts a `Mesh` and optional named parameters to set the weight of the line policy relative to the plane policy, set the boundary cost multiplier or provide vertex normals.
+-  **Breaking change**: `CGAL::Surface_mesh_simplification::GarlandHeckbert_policies.h` is now an alias of `CGAL::Surface_mesh_simplification::GarlandHeckbert_plane_and_line_policies.h` and is no longer deprecated.
+
+
 ## [Release 6.1](https://github.com/CGAL/cgal/releases/tag/v6.1)
 
 Release date: Sept 2025

STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h

@@ -185,6 +185,8 @@ CGAL_add_named_parameter(get_cost_policy_t, get_cost_policy, get_cost)
 CGAL_add_named_parameter(get_placement_policy_t, get_placement_policy, get_placement)
 CGAL_add_named_parameter(filter_t, filter, filter)
 CGAL_add_named_parameter(use_relaxed_order_t, use_relaxed_order, use_relaxed_order)
+CGAL_add_named_parameter(line_policies_weight_t, line_policies_weight, line_policies_weight)
+CGAL_add_named_parameter(discontinuity_multiplier_t, discontinuity_multiplier, discontinuity_multiplier)
 
 //to be documented
 CGAL_add_named_parameter(face_normal_t, face_normal, face_normal_map)

Surface_mesh_simplification/benchmark/Surface_mesh_simplification/CMakeLists.txt

@@ -0,0 +1,17 @@
+# Created by the script cgal_create_cmake_script_with_options
+# This is the CMake script for compiling a set of CGAL applications.
+
+cmake_minimum_required(VERSION 3.12...3.31)
+project(Surface_mesh_simplification_Examples)
+
+# CGAL and its components
+find_package(CGAL REQUIRED)
+
+find_package(Eigen3 3.1.0 QUIET) #(3.1.0 or greater)
+include(CGAL_Eigen3_support)
+if(TARGET CGAL::Eigen3_support)
+  create_single_source_cgal_program("bench_garland_heckbert.cpp")
+  target_link_libraries(bench_garland_heckbert PRIVATE CGAL::Eigen3_support)
+else()
+  message(STATUS "NOTICE: Garland-Heckbert polices require the Eigen library, which has not been found; related examples will not be compiled.")
+endif()

Surface_mesh_simplification/benchmark/Surface_mesh_simplification/bench_garland_heckbert.cpp

@@ -0,0 +1,156 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_count_ratio_stop_predicate.h>
+#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounded_normal_change_filter.h>
+#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_policies.h>
+
+#include <CGAL/Polygon_mesh_processing/autorefinement.h>
+
+#include <CGAL/Polygon_mesh_processing/manifoldness.h>
+#include <CGAL/Polygon_mesh_processing/distance.h>
+
+#include <CGAL/Surface_mesh_simplification/edge_collapse.h>
+
+#include <CGAL/IO/polygon_mesh_io.h>
+
+#include <chrono>
+#include <fstream>
+#include <iostream>
+#include <filesystem>
+#include <vector>
+
+#define TAG CGAL::Parallel_if_available_tag
+
+typedef CGAL::Simple_cartesian<double>                          Kernel;
+typedef Kernel::FT                                              FT;
+typedef Kernel::Point_3                                         Point_3;
+typedef CGAL::Surface_mesh<Point_3>                             Surface_mesh;
+
+namespace SMS = CGAL::Surface_mesh_simplification;
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+typedef SMS::GarlandHeckbert_plane_policies<Surface_mesh, Kernel>                  Classic_plane;
+typedef SMS::GarlandHeckbert_probabilistic_plane_policies<Surface_mesh, Kernel>    Prob_plane;
+typedef SMS::GarlandHeckbert_triangle_policies<Surface_mesh, Kernel>               Classic_tri;
+typedef SMS::GarlandHeckbert_probabilistic_triangle_policies<Surface_mesh, Kernel> Prob_tri;
+typedef SMS::GarlandHeckbert_policies<Surface_mesh, Kernel>                        Classic_plane_and_line;
+
+// Old policies composed with line policies
+typedef SMS::internal::GarlandHeckbert_line_policies<Surface_mesh, Kernel> Line_quadric;
+typedef SMS::internal::GarlandHeckbert_composed_policies<Surface_mesh, Kernel, Prob_plane, Line_quadric> Proba_plane_and_line;
+typedef SMS::internal::GarlandHeckbert_composed_policies<Surface_mesh, Kernel, Classic_tri, Line_quadric> Classic_tri_plus_line;
+typedef SMS::internal::GarlandHeckbert_composed_policies<Surface_mesh, Kernel, Prob_tri, Line_quadric> Proba_tri_plus_line;
+
+double mean_aspect_ratio(Surface_mesh& mesh)
+{
+  double total_aspect_ratio = 0;
+  double total_area = 0;
+  for(auto f: faces(mesh)){
+    // double a = PMP::face_area(f, mesh);
+    double a = 1.;
+    total_aspect_ratio += a*CGAL::to_double(PMP::face_aspect_ratio(f, mesh));
+    total_area += a;
+  }
+  return total_aspect_ratio/total_area;
+}
+
+template <typename GHPolicies>
+double collapse_gh(Surface_mesh& mesh,
+                   const double ratio,
+                   GHPolicies gh_policies)
+{
+  std::chrono::steady_clock::time_point start_time = std::chrono::steady_clock::now();
+
+  SMS::Edge_count_ratio_stop_predicate<Surface_mesh> stop(ratio);
+
+  // Garland&Heckbert simplification policies
+
+  typedef typename GHPolicies::Get_cost                                        GH_cost;
+  typedef typename GHPolicies::Get_placement                                   GH_placement;
+  typedef SMS::Bounded_normal_change_filter<>                                  Filter;
+
+  const GH_cost& gh_cost = gh_policies.get_cost();
+  const GH_placement& gh_placement = gh_policies.get_placement();
+
+  SMS::edge_collapse(mesh, stop, CGAL::parameters::get_cost(gh_cost).filter(Filter()).get_placement(gh_placement));
+
+  std::chrono::steady_clock::time_point end_time = std::chrono::steady_clock::now();
+
+  return (std::chrono::duration_cast<std::chrono::milliseconds>(end_time - start_time).count()) / 1000.;
+}
+
+std::string getFileNameWithoutExtension(const std::string& filePath)
+{
+  std::filesystem::path path(filePath);
+  return path.stem().string();  // 'stem' gets the filename without extension
+}
+
+int main(int argc, char** argv)
+{
+  Surface_mesh mesh;
+  const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/cube-meshed.off");
+  PMP::duplicate_non_manifold_vertices(mesh);
+  PMP::autorefine(mesh);
+
+  if(!CGAL::IO::read_polygon_mesh(filename, mesh))
+  {
+    std::cerr << "Failed to read input mesh: " << filename << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  if(!CGAL::is_triangle_mesh(mesh))
+  {
+    std::cerr << "Input geometry is not triangulated." << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  const double ratio = (argc > 2) ? std::stod(argv[2]) : 0.2;
+  std::cout << ratio << std::endl;
+  std::cout << "\n\nPolicy Approx_Hausdorff Dist_points_to_input mean_aspect_ratio running_time" << std::endl;
+  auto collapse=[&](std::string policy){
+    using CGAL::Surface_mesh_simplification::make_GarlandHeckbert_plane_and_line_policies;
+    Surface_mesh tmesh(mesh);
+    double time;
+    if(policy == "Classic_plane")
+      time=collapse_gh(tmesh, ratio, Classic_plane(tmesh));
+    else if(policy == "Classic_triangle")
+      time=collapse_gh(tmesh, ratio, Classic_tri(tmesh));
+    else if(policy == "Prob_plane")
+      time=collapse_gh(tmesh, ratio, Prob_plane(tmesh));
+    else if(policy == "Prob_triangle")
+      time=collapse_gh(tmesh, ratio, Prob_tri(tmesh));
+    // else if(policy == "Plane_plus_line_0.1")
+    //   time=collapse_gh(tmesh, ratio, make_GarlandHeckbert_plane_and_line_policies(tmesh, CGAL::parameters::line_policies_weight(0.1)));
+    else if(policy == "Plane_plus_line_0.01")
+      time=collapse_gh(tmesh, ratio, make_GarlandHeckbert_plane_and_line_policies(tmesh, CGAL::parameters::line_policies_weight(0.01)));
+    else if(policy == "Plane_plus_line_0.001")
+      time=collapse_gh(tmesh, ratio, make_GarlandHeckbert_plane_and_line_policies(tmesh, CGAL::parameters::line_policies_weight(0.001)));
+    // else if(policy == "Classic_tri_line")
+    //   time=collapse_gh(tmesh, ratio, Classic_tri_plus_line(tmesh, 100));
+    // else if(policy == "Proba_plane_line")
+    //   time=collapse_gh(tmesh, ratio, Proba_plane_and_line(tmesh, Prob_plane(tmesh), Line_quadric(tmesh), 100));
+    // else if(policy == "Proba_tri_line")
+    //   time=collapse_gh(tmesh, ratio, Proba_tri_plus_line(tmesh, Prob_tri(tmesh), Line_quadric(tmesh), 100));
+
+    std::cout << policy << " " << PMP::approximate_symmetric_Hausdorff_distance<TAG>(tmesh, mesh)
+                        << " " << PMP::max_distance_to_triangle_mesh<TAG>(tmesh.points(),mesh)
+                        << " " << mean_aspect_ratio(tmesh)
+                        << " " << time << std::endl;
+    CGAL::IO::write_polygon_mesh("out_"+policy+"_"+std::to_string(ratio)+".off", tmesh, CGAL::parameters::stream_precision(17));
+  };
+
+  collapse("Classic_plane");
+  collapse("Classic_triangle");
+  collapse("Prob_plane");
+  collapse("Prob_triangle");
+  // collapse("Plane_plus_line_0.1");
+  collapse("Plane_plus_line_0.01");
+  collapse("Plane_plus_line_0.001");
+  // collapse("Classic_tri_line");
+  // collapse("Proba_plane_line");
+  // collapse("Proba_tri_line");
+  // collapse("Plane_plus_line_0.001");
+
+  return EXIT_SUCCESS;
+}

Surface_mesh_simplification/doc/Surface_mesh_simplification/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_plane_and_line_policies.h

@@ -0,0 +1,96 @@
+namespace CGAL {
+namespace Surface_mesh_simplification {
+
+/*!
+\ingroup PkgSurfaceMeshSimplificationRef
+
+The class `GarlandHeckbert_plane_and_line_policies` regroups the cost and placement policies
+based on the Garland-Heckbert "Plane and line" strategy of Liu and colleagues \cgalCite{liu2025linequadrics}.
+
+This policy enhances the original Garland-Heckbert quadric error metrics
+by adding to the cost the distance to the line passing through the input vertices and aligned with their normals.
+Compared to the "classic" plane strategy, this strategy improves the speed and the quality of the result
+(see Section \ref SurfaceMeshSimplificationGarlandHeckbertStrategy).
+
+\note Both the cost and the placement policies must be used together as they internally use
+and share information associating quadrics to vertices.
+Note however, that they may still be wrapped with behavior-modifying classes
+such as `Constrained_placement` or `Bounded_normal_change_placement`.
+
+\tparam TriangleMesh is the type of surface mesh being simplified, and must be a model
+                     of the `MutableFaceGraph` and `HalfedgeListGraph` concepts.
+\tparam GeomTraits must be a model of `Kernel`. If you have passed a traits class in the optional
+                   named parameters in the call to `edge_collapse()`, the types must be identical.
+
+These policies depend on the third party \ref thirdpartyEigen library.
+
+\sa `GarlandHeckbert_policies`
+\sa `GarlandHeckbert_plane_policies`
+\sa `GarlandHeckbert_probabilistic_plane_policies`
+\sa `GarlandHeckbert_triangle_policies`
+\sa `GarlandHeckbert_probabilistic_triangle_policies`
+*/
+template <typename TriangleMesh, typename GeomTraits>
+class GarlandHeckbert_plane_and_line_policies
+{
+public:
+  /// The type of the Garland-Heckbert cost functor, a model of the concept `GetCost`
+  typedef unspecified_type Get_cost;
+
+  /// The type of the Garland-Heckbert placement functor, a model of the concept `GetPlacement`
+  typedef unspecified_type Get_placement;
+
+  /// \name Creation
+  /// @{
+
+  /*!
+  initializes the Garland-Heckbert plane and line policies.
+
+  \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+
+  \param tmesh the triangle mesh
+  \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+
+  \cgalNamedParamsBegin
+    \cgalParamNBegin{vertex_normal_map}
+      \cgalParamDescription{a property map associating to each vertex of `tmesh` a normal direction for that vertex}
+      \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<PolygonMesh>::%vertex_descriptor`
+                     as key type and `Vector_3` as value type}
+      \cgalParamDefault{an internal map filled using `CGAL::Polygon_mesh_processing::compute_vertex_normals`}
+    \cgalParamNEnd
+    \cgalParamNBegin{discontinuity_multiplier}
+      \cgalParamDescription{a multiplier of the error value for boundary edges to preserve the boundaries}
+      \cgalParamType{double}
+      \cgalParamDefault{100}
+    \cgalParamNEnd
+    \cgalParamNBegin{line_policies_weight}
+      \cgalParamDescription{a value that defines the weight of the line policies compared to the plane policies}
+      \cgalParamType{double}
+      \cgalParamDefault{0.01}
+    \cgalParamNEnd
+  \cgalNamedParamsEnd
+  */
+  template<typename NamedParameters = parameters::Default_named_parameters>
+  GarlandHeckbert_plane_and_line_policies(TriangleMesh& tmesh, const NamedParameters &np = parameters::default_values());
+
+  /// @}
+
+  /// \name Accessors
+  /// @{
+
+  const Get_cost& get_cost() const;
+  const Get_placement& get_placement() const;
+
+  /// @}
+
+};
+
+/*!
+creates a Garland-Heckbert plane and line policies object.
+*/
+template<typename TriangleMesh, typename NamedParameters = parameters::Default_named_parameters>
+GarlandHeckbert_plane_and_line_policies make_GarlandHeckbert_plane_and_line_policies(TriangleMesh& tmesh,
+                                                                                     const NamedParameters& np = parameters::default_values());
+
+} // namespace Surface_mesh_simplification
+} // namespace CGAL

Surface_mesh_simplification/doc/Surface_mesh_simplification/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_plane_policies.h

@@ -25,6 +25,7 @@ These policies depend on the third party \ref thirdpartyEigen library.
 \sa `GarlandHeckbert_probabilistic_plane_policies`
 \sa `GarlandHeckbert_triangle_policies`
 \sa `GarlandHeckbert_probabilistic_triangle_policies`
+\sa `GarlandHeckbert_plane_and_line_policies`
 */
 template <typename TriangleMesh, typename GeomTraits>
 class GarlandHeckbert_plane_policies

Surface_mesh_simplification/doc/Surface_mesh_simplification/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_policies.h

@@ -4,55 +4,19 @@ namespace Surface_mesh_simplification {
 /*!
 \ingroup PkgSurfaceMeshSimplificationRef
 
-\deprecated This class is deprecated since \cgal 5.5 and the introduction of variations
-of Garland-Heckbert policies (Section \ref SurfaceMeshSimplificationGarlandHeckbertStrategy).
-The class `GarlandHeckbert_plane_policies` is the modern equivalent to this class.
+This class is an alias for the current state-of-the-art Garland-Heckbert policies
+(see Section \ref SurfaceMeshSimplificationGarlandHeckbertStrategy).
+It is currently an alias of `GarlandHeckbert_plane_and_line_policies`.
 
-The class `GarlandHeckbert_policies` regroups the cost and placement policies
-based on the Garland-Heckbert strategy (Section \ref SurfaceMeshSimplificationGarlandHeckbertStrategy),
-both oracles must indeed be used together as they internally use and share information
-associating quadrics to vertices.
-
-Note however, that they may still be wrapped with behavior modifying classes
-such as `Constrained_placement` or `Bounded_normal_change_placement`.
-
-\tparam TriangleMesh is the type of surface mesh being simplified, and must be a model
-                     of the `MutableFaceGraph` and `HalfedgeListGraph` concepts.
-\tparam GeomTraits must be a model of `Kernel`. If you have passed a traits class in the optional
-                   named parameters in the call to `edge_collapse()`, the types must be identical.
-
-These policies depend on the third party \ref thirdpartyEigen library.
+\sa `GarlandHeckbert_plane_policies`
+\sa `GarlandHeckbert_probabilistic_plane_policies`
+\sa `GarlandHeckbert_triangle_policies`
+\sa `GarlandHeckbert_probabilistic_triangle_policies`
+\sa `GarlandHeckbert_plane_and_line_policies`
 
 */
-template <typename TriangleMesh, typename GeomTraits>
-class GarlandHeckbert_policies
-{
-public:
-  /// The type of the Garland Heckbert cost functor, a model of the concept `GetCost`
-  typedef unspecified_type Get_cost;
-
-  /// The type of the Garland Heckbert placement functor, a model of the concept `GetPlacement`
-  typedef unspecified_type Get_placement;
-
-  /// \name Creation
-  /// @{
-
-  /*!
-  Initializes the Garland-Heckbert policies.
-  */
-  GarlandHeckbert_policies(TriangleMesh& tmesh);
-
-  /// @}
-
-  /// \name Accessors
-  /// @{
-
-  const Get_cost& get_cost() const;
-  const Get_placement& get_placement() const;
-
-  /// @}
-
-};
+template<typename TriangleMesh, typename GeomTraits>
+using GarlandHeckbert_policies = GarlandHeckbert_plane_and_line_policies<TriangleMesh, GeomTraits>;
 
 } // namespace Surface_mesh_simplification
 } // namespace CGAL

Surface_mesh_simplification/doc/Surface_mesh_simplification/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_probabilistic_plane_policies.h

@@ -26,6 +26,7 @@ These policies depend on the third party \ref thirdpartyEigen library.
 \sa `GarlandHeckbert_plane_policies`
 \sa `GarlandHeckbert_triangle_policies`
 \sa `GarlandHeckbert_probabilistic_triangle_policies`
+\sa `GarlandHeckbert_plane_and_line_policies`
 */
 template <typename TriangleMesh, typename GeomTraits>
 class GarlandHeckbert_probabilistic_plane_policies

Surface_mesh_simplification/doc/Surface_mesh_simplification/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_probabilistic_triangle_policies.h

@@ -26,6 +26,7 @@ These policies depend on the third party \ref thirdpartyEigen library.
 \sa `GarlandHeckbert_plane_policies`
 \sa `GarlandHeckbert_probabilistic_plane_policies`
 \sa `GarlandHeckbert_triangle_policies`
+\sa `GarlandHeckbert_plane_and_line_policies`
 */
 template <typename TriangleMesh, typename GeomTraits>
 class GarlandHeckbert_probabilistic_triangle_policies

Surface_mesh_simplification/doc/Surface_mesh_simplification/Doxyfile.in

@@ -2,4 +2,5 @@
 
 PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - Triangulated Surface Mesh Simplification"
 
-HTML_EXTRA_FILES           =  ${CGAL_PACKAGE_DOC_DIR}/fig/sappho.jpg
+HTML_EXTRA_FILES           =  ${CGAL_PACKAGE_DOC_DIR}/fig/sappho.png \
+                              ${CGAL_PACKAGE_DOC_DIR}/fig/octocat.png

Surface_mesh_simplification/doc/Surface_mesh_simplification/PackageDescription.txt

@@ -49,7 +49,8 @@
 - `CGAL::Surface_mesh_simplification::GarlandHeckbert_probabilistic_plane_policies<TriangleMesh, GeomTraits>`
 - `CGAL::Surface_mesh_simplification::GarlandHeckbert_triangle_policies<TriangleMesh, GeomTraits>`
 - `CGAL::Surface_mesh_simplification::GarlandHeckbert_probabilistic_triangle_policies<TriangleMesh, GeomTraits>`
-- `CGAL::Surface_mesh_simplification::GarlandHeckbert_policies<TriangleMesh, GeomTraits>` (deprecated)
+- `CGAL::Surface_mesh_simplification::GarlandHeckbert_plane_and_line_policies<TriangleMesh, GeomTraits>`
+- `CGAL::Surface_mesh_simplification::GarlandHeckbert_policies<TriangleMesh, GeomTraits>`
 
 \cgalCRPSection{Policy Enhancements}
 - `CGAL::Surface_mesh_simplification::Constrained_placement<Placement, TriangleMesh>`

Surface_mesh_simplification/doc/Surface_mesh_simplification/Surface_mesh_simplification.txt

@@ -188,21 +188,38 @@ face normals). This variance naturally deteriorates the tightness of the result,
 hand it enables creating more uniform triangulations and the approach is more tolerant to noise,
 while still maintaining feature sensitivity.
 
+Another extension of the Garland–Heckbert approach was introduced by Liu et al. \cgalCite{liu2025linequadrics},
+who enhance the classic plane-based quadric error metrics with line quadrics.
+More precisely, line quadrics, which encode the distance to a line, are defined at all input vertices using the supporting lines of the normals.
+These quadrics are then combined with the classic plane quadrics, but with a small weight (0.01 by default).
+This results in more uniform triangulations, while preserving the tightness of the approximation.
+
 \cgalFigureAnchor{SurfaceMeshSimplification_GH}
 <center>
-<img src="sappho.jpg" style="max-width:70%;"/>
+<img src="sappho.png" style="max-width:70%;"/>
 </center>
 \cgalFigureCaptionBegin{SurfaceMeshSimplification_GH}
 <a href="https://ten-thousand-models.appspot.com/detail.html?file_id=72286">Sappho's Head</a> model (leftmost, 34882 vertices).
 From left to right, simplified output (1745 vertices) and symmetric Hausdorff distance
-for the four Garland-Heckbert variations: plane (0.217912), probabilistic
-plane (0.256801), triangle (0.268872), and probabilistic triangle (0.490846).
+for the five Garland-Heckbert variations: plane (0.217912), probabilistic
+plane (0.256801), triangle (0.268872), probabilistic triangle (0.490846) and plane and line (0.173344).
+\cgalFigureCaptionEnd
+
+\cgalFigureAnchor{SurfaceMeshSimplification_GH_octocat}
+<center>
+<img src="octocat.png" style="max-width:70%;"/>
+</center>
+\cgalFigureCaptionBegin{SurfaceMeshSimplification_GH_octocat}
+<a href="https://ten-thousand-models.appspot.com/detail.html?file_id=32770">Octocat</a> model (top-left, 18944 vertices).
+From top-left to bottom-right, simplified output (1895 vertices) and symmetric Hausdorff distance
+for the five Garland-Heckbert variations: plane (0.19631), probabilistic
+plane (0.59446), triangle (0.377533), probabilistic triangle (0.891286) and plane and line (0.304001).
 \cgalFigureCaptionEnd
 
 \subsection Surface_mesh_simplificationCostStrategyPolicies Cost Strategy Policies
 
 The cost strategy used by the algorithm is selected by means of three policies:
-`GetPlacement`, `GetCost`, and `Filter`.
+`GetPlacement`, `GetCost`, and `PlacementFilter`.
 
 The `GetPlacement` policy is called to compute the new position
 for the remaining vertex after the halfedge-collapse. It returns
@@ -400,8 +417,8 @@ steps of the simplification algorithm.
 
 Each Garland-Heckbert simplification strategy is implemented with a single class that regroups
 both the cost and the placement policies, which must be used together as they share vertex quadric data.
-The classic strategy of using plane-based quadric error metrics is implemented with the class
-`Surface_mesh_simplification::GarlandHeckbert_plane_policies`.
+The state-of-the-art strategy of combining plane and line-based quadric error metrics is implemented with the class
+`Surface_mesh_simplification::GarlandHeckbert_policies`.
 Although both policies must be used together, it is still possible to wrap either policy
 using behavior modifiers such as `Surface_mesh_simplification::Bounded_normal_change_placement`.
 
@@ -420,5 +437,7 @@ The implementation of the Garland-Heckbert simplification policies is the result
 of Baskın Şenbaşlar (Google Summer of Code 2019), and Julian Komaromy (Google Summer of Code 2021)
 They both were mentored by Mael Rouxel-Labbé, who also contributed to the code and to the documentation.
 
+Léo Valque added the `Surface_mesh_simplification::GarlandHeckbert_plane_and_line_policies` functionality in \cgal 6.2.
+
 */
 } /* namespace CGAL */

Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_garland_heckbert.cpp

@@ -2,7 +2,7 @@
 #include <CGAL/Surface_mesh.h>
 
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_count_ratio_stop_predicate.h>
-#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounded_normal_change_placement.h>
+#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounded_normal_change_filter.h>
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_policies.h>
 #include <CGAL/Surface_mesh_simplification/edge_collapse.h>
 
@@ -24,6 +24,7 @@ typedef SMS::GarlandHeckbert_plane_policies<Surface_mesh, Kernel>
 typedef SMS::GarlandHeckbert_probabilistic_plane_policies<Surface_mesh, Kernel>    Prob_plane;
 typedef SMS::GarlandHeckbert_triangle_policies<Surface_mesh, Kernel>               Classic_tri;
 typedef SMS::GarlandHeckbert_probabilistic_triangle_policies<Surface_mesh, Kernel> Prob_tri;
+typedef SMS::GarlandHeckbert_plane_and_line_policies<Surface_mesh, Kernel>         Plane_and_line;
 
 template <typename GHPolicies>
 void collapse_gh(Surface_mesh& mesh,
@@ -37,16 +38,17 @@ void collapse_gh(Surface_mesh& mesh,
 
   typedef typename GHPolicies::Get_cost                                        GH_cost;
   typedef typename GHPolicies::Get_placement                                   GH_placement;
-  typedef SMS::Bounded_normal_change_placement<GH_placement>                   Bounded_GH_placement;
+  typedef SMS::Bounded_normal_change_filter<>                                  Filter;
 
   GHPolicies gh_policies(mesh);
   const GH_cost& gh_cost = gh_policies.get_cost();
   const GH_placement& gh_placement = gh_policies.get_placement();
-  Bounded_GH_placement placement(gh_placement);
+  Filter filter;
 
   int r = SMS::edge_collapse(mesh, stop,
                              CGAL::parameters::get_cost(gh_cost)
-                                              .get_placement(placement));
+                                              .filter(filter)
+                                              .get_placement(gh_placement));
 
   std::chrono::steady_clock::time_point end_time = std::chrono::steady_clock::now();
 
@@ -59,7 +61,7 @@ void collapse_gh(Surface_mesh& mesh,
 
 // Usage:
 // ./command [input] [ratio] [policy] [output]
-// policy can be "cp" (classic plane), "ct" (classic triangle), "pp" (probabilistic plane), "pt" (probabilistic triangle)
+// policy can be "cp" (classic plane), "ct" (classic triangle), "pp" (probabilistic plane), "pt" (probabilistic triangle), "pl" (plane and line)
 int main(int argc, char** argv)
 {
   Surface_mesh mesh;
@@ -91,8 +93,10 @@ int main(int argc, char** argv)
     collapse_gh<Classic_tri>(mesh, ratio);
   else if(policy == "pp")
     collapse_gh<Prob_plane>(mesh, ratio);
-  else
+  else if(policy == "pt")
     collapse_gh<Prob_tri>(mesh, ratio);
+  else
+    collapse_gh<Plane_and_line>(mesh, ratio);
 
   CGAL::IO::write_polygon_mesh((argc > 4) ? argv[4] : "out.off", mesh, CGAL::parameters::stream_precision(17));
 

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

@@ -0,0 +1,126 @@
+// 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)     : Leo Valque
+
+#ifndef CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_GARLANDHECKBERT_PLANE_AND_LINE_POLICIES_H
+#define CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_GARLANDHECKBERT_PLANE_AND_LINE_POLICIES_H
+
+#include <CGAL/license/Surface_mesh_simplification.h>
+
+#include <CGAL/Surface_mesh_simplification/internal/Common.h>
+#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_policy_base.h>
+#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_functions.h>
+#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_composed_policies.h>
+#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_line_policies.h>
+
+#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_plane_policies.h>
+
+#include <CGAL/Named_function_parameters.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+#include <boost/property_map/function_property_map.hpp>
+
+namespace CGAL {
+namespace Surface_mesh_simplification {
+
+namespace internal {
+
+// Helpers to perform computation in the constructor GH_plane_and_line_policies before calling the constructor of Base
+template< typename TriangleMesh = void, typename NamedParameters = void>
+struct GH_helper{
+  typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type GT;
+  typedef typename GT::Vector_3 Vector_3;
+  typedef typename GT::FT FT;
+
+  typedef dynamic_vertex_property_t<Vector_3>                                       Vertex_normal_tag;
+  typedef typename boost::property_map<TriangleMesh, Vertex_normal_tag>::type       Vertex_normal_dmap;
+  typedef typename internal_np::Lookup_named_param_def<internal_np::vertex_normal_map_t,
+                                                       NamedParameters,
+                                                       Vertex_normal_dmap>::type    Vertex_normal_map;
+
+  NamedParameters np;
+  GH_helper(const NamedParameters &np_):np(np_){ }
+
+  Vertex_normal_map vnm(TriangleMesh& tmesh) const{
+    using parameters::choose_parameter;
+    using parameters::is_default_parameter;
+    using parameters::get_parameter;
+    constexpr bool must_compute_vertex_normals = is_default_parameter<NamedParameters, internal_np::vertex_normal_map_t>::value;
+
+    Vertex_normal_map vertex_normals = choose_parameter(get_parameter(np, internal_np::vertex_normal_map),
+                                                        get(Vertex_normal_tag(), tmesh));
+    if constexpr(must_compute_vertex_normals){
+      Polygon_mesh_processing::compute_vertex_normals(tmesh, vertex_normals, np);
+    }
+    return vertex_normals;
+  }
+
+  double lw() const{
+    using parameters::choose_parameter;
+    using parameters::get_parameter;
+    return choose_parameter(get_parameter(np, internal_np::line_policies_weight), 0.01);
+  }
+
+  double dm() const{
+    using parameters::choose_parameter;
+    using parameters::get_parameter;
+
+    return choose_parameter(get_parameter(np, internal_np::discontinuity_multiplier), 100);
+  }
+};
+
+} // namespace internal
+
+template<typename TriangleMesh, typename GeomTraits>
+class GarlandHeckbert_plane_and_line_policies
+  : public internal::GarlandHeckbert_composed_policies<TriangleMesh, GeomTraits,
+                                             GarlandHeckbert_plane_policies<TriangleMesh, GeomTraits>,
+                                             internal::GarlandHeckbert_line_policies<TriangleMesh, GeomTraits>,
+                                             true>
+{
+  typedef GarlandHeckbert_plane_policies<TriangleMesh, GeomTraits>                            GH_plane_polices;
+  typedef internal::GarlandHeckbert_line_policies<TriangleMesh, GeomTraits>                   GH_line_polices;
+  typedef internal::GarlandHeckbert_composed_policies<TriangleMesh, GeomTraits,
+                                                      GH_plane_polices,
+                                                      GH_line_polices,
+                                                      true>                                   Base;
+
+public:
+  typedef typename Base::Quadric_calculator Quadric_calculator;
+
+  typedef typename GeomTraits::FT                                              FT;
+
+public:
+  template<typename NP = parameters::Default_named_parameters>
+  GarlandHeckbert_plane_and_line_policies(TriangleMesh& tmesh, const NP& np = parameters::default_values()):
+    Base(tmesh,
+      GH_plane_polices(tmesh),
+      GH_line_polices(tmesh, internal::GH_helper<TriangleMesh,NP>(np).vnm(tmesh)),
+      FT(1.)/internal::GH_helper<TriangleMesh,NP>(np).lw(),
+      FT(1.),
+      internal::GH_helper<TriangleMesh,NP>(np).dm())
+  { }
+
+public:
+  using Base::operator();
+  using Base::get_cost;
+  using Base::get_placement;
+};
+
+template<typename TriangleMesh, typename NamedParameters = parameters::Default_named_parameters>
+GarlandHeckbert_plane_and_line_policies<TriangleMesh, typename GetGeomTraits<TriangleMesh,NamedParameters>::type>
+  make_GarlandHeckbert_plane_and_line_policies(TriangleMesh& tmesh,
+                                               const NamedParameters& np = parameters::default_values()){
+  typedef typename GetGeomTraits<TriangleMesh,NamedParameters>::type GT;
+  return GarlandHeckbert_plane_and_line_policies<TriangleMesh, GT>(tmesh, np);
+}
+
+} // namespace Surface_mesh_simplification
+} // namespace CGAL
+
+#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_GARLANDHECKBERT_PLANE_AND_LINE_POLICIES_H

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

@@ -32,6 +32,15 @@ class Plane_quadric_calculator
 public:
   Plane_quadric_calculator() { }
 
+  template <typename VertexPointMap>
+  Mat_4 construct_quadric_from_vertex(typename boost::graph_traits<TriangleMesh>::vertex_descriptor /*v*/,
+                                      const TriangleMesh& /*tmesh*/,
+                                      const VertexPointMap /*point_map*/,
+                                      const GeomTraits& /*gt*/) const
+  {
+    return Mat_4::Zero();
+  }
+
   template <typename VertexPointMap>
   Mat_4 construct_quadric_from_edge(typename boost::graph_traits<TriangleMesh>::halfedge_descriptor he,
                                     const TriangleMesh& tmesh,

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

@@ -21,12 +21,13 @@
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_triangle_policies.h>
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_probabilistic_plane_policies.h>
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_probabilistic_triangle_policies.h>
+#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_plane_and_line_policies.h>
 
 namespace CGAL {
 namespace Surface_mesh_simplification {
 
 template<typename TriangleMesh, typename GeomTraits>
-using GarlandHeckbert_policies CGAL_DEPRECATED = GarlandHeckbert_plane_policies<TriangleMesh, GeomTraits>;
+using GarlandHeckbert_policies = GarlandHeckbert_plane_and_line_policies<TriangleMesh, GeomTraits>;
 
 } // namespace Surface_mesh_simplification
 } // namespace CGAL

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

@@ -74,6 +74,15 @@ public:
   }
 
 public:
+  template <typename VertexPointMap>
+  Mat_4 construct_quadric_from_vertex(typename boost::graph_traits<TriangleMesh>::vertex_descriptor /*v*/,
+                                      const TriangleMesh& /*tmesh*/,
+                                      const VertexPointMap /*point_map*/,
+                                      const GeomTraits& /*gt*/) const
+  {
+    return Mat_4::Zero();
+  }
+
   template <typename VertexPointMap>
   Mat_4 construct_quadric_from_edge(const halfedge_descriptor he,
                                     const TriangleMesh& tmesh,

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

@@ -48,7 +48,7 @@ private:
 
   // this is only used when we fall back to probabilistic planes for the discontinuous edges,
   // the actual triangle quadric calculation only uses the normal variance
-  static constexpr FT position_variance_factor = 1;
+  static constexpr FT position_variance_factor = 0.1;
 
   Face_variance_map m_face_variance_map;
 
@@ -74,6 +74,15 @@ public:
   }
 
 public:
+  template <typename VertexPointMap>
+  Mat_4 construct_quadric_from_vertex(typename boost::graph_traits<TriangleMesh>::vertex_descriptor /*v*/,
+                                      const TriangleMesh& /*tmesh*/,
+                                      const VertexPointMap /*point_map*/,
+                                      const GeomTraits& /*gt*/) const
+  {
+    return Mat_4::Zero();
+  }
+
   // we don't have a sensible way to construct a triangle quadric
   // from an edge, so we fall back to probabilistic plane quadrics here
   template<typename VertexPointMap>

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

@@ -32,6 +32,15 @@ class Triangle_quadric_calculator
 public:
   Triangle_quadric_calculator() { }
 
+  template <typename VertexPointMap>
+  Mat_4 construct_quadric_from_vertex(typename boost::graph_traits<TriangleMesh>::vertex_descriptor /*v*/,
+                                      const TriangleMesh& /*tmesh*/,
+                                      const VertexPointMap /*point_map*/,
+                                      const GeomTraits& /*gt*/) const
+  {
+    return Mat_4::Zero();
+  }
+
   template <typename VertexPointMap>
   Mat_4 construct_quadric_from_edge(typename boost::graph_traits<TriangleMesh>::halfedge_descriptor he,
                                     const TriangleMesh& tmesh,

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

@@ -0,0 +1,144 @@
+// 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)     : Leo Valque
+
+#ifndef CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_GARLANDHECKBERT_COMPOSED_POLICIES_H
+#define CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_GARLANDHECKBERT_COMPOSED_POLICIES_H
+
+#include <CGAL/license/Surface_mesh_simplification.h>
+
+#include <CGAL/Surface_mesh_simplification/internal/Common.h>
+#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_policy_base.h>
+#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_functions.h>
+
+namespace CGAL {
+namespace Surface_mesh_simplification {
+namespace internal {
+
+/*
+This policy was created to combine the line policy with the existing ones.
+The composition of some quadrics is always invertible, but this is not necessarily true for all composed quadrics.
+`invertible` boolean parameter reflects this distinction.
+*/
+template <typename TriangleMesh, typename GeomTraits, typename Quadric_calculator_1, typename Quadric_calculator_2, bool invertible=false>
+class Composed_quadric_calculator
+{
+  typedef typename GarlandHeckbert_matrix_types<GeomTraits>::Mat_4             Mat_4;
+  typedef typename GarlandHeckbert_matrix_types<GeomTraits>::Col_4             Col_4;
+  typedef typename GarlandHeckbert_matrix_types<GeomTraits>::Row_4             Row_4;
+
+  Quadric_calculator_1 quadric_calculator_1;
+  Quadric_calculator_2 quadric_calculator_2;
+  double weight_1;
+  double weight_2;
+
+public:
+  Composed_quadric_calculator(const Quadric_calculator_1& qc1, const Quadric_calculator_2& qc2, double w1 = 1., double w2 = 1.)
+      : quadric_calculator_1(qc1)
+      , quadric_calculator_2(qc2)
+      , weight_1(w1)
+      , weight_2(w2) {}
+  Composed_quadric_calculator(double w1 = 1., double w2 = 1.)
+      : weight_1(w1)
+      , weight_2(w2) {}
+
+  template <typename VertexPointMap>
+  Mat_4 construct_quadric_from_vertex(typename boost::graph_traits<TriangleMesh>::vertex_descriptor v,
+                                      const TriangleMesh& tmesh,
+                                      const VertexPointMap point_map,
+                                      const GeomTraits& gt) const
+  {
+    return weight_1 * quadric_calculator_1.construct_quadric_from_vertex(v, tmesh, point_map, gt) +
+           weight_2 * quadric_calculator_2.construct_quadric_from_vertex(v, tmesh, point_map, gt);
+  }
+
+  template <typename VertexPointMap>
+  Mat_4 construct_quadric_from_edge(typename boost::graph_traits<TriangleMesh>::halfedge_descriptor he,
+                                    const TriangleMesh& tmesh,
+                                    const VertexPointMap point_map,
+                                    const GeomTraits& gt) const
+  {
+    return weight_1 * quadric_calculator_1.construct_quadric_from_edge(he, tmesh, point_map, gt) +
+           weight_2 * quadric_calculator_2.construct_quadric_from_edge(he, tmesh, point_map, gt);
+  }
+
+  template <typename VertexPointMap>
+  Mat_4 construct_quadric_from_face(typename boost::graph_traits<TriangleMesh>::face_descriptor f,
+                                    const TriangleMesh& tmesh,
+                                    const VertexPointMap point_map,
+                                    const GeomTraits& gt) const
+  {
+    return weight_1 * quadric_calculator_1.construct_quadric_from_face(f, tmesh, point_map, gt) +
+           weight_2 * quadric_calculator_2.construct_quadric_from_face(f, tmesh, point_map, gt);
+  }
+
+
+  Col_4 construct_optimal_point(const Mat_4& quadric,
+                                const Col_4& p0,
+                                const Col_4& p1) const
+  {
+    if constexpr(invertible)
+      return construct_optimal_point_invertible<GeomTraits>(quadric);
+    else
+      return construct_optimal_point_singular<GeomTraits>(quadric, p0, p1);
+  }
+};
+
+template<typename TriangleMesh, typename GeomTraits, typename GH_policies_1, typename GH_policies_2, bool invertible=false>
+class GarlandHeckbert_composed_policies
+  : public internal::GarlandHeckbert_cost_and_placement<
+             internal::Composed_quadric_calculator<TriangleMesh, GeomTraits,
+                                                   typename GH_policies_1::Quadric_calculator,
+                                                   typename GH_policies_2::Quadric_calculator,
+                                                   invertible>,
+             TriangleMesh, GeomTraits>
+{
+public:
+  typedef internal::Composed_quadric_calculator<TriangleMesh, GeomTraits,
+                                               typename GH_policies_1::Quadric_calculator,
+                                               typename GH_policies_2::Quadric_calculator,
+                                               invertible>
+          Quadric_calculator;
+
+private:
+  typedef internal::GarlandHeckbert_cost_and_placement<
+            Quadric_calculator, TriangleMesh, GeomTraits>                      Base;
+  typedef GarlandHeckbert_composed_policies<TriangleMesh, GeomTraits, GH_policies_1, GH_policies_2, invertible> Self;
+
+public:
+  typedef Self                                                                 Get_cost;
+  typedef Self                                                                 Get_placement;
+
+  typedef typename GeomTraits::FT                                              FT;
+
+public:
+  GarlandHeckbert_composed_policies(TriangleMesh& tmesh,
+                                    double w1=1., double w2=1.,const double dm = 100)
+    : Base(tmesh, Quadric_calculator(w1, w2), FT(dm))
+  { }
+
+  GarlandHeckbert_composed_policies(TriangleMesh& tmesh,
+                                    GH_policies_1 ghp1,
+                                    GH_policies_2 ghp2,
+                                    double w1=1., double w2=1.,const double dm = 100)
+    : Base(tmesh, Quadric_calculator(ghp1.quadric_calculator(), ghp2.quadric_calculator(), w1, w2), FT(dm))
+  { }
+
+public:
+  const Get_cost& get_cost() const { return *this; }
+  const Get_placement& get_placement() const { return *this; }
+
+  using Base::operator();
+};
+
+} // namespace internal
+} // namespace Surface_mesh_simplification
+} // namespace CGAL
+
+#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_GARLANDHECKBERT_COMPOSED_POLICIES_H

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

@@ -17,8 +17,13 @@
 #include <CGAL/Surface_mesh_simplification/internal/Common.h>
 #include <CGAL/boost/graph/helpers.h>
 
+#include <CGAL/Polygon_mesh_processing/compute_normal.h>
+
 #include <CGAL/Origin.h>
 
+#include <Eigen/Dense>
+#include <Eigen/Eigenvalues>
+
 #include <iostream>
 
 namespace CGAL {
@@ -236,8 +241,8 @@ construct_optimal_point_singular(const typename GarlandHeckbert_matrix_types<Geo
   }
   else
   {
+#if 1
     Col_4 opt_pt;
-
     const Col_4 p1mp0 = p1 - p0;
     const FT a = (p1mp0.transpose() * quadric * p1mp0)(0, 0);
     const FT b = 2 * (p0.transpose() * quadric * p1mp0)(0, 0);
@@ -272,6 +277,26 @@ construct_optimal_point_singular(const typename GarlandHeckbert_matrix_types<Geo
       }
     }
 
+#else
+    /*
+    Note:
+    A simpler code for the case of a non-invertible matrix
+    Experiments show that the results are identical to the version above in the majority of the cases, and nearly identical runtime
+    Therefore old case was kept.
+    */
+
+    // low rank -> svd pseudo-inverse
+    Eigen::JacobiSVD<Eigen::MatrixXd> svd_decomp(mat, Eigen::ComputeThinU | Eigen::ComputeThinV);
+    svd_decomp.setThreshold(1e-5);
+
+    opt_pt(0) = mean.x();
+    opt_pt(1) = mean.y();
+    opt_pt(2) = mean.z();
+    opt_pt(3) = 1.;
+
+    opt_pt += svd_decomp.solve(Col_4(0,0,0,1) - mat * opt_pt);
+#endif
+
     return opt_pt;
   }
 }
@@ -508,6 +533,45 @@ construct_prob_triangle_quadric_from_face(typename boost::graph_traits<TriangleM
   return ret;
 }
 
+////////////////////////////////////////////////////////////////////////////////////////////////////
+/// LINE QUADRICS (Liu et al. 2025)
+////////////////////////////////////////////////////////////////////////////////////////////////////
+
+template <typename GeomTraits>
+typename GarlandHeckbert_matrix_types<GeomTraits>::Mat_4
+construct_line_quadric_from_normal(const typename GeomTraits::Vector_3& normal,
+                                   const typename GeomTraits::Point_3& point,
+                                   const GeomTraits& gt)
+{
+  typedef typename GeomTraits::Vector_3 Vector_3;
+
+  auto cp = gt.construct_cross_product_vector_3_object();
+  auto plane = gt.construct_plane_3_object();
+  auto c_point = gt.construct_point_3_object();
+  auto base= gt.construct_base_vector_3_object();
+
+  Vector_3 x = base(plane(c_point(0,0,0),normal), 1);
+  CGAL::Polygon_mesh_processing::internal::normalize(x, gt);
+  Vector_3 y = cp(x,normal);
+
+  return construct_classic_plane_quadric_from_normal(x, point, gt)+construct_classic_plane_quadric_from_normal(y, point, gt);
+}
+
+/*
+template <typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
+typename GarlandHeckbert_matrix_types<GeomTraits>::Mat_4
+construct_line_quadric_from_vertex(const typename boost::graph_traits<TriangleMesh>::vertex_descriptor v,
+                                   const TriangleMesh& mesh,
+                                   const VertexPointMap vpm,
+                                   const GeomTraits& gt)
+{
+  typedef typename GeomTraits::Vector_3 Vector_3;
+
+  const Vector_3 normal = Polygon_mesh_processing::compute_vertex_normal(v, mesh, parameters::geom_traits(gt).vertex_point_map(vpm));
+  return construct_line_quadric_from_normal(normal, get(vpm, v), gt);
+}
+*/
+
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 /// PROB VARIANCE
 ////////////////////////////////////////////////////////////////////////////////////////////////////

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

@@ -0,0 +1,132 @@
+// 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)     : Leo Valque
+
+#ifndef CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_GARLANDHECKBERT_LINE_POLICIES_H
+#define CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_GARLANDHECKBERT_LINE_POLICIES_H
+
+#include <CGAL/license/Surface_mesh_simplification.h>
+
+#include <CGAL/Surface_mesh_simplification/internal/Common.h>
+#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_policy_base.h>
+#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_functions.h>
+
+#include <CGAL/Polygon_mesh_processing/compute_normal.h>
+
+#include <CGAL/Default.h>
+#include <boost/property_map/property_map.hpp>
+#include <boost/property_map/function_property_map.hpp>
+
+namespace CGAL {
+namespace Surface_mesh_simplification {
+namespace internal {
+
+/*
+This policy is not useful on its own; it is designed to be combined with another policy using a small weight.
+Therefore, it is kept internal.
+*/
+template <typename TriangleMesh, typename GeomTraits>
+class Line_quadric_calculator
+{
+  typedef typename GarlandHeckbert_matrix_types<GeomTraits>::Mat_4             Mat_4;
+  typedef typename GarlandHeckbert_matrix_types<GeomTraits>::Col_4             Col_4;
+  typedef typename GarlandHeckbert_matrix_types<GeomTraits>::Row_4             Row_4;
+
+  typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor        vertex_descriptor;
+  typedef typename GeomTraits::Vector_3                                        Vector_3;
+
+private:
+  std::function<Vector_3(vertex_descriptor)> m_vertex_normal_map;
+
+public:
+  Line_quadric_calculator() = delete;
+
+  template <typename VNM>
+  Line_quadric_calculator(const VNM vnm)
+      : m_vertex_normal_map([vnm](vertex_descriptor v) -> Vector_3{ return get(vnm, v); })
+  { }
+
+  template <typename VertexPointMap>
+  Mat_4 construct_quadric_from_vertex(typename boost::graph_traits<TriangleMesh>::vertex_descriptor v,
+                                      const TriangleMesh& /*tmesh*/,
+                                      const VertexPointMap point_map,
+                                      const GeomTraits& gt) const
+  {
+    return construct_line_quadric_from_normal(m_vertex_normal_map(v), get(point_map, v), gt);
+  }
+
+  template <typename VertexPointMap>
+  Mat_4 construct_quadric_from_edge(typename boost::graph_traits<TriangleMesh>::halfedge_descriptor /*he*/,
+                                    const TriangleMesh& /*tmesh*/,
+                                    const VertexPointMap /*point_map*/,
+                                    const GeomTraits& /*gt*/) const
+  {
+    return Mat_4::Zero();
+  }
+
+  template <typename VertexPointMap>
+  Mat_4 construct_quadric_from_face(typename boost::graph_traits<TriangleMesh>::face_descriptor /*f*/,
+                                    const TriangleMesh& /*tmesh*/,
+                                    const VertexPointMap /*point_map*/,
+                                    const GeomTraits& /*gt*/) const
+  {
+    return Mat_4::Zero();
+  }
+
+  /*
+  // Since these policies are never used alone but always composed with another one, this code is not used
+  Col_4 construct_optimal_point(const Mat_4& quadric,
+                                const Col_4& p0,
+                                const Col_4& p1) const
+  {
+    return construct_optimal_point_singular<GeomTraits>(quadric, p0, p1);
+  }
+  */
+};
+
+template<typename TriangleMesh,
+         typename GeomTraits,
+         typename VertexNormalMap = typename boost::property_map<TriangleMesh,
+                                                CGAL::dynamic_vertex_property_t<typename GeomTraits::Vector_3> >::type >
+class GarlandHeckbert_line_policies
+  : public internal::GarlandHeckbert_cost_and_placement<
+             internal::Line_quadric_calculator<TriangleMesh, GeomTraits>, TriangleMesh, GeomTraits>
+{
+public:
+  typedef internal::Line_quadric_calculator<TriangleMesh, GeomTraits> Quadric_calculator;
+
+private:
+  typedef internal::GarlandHeckbert_cost_and_placement<
+            Quadric_calculator, TriangleMesh, GeomTraits>                           Base;
+  typedef GarlandHeckbert_line_policies<TriangleMesh, GeomTraits, VertexNormalMap>  Self;
+
+public:
+  typedef Self                                                                      Get_cost;
+  typedef Self                                                                      Get_placement;
+
+  typedef typename GeomTraits::FT                                                   FT;
+
+public:
+  template<typename VNM>
+  GarlandHeckbert_line_policies(TriangleMesh& tmesh, const VNM vnm)
+    : Base(tmesh, Quadric_calculator(vnm), 0)
+  { }
+
+public:
+  const Get_cost& get_cost() const { return *this; }
+  const Get_placement& get_placement() const { return *this; }
+
+  using Base::operator();
+};
+
+} // namespace internal
+} // namespace Surface_mesh_simplification
+} // namespace CGAL
+
+#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_GARLANDHECKBERT_LINE_POLICIES_H

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

@@ -72,6 +72,10 @@ public:
   {
     m_cost_matrices = get(Cost_property(), tmesh);
   }
+
+  Quadric_calculator quadric_calculator() const{
+    return m_quadric_calculator;
+  }
 };
 
 template <typename QuadricCalculator, typename TriangleMesh, typename GeomTraits>
@@ -129,6 +133,15 @@ public:
   }
 
 public:
+  template <typename VertexPointMap>
+  Mat_4 construct_quadric(const vertex_descriptor v,
+                          const TriangleMesh& tmesh,
+                          const VertexPointMap vpm,
+                          const GeomTraits& gt) const
+  {
+    return quadric_calculator().construct_quadric_from_vertex(v, tmesh, vpm, gt);
+  }
+
   template <typename VertexPointMap>
   Mat_4 construct_quadric(const halfedge_descriptor he,
                           const TriangleMesh& tmesh,
@@ -154,10 +167,8 @@ public:
                   const VertexPointMap vpm,
                   const GeomTraits& gt) const
   {
-    Mat_4 zero_mat = Mat_4::Zero();
-
     for(vertex_descriptor v : vertices(tmesh))
-      put(vcm(), v, zero_mat);
+      put(vcm(), v, construct_quadric(v, tmesh, vpm, gt));
 
     for(face_descriptor f : faces(tmesh))
     {

Surface_mesh_simplification/package_info/Surface_mesh_simplification/dependencies

@@ -18,3 +18,4 @@ STL_Extension
 Spatial_searching
 Stream_support
 Surface_mesh_simplification
+Polygon_mesh_processing

Surface_mesh_simplification/test/Surface_mesh_simplification/edge_collapse_garland_heckbert_variations.cpp

@@ -32,6 +32,7 @@ typedef SMS::GarlandHeckbert_plane_policies<Surface_mesh, Kernel>
 typedef SMS::GarlandHeckbert_probabilistic_plane_policies<Surface_mesh, Kernel>    Prob_plane;
 typedef SMS::GarlandHeckbert_triangle_policies<Surface_mesh, Kernel>               Classic_tri;
 typedef SMS::GarlandHeckbert_probabilistic_triangle_policies<Surface_mesh, Kernel> Prob_tri;
+typedef SMS::GarlandHeckbert_plane_and_line_policies<Surface_mesh, Kernel>         Classic_plane_and_line;
 
 // settings for benchmarking - throw away the first n_burns results and keep the n_samples samples
 constexpr int n_burns = 1;
@@ -41,6 +42,7 @@ constexpr std::size_t classic_plane_index = 0;
 constexpr std::size_t prob_plane_index = 1;
 constexpr std::size_t classic_tri_index = 2;
 constexpr std::size_t prob_tri_index = 3;
+constexpr std::size_t plane_and_line_index = 4;
 
 // =================================================================================================
 // =================================================================================================
@@ -89,7 +91,8 @@ get_all_meshes(const std::vector<std::string>& filenames)
 
 template <typename Policy>
 Surface_mesh edge_collapse(Surface_mesh& mesh,
-                           const double ratio = 0.2)
+                           const double ratio,
+                           const Policy p)
 {
   typedef typename Policy::Get_cost Cost;
   typedef typename Policy::Get_placement Placement;
@@ -97,8 +100,6 @@ Surface_mesh edge_collapse(Surface_mesh& mesh,
 
   std::cout << "Edge collapse mesh of " << num_edges(mesh) << " edges. Policy: " << typeid(Policy).name() << std::endl;
 
-  const Policy p { mesh, 100 };
-
   const Cost& cost = p.get_cost();
   const Placement& unbounded_placement = p.get_placement();
   Bounded_placement bounded_placement(unbounded_placement);
@@ -117,6 +118,12 @@ Surface_mesh edge_collapse(Surface_mesh& mesh,
   return mesh;
 }
 
+template <typename Policy>
+Surface_mesh edge_collapse(Surface_mesh& mesh,
+                           const double ratio = 0.2){
+  return edge_collapse(mesh, ratio, Policy(mesh));
+}
+
 // =================================================================================================
 // =================================================================================================
 // =================================================================================================
@@ -145,21 +152,6 @@ void time_mesh(const TriangleMesh mesh,
       << "Average: " << elapsed_ns / n_samples << " (ms)" << std::endl;
 }
 
-template <typename TriangleMesh>
-void time_policy(const TriangleMesh& mesh,
-                 std::ostream& out,
-                 const std::string& policy)
-{
-  if(policy == "classic_plane")
-    time_mesh<Classic_plane>(mesh, out);
-  else if(policy == "classic_tri")
-    time_mesh<Classic_tri>(mesh, out);
-  else if(policy == "prob_plane")
-    time_mesh<Prob_plane>(mesh, out);
-  else if(policy == "prob_tri")
-    time_mesh<Prob_tri>(mesh, out);
-}
-
 template <typename TriangleMesh>
 void time_all_policies(const TriangleMesh& mesh,
                        std::ostream& out)
@@ -170,6 +162,7 @@ void time_all_policies(const TriangleMesh& mesh,
   time_mesh<Classic_tri>(mesh, out);
   time_mesh<Prob_plane>(mesh, out);
   time_mesh<Prob_tri>(mesh, out);
+  time_mesh<Classic_plane_and_line>(mesh, out);
 }
 
 // =================================================================================================
@@ -196,15 +189,16 @@ double hausdorff_error(const TriangleMesh& mesh,
 
 // calculate approximate Hausdorff errors for all different policies at the same decimation ratio
 template <typename TriangleMesh>
-std::array<FT, 4> hausdorff_errors(const TriangleMesh& mesh,
+std::array<FT, 5> hausdorff_errors(const TriangleMesh& mesh,
                                    double ratio)
 {
-  std::array<FT, 4> ret { {0, 0, 0, 0} };
+  std::array<FT, 5> ret { {0, 0, 0, 0, 0} };
 
   ret[classic_plane_index] = hausdorff_error<Classic_plane>(mesh, ratio);
   ret[prob_plane_index] = hausdorff_error<Prob_plane>(mesh, ratio);
   ret[classic_tri_index] = hausdorff_error<Classic_tri>(mesh, ratio);
   ret[prob_tri_index] = hausdorff_error<Prob_tri>(mesh, ratio);
+  ret[plane_and_line_index] = hausdorff_error<Classic_plane_and_line>(mesh, ratio);
 
   return ret;
 }
@@ -218,13 +212,14 @@ void hausdorff_errors(const TriangleMesh& mesh,
 
   for(InputIt it=begin; it!=end; ++it)
   {
-    std::array<double, 4> errs = hausdorff_errors(mesh, *it);
+    std::array<double, 5> errs = hausdorff_errors(mesh, *it);
     out << " -- Hausdorff error for collapse with ratio: " << *it << '\n';
 
     out << "classic plane: " << errs[classic_plane_index] << std::endl;
     out << "prob plane   : " << errs[prob_plane_index] << std::endl;
     out << "classic tri  : " << errs[classic_tri_index] << std::endl;
     out << "prob tri     : " << errs[prob_tri_index] << std::endl;
+    out << "plane + line : " << errs[plane_and_line_index] << std::endl;
   }
 }
 
@@ -252,11 +247,13 @@ void gather_face_aspect_ratio(const TriangleMesh& mesh,
   Surface_mesh pp = mesh;
   Surface_mesh ct = mesh;
   Surface_mesh pt = mesh;
+  Surface_mesh pl = mesh;
 
   edge_collapse<Classic_plane>(cp);
   edge_collapse<Prob_plane>(pp);
   edge_collapse<Classic_tri>(ct);
   edge_collapse<Prob_tri>(pt);
+  edge_collapse<Classic_plane_and_line>(pl);
 
   out << "Face aspect-ratio: classic plane\n";
   write_aspect_ratios(cp, out);
@@ -266,6 +263,8 @@ void gather_face_aspect_ratio(const TriangleMesh& mesh,
   write_aspect_ratios(ct, out);
   out << "Face aspect-ratio: prob triangle\n";
   write_aspect_ratios(pt, out);
+  out << "Face aspect-ratio: classic plane and line\n";
+  write_aspect_ratios(pl, out);
 }
 
 template <typename TriangleMesh>
@@ -287,6 +286,21 @@ void run(const std::pair<TriangleMesh, std::string>& input)
   gather_face_aspect_ratio(input.first, out);
 }
 
+template<typename TriangleMesh>
+void test_parameters_plane_and_line(const TriangleMesh& mesh){
+  using CGAL::Surface_mesh_simplification::make_GarlandHeckbert_plane_and_line_policies;
+  using PMap = boost::associative_property_map<std::map<typename boost::graph_traits<TriangleMesh>::vertex_descriptor, typename Kernel::Vector_3>>;
+  std::map<typename boost::graph_traits<TriangleMesh>::vertex_descriptor, typename Kernel::Vector_3> map;
+  PMap pmap(map);
+  CGAL::Polygon_mesh_processing::compute_vertex_normals(mesh, pmap);
+  TriangleMesh cp = mesh;
+  edge_collapse(cp, 0.2, make_GarlandHeckbert_plane_and_line_policies(cp,
+                                  CGAL::parameters::line_policies_weight(0.001)
+                                    .discontinuity_multiplier(50)
+                                    .geom_traits(Kernel())
+                                    .vertex_normal_map(pmap)));
+}
+
 int main(int argc, char** argv)
 {
   std::vector<std::string> default_data = { "data/helmet.off",
@@ -301,8 +315,10 @@ int main(int argc, char** argv)
     data = default_data;
 
   std::vector<std::pair<Surface_mesh, std::string> > named_meshes = get_all_meshes(data);
-  for(const auto& e : named_meshes)
+  for(const auto& e : named_meshes){
     run(e);
+    test_parameters_plane_and_line(e.first);
+  }
 
   std::cout << "Done!" << std::endl;