From fef9a2930b037356e807f38070558e2438c2c53d Mon Sep 17 00:00:00 2001
From: JM Etancelin <jean-matthieu.etancelin@univ-pau.fr>
Date: Thu, 9 Nov 2023 09:37:52 +0100
Subject: [PATCH 1/7] add image feature detection in mesh generation from inr.

---
 pygalmesh/main.py         |  7 +++-
 src/generate_from_inr.cpp | 77 ++++++++++++++++++++++++++++-----------
 src/generate_from_inr.hpp |  1 +
 src/pybind11.cpp          |  1 +
 4 files changed, 62 insertions(+), 24 deletions(-)

diff --git a/pygalmesh/main.py b/pygalmesh/main.py
index 50b3777..6540f10 100644
--- a/pygalmesh/main.py
+++ b/pygalmesh/main.py
@@ -24,7 +24,7 @@ class Wrapper(SizingFieldBase):
     def __init__(self, f):
         self.f = f
         super().__init__()
-
+ 
     def eval(self, x):
         return self.f(x)
 
@@ -295,6 +295,7 @@ def generate_from_inr(
     odt: bool = False,
     perturb: bool = True,
     exude: bool = True,
+    with_features: bool = False,
     max_edge_size_at_feature_edges: float = 0.0,
     min_facet_angle: float = 0.0,
     max_radius_surface_delaunay_ball: float = 0.0,
@@ -308,7 +309,6 @@ def generate_from_inr(
 ):
     fh, outfile = tempfile.mkstemp(suffix=".mesh")
     os.close(fh)
-
     if isinstance(max_cell_circumradius, float):
         _generate_from_inr(
             inr_filename,
@@ -317,6 +317,7 @@ def generate_from_inr(
             odt=odt,
             perturb=perturb,
             exude=exude,
+            with_features=with_features,
             max_edge_size_at_feature_edges=max_edge_size_at_feature_edges,
             min_facet_angle=min_facet_angle,
             max_radius_surface_delaunay_ball=max_radius_surface_delaunay_ball,
@@ -446,6 +447,7 @@ def generate_from_array(
     odt: bool = False,
     perturb: bool = True,
     exude: bool = True,
+    with_features: bool = False,
     max_edge_size_at_feature_edges: float = 0.0,
     min_facet_angle: float = 0.0,
     max_radius_surface_delaunay_ball: float = 0.0,
@@ -465,6 +467,7 @@ def generate_from_array(
         odt,
         perturb,
         exude,
+        with_features,
         max_edge_size_at_feature_edges,
         min_facet_angle,
         max_radius_surface_delaunay_ball,
diff --git a/src/generate_from_inr.cpp b/src/generate_from_inr.cpp
index 1d63fe6..0ce236b 100644
--- a/src/generate_from_inr.cpp
+++ b/src/generate_from_inr.cpp
@@ -14,12 +14,16 @@
 #include <CGAL/Implicit_mesh_domain_3.h>
 #include <CGAL/Mesh_domain_with_polyline_features_3.h>
 #include <CGAL/make_mesh_3.h>
+#include <CGAL/Labeled_mesh_domain_3.h>
+
+#include <CGAL/Mesh_3/Detect_features_in_image.h>
 
 namespace pygalmesh {
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
 
 typedef CGAL::Labeled_mesh_domain_3<K> Mesh_domain;
+typedef CGAL::Mesh_domain_with_polyline_features_3<Mesh_domain> Mesh_domain_with_features;
 
 // Triangulation
 typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr;
@@ -33,6 +37,7 @@ typedef Mesh_criteria::Cell_criteria Cell_criteria;
 typedef CGAL::Mesh_constant_domain_field_3<Mesh_domain::R,
                                            Mesh_domain::Index> Sizing_field_cell;
 
+
 void
 generate_from_inr(
     const std::string & inr_filename,
@@ -41,6 +46,7 @@ generate_from_inr(
     const bool odt,
     const bool perturb,
     const bool exude,
+    const bool with_features,
     const double max_edge_size_at_feature_edges,
     const double min_facet_angle,
     const double max_radius_surface_delaunay_ball,
@@ -60,7 +66,6 @@ generate_from_inr(
   if (!success) {
     throw "Could not read image file";
   }
-  Mesh_domain cgal_domain = Mesh_domain::create_labeled_image_mesh_domain(image);
 
   Mesh_criteria criteria(
       CGAL::parameters::edge_size=max_edge_size_at_feature_edges,
@@ -68,7 +73,8 @@ generate_from_inr(
       CGAL::parameters::facet_size=max_radius_surface_delaunay_ball,
       CGAL::parameters::facet_distance=max_facet_distance,
       CGAL::parameters::cell_radius_edge_ratio=max_circumradius_edge_ratio,
-      CGAL::parameters::cell_size=max_cell_circumradius
+      CGAL::parameters::cell_size=max_cell_circumradius,
+      CGAL::parameters::facet_topology = with_features ? CGAL::FACET_VERTICES_ON_SAME_SURFACE_PATCH : CGAL::FACET_VERTICES_ON_SURFACE
       );
 
   // Mesh generation
@@ -76,27 +82,54 @@ generate_from_inr(
     // suppress output
     std::cerr.setstate(std::ios_base::failbit);
   }
-  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(
-      cgal_domain,
-      criteria,
-      lloyd ? CGAL::parameters::lloyd() : CGAL::parameters::no_lloyd(),
-      odt ? CGAL::parameters::odt() : CGAL::parameters::no_odt(),
-      perturb ? CGAL::parameters::perturb() : CGAL::parameters::no_perturb(),
-      exude ?
-        CGAL::parameters::exude(
-          CGAL::parameters::time_limit = exude_time_limit,
-          CGAL::parameters::sliver_bound = exude_sliver_bound
-        ) :
-        CGAL::parameters::no_exude()
-      );
-  if (!verbose) {
-    std::cerr.clear();
-  }
 
-  // Output
-  std::ofstream medit_file(outfile);
-  c3t3.output_to_medit(medit_file);
-  medit_file.close();
+  if (with_features) {
+      C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(
+        Mesh_domain_with_features::create_labeled_image_mesh_domain(
+          image,
+          CGAL::parameters::features_detector = CGAL::Mesh_3::Detect_features_in_image()),
+        criteria,
+        lloyd ? CGAL::parameters::lloyd() : CGAL::parameters::no_lloyd(),
+        odt ? CGAL::parameters::odt() : CGAL::parameters::no_odt(),
+        perturb ? CGAL::parameters::perturb() : CGAL::parameters::no_perturb(),
+        exude ?
+          CGAL::parameters::exude(
+            CGAL::parameters::time_limit = exude_time_limit,
+            CGAL::parameters::sliver_bound = exude_sliver_bound
+            ) :
+          CGAL::parameters::no_exude()
+        );
+      if (!verbose) {
+          std::cerr.clear();
+      }
+
+      // Output
+      std::ofstream medit_file(outfile);
+      c3t3.output_to_medit(medit_file);
+      medit_file.close();
+  } else {
+      C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(
+        Mesh_domain::create_labeled_image_mesh_domain(image),
+        criteria,
+        lloyd ? CGAL::parameters::lloyd() : CGAL::parameters::no_lloyd(),
+        odt ? CGAL::parameters::odt() : CGAL::parameters::no_odt(),
+        perturb ? CGAL::parameters::perturb() : CGAL::parameters::no_perturb(),
+        exude ?
+          CGAL::parameters::exude(
+            CGAL::parameters::time_limit = exude_time_limit,
+            CGAL::parameters::sliver_bound = exude_sliver_bound
+            ) :
+          CGAL::parameters::no_exude()
+        );
+      if (!verbose) {
+          std::cerr.clear();
+      }
+
+      // Output
+      std::ofstream medit_file(outfile);
+      c3t3.output_to_medit(medit_file);
+      medit_file.close();
+  }
   return;
 }
 
diff --git a/src/generate_from_inr.hpp b/src/generate_from_inr.hpp
index 79b04b7..97f33a8 100644
--- a/src/generate_from_inr.hpp
+++ b/src/generate_from_inr.hpp
@@ -13,6 +13,7 @@ void generate_from_inr(
     const bool odt = false,
     const bool perturb = true,
     const bool exude = true,
+    const bool with_features = false,
     const double max_edge_size_at_feature_edges = 0.0,  // std::numeric_limits<double>::max(),
     const double min_facet_angle = 0.0,
     const double max_radius_surface_delaunay_ball = 0.0,
diff --git a/src/pybind11.cpp b/src/pybind11.cpp
index 9c8afc1..b0ea686 100644
--- a/src/pybind11.cpp
+++ b/src/pybind11.cpp
@@ -349,6 +349,7 @@ PYBIND11_MODULE(_pygalmesh, m) {
         py::arg("odt") = false,
         py::arg("perturb") = true,
         py::arg("exude") = true,
+        py::arg("with_features") = false,
         py::arg("max_edge_size_at_feature_edges") = 0.0,
         py::arg("min_facet_angle") = 0.0,
         py::arg("max_radius_surface_delaunay_ball") = 0.0,

From 49f2e3d4a9270ae53a81eb9ff49cc33e42524ebb Mon Sep 17 00:00:00 2001
From: JM Etancelin <jean-matthieu.etancelin@univ-pau.fr>
Date: Tue, 14 Nov 2023 17:25:05 +0100
Subject: [PATCH 2/7] cleanup

---
 src/generate_from_inr.cpp | 74 +++++++++++++++------------------------
 1 file changed, 28 insertions(+), 46 deletions(-)

diff --git a/src/generate_from_inr.cpp b/src/generate_from_inr.cpp
index 0ce236b..0a644e1 100644
--- a/src/generate_from_inr.cpp
+++ b/src/generate_from_inr.cpp
@@ -66,6 +66,14 @@ generate_from_inr(
   if (!success) {
     throw "Could not read image file";
   }
+  Mesh_domain_with_features *cgal_domain;
+  if (with_features)
+      cgal_domain = new Mesh_domain_with_features(
+          Mesh_domain_with_features::create_labeled_image_mesh_domain(
+              image,
+              CGAL::parameters::features_detector = CGAL::Mesh_3::Detect_features_in_image()));
+  else
+      cgal_domain = new Mesh_domain_with_features(Mesh_domain_with_features::create_labeled_image_mesh_domain(image));
 
   Mesh_criteria criteria(
       CGAL::parameters::edge_size=max_edge_size_at_feature_edges,
@@ -83,53 +91,27 @@ generate_from_inr(
     std::cerr.setstate(std::ios_base::failbit);
   }
 
-  if (with_features) {
-      C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(
-        Mesh_domain_with_features::create_labeled_image_mesh_domain(
-          image,
-          CGAL::parameters::features_detector = CGAL::Mesh_3::Detect_features_in_image()),
-        criteria,
-        lloyd ? CGAL::parameters::lloyd() : CGAL::parameters::no_lloyd(),
-        odt ? CGAL::parameters::odt() : CGAL::parameters::no_odt(),
-        perturb ? CGAL::parameters::perturb() : CGAL::parameters::no_perturb(),
-        exude ?
-          CGAL::parameters::exude(
-            CGAL::parameters::time_limit = exude_time_limit,
-            CGAL::parameters::sliver_bound = exude_sliver_bound
-            ) :
-          CGAL::parameters::no_exude()
-        );
-      if (!verbose) {
-          std::cerr.clear();
-      }
-
-      // Output
-      std::ofstream medit_file(outfile);
-      c3t3.output_to_medit(medit_file);
-      medit_file.close();
-  } else {
-      C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(
-        Mesh_domain::create_labeled_image_mesh_domain(image),
-        criteria,
-        lloyd ? CGAL::parameters::lloyd() : CGAL::parameters::no_lloyd(),
-        odt ? CGAL::parameters::odt() : CGAL::parameters::no_odt(),
-        perturb ? CGAL::parameters::perturb() : CGAL::parameters::no_perturb(),
-        exude ?
-          CGAL::parameters::exude(
-            CGAL::parameters::time_limit = exude_time_limit,
-            CGAL::parameters::sliver_bound = exude_sliver_bound
-            ) :
-          CGAL::parameters::no_exude()
-        );
-      if (!verbose) {
-          std::cerr.clear();
-      }
-
-      // Output
-      std::ofstream medit_file(outfile);
-      c3t3.output_to_medit(medit_file);
-      medit_file.close();
+  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(
+    *cgal_domain,
+    criteria,
+    lloyd ? CGAL::parameters::lloyd() : CGAL::parameters::no_lloyd(),
+    odt ? CGAL::parameters::odt() : CGAL::parameters::no_odt(),
+    perturb ? CGAL::parameters::perturb() : CGAL::parameters::no_perturb(),
+    exude ?
+      CGAL::parameters::exude(
+        CGAL::parameters::time_limit = exude_time_limit,
+        CGAL::parameters::sliver_bound = exude_sliver_bound
+        ) :
+      CGAL::parameters::no_exude()
+    );
+  if (!verbose) {
+      std::cerr.clear();
   }
+
+  // Output
+  std::ofstream medit_file(outfile);
+  c3t3.output_to_medit(medit_file);
+  medit_file.close();
   return;
 }
 

From 925721fc4096c0a52237f6b817a32e7f4ddda22a Mon Sep 17 00:00:00 2001
From: JM Etancelin <jean-matthieu.etancelin@univ-pau.fr>
Date: Thu, 11 Jan 2024 08:27:21 +0100
Subject: [PATCH 3/7] cleanup

---
 pygalmesh/main.py         | 3 ++-
 src/generate_from_inr.cpp | 3 ---
 2 files changed, 2 insertions(+), 4 deletions(-)

diff --git a/pygalmesh/main.py b/pygalmesh/main.py
index 6540f10..57d3650 100644
--- a/pygalmesh/main.py
+++ b/pygalmesh/main.py
@@ -24,7 +24,7 @@ class Wrapper(SizingFieldBase):
     def __init__(self, f):
         self.f = f
         super().__init__()
- 
+
     def eval(self, x):
         return self.f(x)
 
@@ -309,6 +309,7 @@ def generate_from_inr(
 ):
     fh, outfile = tempfile.mkstemp(suffix=".mesh")
     os.close(fh)
+
     if isinstance(max_cell_circumradius, float):
         _generate_from_inr(
             inr_filename,
diff --git a/src/generate_from_inr.cpp b/src/generate_from_inr.cpp
index 0a644e1..3507d3b 100644
--- a/src/generate_from_inr.cpp
+++ b/src/generate_from_inr.cpp
@@ -15,7 +15,6 @@
 #include <CGAL/Mesh_domain_with_polyline_features_3.h>
 #include <CGAL/make_mesh_3.h>
 #include <CGAL/Labeled_mesh_domain_3.h>
-
 #include <CGAL/Mesh_3/Detect_features_in_image.h>
 
 namespace pygalmesh {
@@ -37,7 +36,6 @@ typedef Mesh_criteria::Cell_criteria Cell_criteria;
 typedef CGAL::Mesh_constant_domain_field_3<Mesh_domain::R,
                                            Mesh_domain::Index> Sizing_field_cell;
 
-
 void
 generate_from_inr(
     const std::string & inr_filename,
@@ -90,7 +88,6 @@ generate_from_inr(
     // suppress output
     std::cerr.setstate(std::ios_base::failbit);
   }
-
   C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(
     *cgal_domain,
     criteria,

From 3dece32bc380f1622833b628fc7ab82b496b626a Mon Sep 17 00:00:00 2001
From: JM Etancelin <jean-matthieu.etancelin@univ-pau.fr>
Date: Fri, 15 Mar 2024 14:51:20 +0100
Subject: [PATCH 4/7] add relative error bound parameter for generate and
 generate_from_inr

---
 pygalmesh/main.py         | 4 ++++
 src/generate.cpp          | 4 +++-
 src/generate.hpp          | 1 +
 src/generate_from_inr.cpp | 5 +++--
 src/generate_from_inr.hpp | 1 +
 src/pybind11.cpp          | 2 ++
 6 files changed, 14 insertions(+), 3 deletions(-)

diff --git a/pygalmesh/main.py b/pygalmesh/main.py
index 57d3650..c3a3c9a 100644
--- a/pygalmesh/main.py
+++ b/pygalmesh/main.py
@@ -45,6 +45,7 @@ def generate_mesh(
     max_cell_circumradius: float | Callable[..., float] = 0.0,
     exude_time_limit: float = 0.0,
     exude_sliver_bound: float = 0.0,
+    relative_error_bound: float = 1e-3,
     verbose: bool = True,
     seed: int = 0,
 ):
@@ -125,6 +126,7 @@ def _select(obj):
         max_cell_circumradius_field=max_cell_circumradius_field,
         exude_time_limit=exude_time_limit,
         exude_sliver_bound=exude_sliver_bound,
+        relative_error_bound=relative_error_bound,
         verbose=verbose,
         seed=seed,
     )
@@ -304,6 +306,7 @@ def generate_from_inr(
     max_cell_circumradius: float | dict[int | str, float] = 0.0,
     exude_time_limit: float = 0.0,
     exude_sliver_bound: float = 0.0,
+    relative_error_bound: float = 1e-3,
     verbose: bool = True,
     seed: int = 0,
 ):
@@ -327,6 +330,7 @@ def generate_from_inr(
             max_cell_circumradius=max_cell_circumradius,
             exude_time_limit=exude_time_limit,
             exude_sliver_bound=exude_sliver_bound,
+            relative_error_bound=relative_error_bound,
             verbose=verbose,
             seed=seed,
         )
diff --git a/src/generate.cpp b/src/generate.cpp
index 64b0ba5..3b4ffa3 100644
--- a/src/generate.cpp
+++ b/src/generate.cpp
@@ -74,6 +74,7 @@ generate_mesh(
     //
     const double exude_time_limit,
     const double exude_sliver_bound,
+    const double relative_error_bound,
     //
     const bool verbose,
     const int seed
@@ -93,7 +94,8 @@ generate_mesh(
 
   Mesh_domain cgal_domain = Mesh_domain::create_implicit_mesh_domain(
        d,
-       K::Sphere_3(CGAL::ORIGIN, bounding_sphere_radius2)
+       K::Sphere_3(CGAL::ORIGIN, bounding_sphere_radius2),
+       CGAL::parameters::relative_error_bound(relative_error_bound)
        );
 
   // cgal_domain.detect_features();
diff --git a/src/generate.hpp b/src/generate.hpp
index eb188d3..daabe29 100644
--- a/src/generate.hpp
+++ b/src/generate.hpp
@@ -39,6 +39,7 @@ void generate_mesh(
     //
     const double exude_time_limit = 0.0,
     const double exude_sliver_bound = 0.0,
+    const double relative_error_bound = 1e-3,
     //
     const bool verbose = true,
     const int seed = 0
diff --git a/src/generate_from_inr.cpp b/src/generate_from_inr.cpp
index 3507d3b..6649adb 100644
--- a/src/generate_from_inr.cpp
+++ b/src/generate_from_inr.cpp
@@ -53,6 +53,7 @@ generate_from_inr(
     const double max_cell_circumradius,
     const double exude_time_limit,
     const double exude_sliver_bound,
+    const double relative_error_bound,
     const bool verbose,
     const int seed
     )
@@ -69,9 +70,9 @@ generate_from_inr(
       cgal_domain = new Mesh_domain_with_features(
           Mesh_domain_with_features::create_labeled_image_mesh_domain(
               image,
-              CGAL::parameters::features_detector = CGAL::Mesh_3::Detect_features_in_image()));
+              CGAL::parameters::features_detector(CGAL::Mesh_3::Detect_features_in_image()).relative_error_bound(relative_error_bound)));
   else
-      cgal_domain = new Mesh_domain_with_features(Mesh_domain_with_features::create_labeled_image_mesh_domain(image));
+      cgal_domain = new Mesh_domain_with_features(Mesh_domain_with_features::create_labeled_image_mesh_domain(image, CGAL::parameters::relative_error_bound(relative_error_bound)));
 
   Mesh_criteria criteria(
       CGAL::parameters::edge_size=max_edge_size_at_feature_edges,
diff --git a/src/generate_from_inr.hpp b/src/generate_from_inr.hpp
index 97f33a8..ba4f094 100644
--- a/src/generate_from_inr.hpp
+++ b/src/generate_from_inr.hpp
@@ -22,6 +22,7 @@ void generate_from_inr(
     const double max_cell_circumradius = 0.0,
     const double exude_time_limit = 0.0,
     const double exude_sliver_bound = 0.0,
+    const double relative_error_bound = 1e-3,
     const bool verbose = true,
     const int seed = 0
     );
diff --git a/src/pybind11.cpp b/src/pybind11.cpp
index b0ea686..e71a311 100644
--- a/src/pybind11.cpp
+++ b/src/pybind11.cpp
@@ -288,6 +288,7 @@ PYBIND11_MODULE(_pygalmesh, m) {
         py::arg("max_cell_circumradius_field") = nullptr,
         py::arg("exude_time_limit") = 0.0,
         py::arg("exude_sliver_bound") = 0.0,
+        py::arg("relative_error_bound") = 1e-3,
         py::arg("verbose") = true,
         py::arg("seed") = 0
         );
@@ -358,6 +359,7 @@ PYBIND11_MODULE(_pygalmesh, m) {
         py::arg("max_cell_circumradius") = 0.0,
         py::arg("exude_time_limit") = 0.0,
         py::arg("exude_sliver_bound") = 0.0,
+        py::arg("relative_error_bound") = 1e-3,
         py::arg("verbose") = true,
         py::arg("seed") = 0
         );

From b9250d76c68ed4944a4e6de87c5affa7cea35d43 Mon Sep 17 00:00:00 2001
From: JM Etancelin <jean-matthieu.etancelin@univ-pau.fr>
Date: Fri, 15 Mar 2024 15:11:08 +0100
Subject: [PATCH 5/7] fixup missing inr usages

---
 pygalmesh/main.py         | 3 +++
 src/generate_from_inr.cpp | 3 ++-
 src/generate_from_inr.hpp | 1 +
 src/pybind11.cpp          | 1 +
 4 files changed, 7 insertions(+), 1 deletion(-)

diff --git a/pygalmesh/main.py b/pygalmesh/main.py
index c3a3c9a..4fff08a 100644
--- a/pygalmesh/main.py
+++ b/pygalmesh/main.py
@@ -359,6 +359,7 @@ def generate_from_inr(
             max_radius_surface_delaunay_ball=max_radius_surface_delaunay_ball,
             max_facet_distance=max_facet_distance,
             max_circumradius_edge_ratio=max_circumradius_edge_ratio,
+            relative_error_bound=relative_error_bound,
             verbose=verbose,
             seed=seed,
         )
@@ -459,6 +460,7 @@ def generate_from_array(
     max_cell_circumradius: float | dict[int | str, float] = 0.0,
     max_facet_distance: float = 0.0,
     max_circumradius_edge_ratio: float = 0.0,
+    relative_error_bound: float = 1e-3,
     verbose: bool = True,
     seed: int = 0,
 ):
@@ -479,6 +481,7 @@ def generate_from_array(
         max_facet_distance,
         max_circumradius_edge_ratio,
         max_cell_circumradius,
+        relative_error_bound,
         verbose,
         seed,
     )
diff --git a/src/generate_from_inr.cpp b/src/generate_from_inr.cpp
index 6649adb..da2cf14 100644
--- a/src/generate_from_inr.cpp
+++ b/src/generate_from_inr.cpp
@@ -132,6 +132,7 @@ generate_from_inr_with_subdomain_sizing(
     const double max_circumradius_edge_ratio,
     const double exude_time_limit,
     const double exude_sliver_bound,
+    const double relative_error_bound,
     const bool verbose,
     const int seed
     )
@@ -143,7 +144,7 @@ generate_from_inr_with_subdomain_sizing(
   if (!success) {
     throw "Could not read image file";
   }
-  Mesh_domain cgal_domain = Mesh_domain::create_labeled_image_mesh_domain(image);
+  Mesh_domain cgal_domain = Mesh_domain::create_labeled_image_mesh_domain(image, CGAL::parameters::relative_error_bound(relative_error_bound));
 
   Sizing_field_cell max_cell_circumradius(default_max_cell_circumradius);
   const int ndimensions = 3;
diff --git a/src/generate_from_inr.hpp b/src/generate_from_inr.hpp
index ba4f094..5a3f684 100644
--- a/src/generate_from_inr.hpp
+++ b/src/generate_from_inr.hpp
@@ -45,6 +45,7 @@ generate_from_inr_with_subdomain_sizing(
     const double max_circumradius_edge_ratio = 0.0,
     const double exude_time_limit = 0.0,
     const double exude_sliver_bound = 0.0,
+    const double relative_error_bound = 1e-3,
     const bool verbose = true,
     const int seed = 0
     );
diff --git a/src/pybind11.cpp b/src/pybind11.cpp
index e71a311..7e4c35d 100644
--- a/src/pybind11.cpp
+++ b/src/pybind11.cpp
@@ -381,6 +381,7 @@ PYBIND11_MODULE(_pygalmesh, m) {
         py::arg("max_circumradius_edge_ratio") = 0.0,
         py::arg("exude_time_limit") = 0.0,
         py::arg("exude_sliver_bound") = 0.0,
+        py::arg("relative_error_bound") = 1e-3,
         py::arg("verbose") = true,
         py::arg("seed") = 0
         );

From 2c128b53a8956a1f0c29c36ababa6d94dfccd08b Mon Sep 17 00:00:00 2001
From: JM Etancelin <jean-matthieu.etancelin@univ-pau.fr>
Date: Wed, 10 Apr 2024 11:19:50 +0200
Subject: [PATCH 6/7] Add gray image with sigle iso-value mesh generation

---
 pygalmesh/main.py         | 39 +++++++++++++++++++++++----------------
 src/generate_from_inr.cpp | 34 +++++++++++++++++++++++++++-------
 src/generate_from_inr.hpp |  3 +++
 src/pybind11.cpp          |  2 ++
 4 files changed, 55 insertions(+), 23 deletions(-)

diff --git a/pygalmesh/main.py b/pygalmesh/main.py
index 4fff08a..652980e 100644
--- a/pygalmesh/main.py
+++ b/pygalmesh/main.py
@@ -307,6 +307,8 @@ def generate_from_inr(
     exude_time_limit: float = 0.0,
     exude_sliver_bound: float = 0.0,
     relative_error_bound: float = 1e-3,
+    iso_value: float = np.NAN,
+    value_outside: float = 0.,
     verbose: bool = True,
     seed: int = 0,
 ):
@@ -331,6 +333,7 @@ def generate_from_inr(
             exude_time_limit=exude_time_limit,
             exude_sliver_bound=exude_sliver_bound,
             relative_error_bound=relative_error_bound,
+            iso_value=iso_value,
             verbose=verbose,
             seed=seed,
         )
@@ -461,29 +464,33 @@ def generate_from_array(
     max_facet_distance: float = 0.0,
     max_circumradius_edge_ratio: float = 0.0,
     relative_error_bound: float = 1e-3,
+    iso_value: float = np.NAN,
+    value_outside: float = 0.,
     verbose: bool = True,
     seed: int = 0,
 ):
-    assert vol.dtype in ["uint8", "uint16"]
+    assert vol.dtype in ["uint8", "uint16", "float64", "float32"]
     fh, inr_filename = tempfile.mkstemp(suffix=".inr")
     os.close(fh)
     save_inr(vol, voxel_size, inr_filename)
     mesh = generate_from_inr(
-        inr_filename,
-        lloyd,
-        odt,
-        perturb,
-        exude,
-        with_features,
-        max_edge_size_at_feature_edges,
-        min_facet_angle,
-        max_radius_surface_delaunay_ball,
-        max_facet_distance,
-        max_circumradius_edge_ratio,
-        max_cell_circumradius,
-        relative_error_bound,
-        verbose,
-        seed,
+        inr_filename= inr_filename,
+        lloyd=lloyd,
+        odt=odt,
+        perturb=perturb,
+        exude=exude,
+        with_features=with_features,
+        max_edge_size_at_feature_edges=max_edge_size_at_feature_edges,
+        min_facet_angle=min_facet_angle,
+        max_radius_surface_delaunay_ball=max_radius_surface_delaunay_ball,
+        max_facet_distance=max_facet_distance,
+        max_circumradius_edge_ratio=max_circumradius_edge_ratio,
+        max_cell_circumradius=max_cell_circumradius,
+        relative_error_bound=relative_error_bound,
+        iso_value=iso_value,
+        value_outside=value_outside,
+        verbose=verbose,
+        seed=seed,
     )
     os.remove(inr_filename)
     return mesh
diff --git a/src/generate_from_inr.cpp b/src/generate_from_inr.cpp
index da2cf14..0b9d6b0 100644
--- a/src/generate_from_inr.cpp
+++ b/src/generate_from_inr.cpp
@@ -16,6 +16,7 @@
 #include <CGAL/make_mesh_3.h>
 #include <CGAL/Labeled_mesh_domain_3.h>
 #include <CGAL/Mesh_3/Detect_features_in_image.h>
+#include <CGAL/Mesh_3/Detect_features_on_image_bbox.h>
 
 namespace pygalmesh {
 
@@ -54,6 +55,8 @@ generate_from_inr(
     const double exude_time_limit,
     const double exude_sliver_bound,
     const double relative_error_bound,
+    const double iso_value,
+    const double value_outside,
     const bool verbose,
     const int seed
     )
@@ -65,14 +68,31 @@ generate_from_inr(
   if (!success) {
     throw "Could not read image file";
   }
+
   Mesh_domain_with_features *cgal_domain;
-  if (with_features)
-      cgal_domain = new Mesh_domain_with_features(
-          Mesh_domain_with_features::create_labeled_image_mesh_domain(
-              image,
-              CGAL::parameters::features_detector(CGAL::Mesh_3::Detect_features_in_image()).relative_error_bound(relative_error_bound)));
-  else
-      cgal_domain = new Mesh_domain_with_features(Mesh_domain_with_features::create_labeled_image_mesh_domain(image, CGAL::parameters::relative_error_bound(relative_error_bound)));
+  if (std::isnan(iso_value)) {
+      if (with_features)
+          cgal_domain = new Mesh_domain_with_features(
+              Mesh_domain_with_features::create_labeled_image_mesh_domain(
+                  image,
+                  CGAL::parameters::features_detector(CGAL::Mesh_3::Detect_features_in_image()).relative_error_bound(relative_error_bound)));
+      else {
+          cgal_domain = new Mesh_domain_with_features(Mesh_domain_with_features::create_labeled_image_mesh_domain(
+                  image,
+                  CGAL::parameters::relative_error_bound(relative_error_bound)));
+      }
+  }
+  else {
+      if (with_features)
+          cgal_domain = new Mesh_domain_with_features(
+              Mesh_domain_with_features::create_gray_image_mesh_domain(
+                  image,
+                  CGAL::parameters::features_detector(CGAL::Mesh_3::Detect_features_on_image_bbox()).relative_error_bound(relative_error_bound).iso_value(iso_value).value_outside(value_outside)));
+      else {
+          cgal_domain =
+              new Mesh_domain_with_features(Mesh_domain::create_gray_image_mesh_domain(image, CGAL::parameters::iso_value(iso_value).value_outside(value_outside)));
+      }
+  }
 
   Mesh_criteria criteria(
       CGAL::parameters::edge_size=max_edge_size_at_feature_edges,
diff --git a/src/generate_from_inr.hpp b/src/generate_from_inr.hpp
index 5a3f684..a8fa2cc 100644
--- a/src/generate_from_inr.hpp
+++ b/src/generate_from_inr.hpp
@@ -3,6 +3,7 @@
 
 #include <string>
 #include <vector>
+#include <cmath>
 
 namespace pygalmesh {
 
@@ -23,6 +24,8 @@ void generate_from_inr(
     const double exude_time_limit = 0.0,
     const double exude_sliver_bound = 0.0,
     const double relative_error_bound = 1e-3,
+    const double iso_value = std::nan(""),
+    const double value_outside = 1e-3,
     const bool verbose = true,
     const int seed = 0
     );
diff --git a/src/pybind11.cpp b/src/pybind11.cpp
index 7e4c35d..d4857dd 100644
--- a/src/pybind11.cpp
+++ b/src/pybind11.cpp
@@ -360,6 +360,8 @@ PYBIND11_MODULE(_pygalmesh, m) {
         py::arg("exude_time_limit") = 0.0,
         py::arg("exude_sliver_bound") = 0.0,
         py::arg("relative_error_bound") = 1e-3,
+        py::arg("iso_value") = std::nan(""),
+        py::arg("value_outside") = 0.,
         py::arg("verbose") = true,
         py::arg("seed") = 0
         );

From f6d74c817ee80577b11f5249b681a4f91a53cdb0 Mon Sep 17 00:00:00 2001
From: JM Etancelin <jean-matthieu.etancelin@univ-pau.fr>
Date: Tue, 21 Jan 2025 15:21:51 +0100
Subject: [PATCH 7/7] add multiple iso-values capabilities from array (and inr)

---
 pygalmesh/main.py         | 13 +++++++---
 src/generate_from_inr.cpp | 53 ++++++++++++++++++++++++++++++++++++---
 src/generate_from_inr.hpp |  3 ++-
 src/pybind11.cpp          |  5 ++--
 4 files changed, 63 insertions(+), 11 deletions(-)

diff --git a/pygalmesh/main.py b/pygalmesh/main.py
index 652980e..06d116d 100644
--- a/pygalmesh/main.py
+++ b/pygalmesh/main.py
@@ -307,18 +307,23 @@ def generate_from_inr(
     exude_time_limit: float = 0.0,
     exude_sliver_bound: float = 0.0,
     relative_error_bound: float = 1e-3,
-    iso_value: float = np.NAN,
+    iso_value: float | tuple = np.NAN,
     value_outside: float = 0.,
     verbose: bool = True,
     seed: int = 0,
 ):
     fh, outfile = tempfile.mkstemp(suffix=".mesh")
     os.close(fh)
-
+    if isinstance(iso_value, float):
+        if np.isnan(iso_value):
+            iso_value = ()
+        else:
+            iso_value = (iso_value, )
     if isinstance(max_cell_circumradius, float):
         _generate_from_inr(
             inr_filename,
             outfile,
+            list(iso_value),
             lloyd=lloyd,
             odt=odt,
             perturb=perturb,
@@ -333,7 +338,7 @@ def generate_from_inr(
             exude_time_limit=exude_time_limit,
             exude_sliver_bound=exude_sliver_bound,
             relative_error_bound=relative_error_bound,
-            iso_value=iso_value,
+            niso_value=len(iso_value),
             verbose=verbose,
             seed=seed,
         )
@@ -464,7 +469,7 @@ def generate_from_array(
     max_facet_distance: float = 0.0,
     max_circumradius_edge_ratio: float = 0.0,
     relative_error_bound: float = 1e-3,
-    iso_value: float = np.NAN,
+    iso_value: float | tuple = np.NAN,
     value_outside: float = 0.,
     verbose: bool = True,
     seed: int = 0,
diff --git a/src/generate_from_inr.cpp b/src/generate_from_inr.cpp
index 0b9d6b0..953dd06 100644
--- a/src/generate_from_inr.cpp
+++ b/src/generate_from_inr.cpp
@@ -34,6 +34,22 @@ typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
 typedef Mesh_criteria::Facet_criteria Facet_criteria;
 typedef Mesh_criteria::Cell_criteria Cell_criteria;
 
+class ImageValueToSubdomainIdx {
+public:
+    // Constructor that takes a boost::flat_set<float>
+    ImageValueToSubdomainIdx(const boost::container::flat_set<float>& inputSet) : mySet(inputSet) {}
+
+    // Overload the () operator
+    int operator()(double v) const {
+        return int(std::distance(mySet.begin(),
+                                 mySet.lower_bound(float(v))));
+    }
+
+private:
+    boost::container::flat_set<float> mySet; // The only member of the class
+};
+
+
 typedef CGAL::Mesh_constant_domain_field_3<Mesh_domain::R,
                                            Mesh_domain::Index> Sizing_field_cell;
 
@@ -41,6 +57,7 @@ void
 generate_from_inr(
     const std::string & inr_filename,
     const std::string & outfile,
+    const std::vector<double> & iso_value,
     const bool lloyd,
     const bool odt,
     const bool perturb,
@@ -55,7 +72,7 @@ generate_from_inr(
     const double exude_time_limit,
     const double exude_sliver_bound,
     const double relative_error_bound,
-    const double iso_value,
+    const int niso_value,
     const double value_outside,
     const bool verbose,
     const int seed
@@ -70,7 +87,7 @@ generate_from_inr(
   }
 
   Mesh_domain_with_features *cgal_domain;
-  if (std::isnan(iso_value)) {
+  if (niso_value==0) {
       if (with_features)
           cgal_domain = new Mesh_domain_with_features(
               Mesh_domain_with_features::create_labeled_image_mesh_domain(
@@ -82,15 +99,43 @@ generate_from_inr(
                   CGAL::parameters::relative_error_bound(relative_error_bound)));
       }
   }
+  else if(niso_value==1) {
+      if (with_features)
+          cgal_domain = new Mesh_domain_with_features(
+              Mesh_domain_with_features::create_gray_image_mesh_domain(
+                  image,
+                  CGAL::parameters::features_detector(CGAL::Mesh_3::Detect_features_on_image_bbox()).
+                  relative_error_bound(relative_error_bound).
+                  iso_value(iso_value[0]).
+                  value_outside(value_outside)));
+      else {
+          cgal_domain =
+              new Mesh_domain_with_features(Mesh_domain::create_gray_image_mesh_domain(
+                  image,
+                  CGAL::parameters::iso_value(iso_value[0]).
+                  value_outside(value_outside)));
+      }
+  }
   else {
+      typedef boost::container::flat_set<float> Flat_set;
+      Flat_set iso_values_set;
+      for(int i=0;i<niso_value;i++)
+          iso_values_set.insert(static_cast<float>(iso_value[i]));
+      ImageValueToSubdomainIdx f(iso_values_set);
       if (with_features)
           cgal_domain = new Mesh_domain_with_features(
               Mesh_domain_with_features::create_gray_image_mesh_domain(
                   image,
-                  CGAL::parameters::features_detector(CGAL::Mesh_3::Detect_features_on_image_bbox()).relative_error_bound(relative_error_bound).iso_value(iso_value).value_outside(value_outside)));
+                  CGAL::parameters::features_detector(CGAL::Mesh_3::Detect_features_on_image_bbox()).
+                  relative_error_bound(relative_error_bound).
+                  image_values_to_subdomain_indices(f).
+                  value_outside(value_outside)));
       else {
           cgal_domain =
-              new Mesh_domain_with_features(Mesh_domain::create_gray_image_mesh_domain(image, CGAL::parameters::iso_value(iso_value).value_outside(value_outside)));
+              new Mesh_domain_with_features(Mesh_domain::create_gray_image_mesh_domain(
+                  image,
+                  CGAL::parameters::image_values_to_subdomain_indices(f).
+                  value_outside(value_outside)));
       }
   }
 
diff --git a/src/generate_from_inr.hpp b/src/generate_from_inr.hpp
index a8fa2cc..455174d 100644
--- a/src/generate_from_inr.hpp
+++ b/src/generate_from_inr.hpp
@@ -10,6 +10,7 @@ namespace pygalmesh {
 void generate_from_inr(
     const std::string & inr_filename,
     const std::string & outfile,
+    const std::vector<double> & iso_value,
     const bool lloyd = false,
     const bool odt = false,
     const bool perturb = true,
@@ -24,7 +25,7 @@ void generate_from_inr(
     const double exude_time_limit = 0.0,
     const double exude_sliver_bound = 0.0,
     const double relative_error_bound = 1e-3,
-    const double iso_value = std::nan(""),
+    const int niso_value = 0,
     const double value_outside = 1e-3,
     const bool verbose = true,
     const int seed = 0
diff --git a/src/pybind11.cpp b/src/pybind11.cpp
index d4857dd..00c1cc7 100644
--- a/src/pybind11.cpp
+++ b/src/pybind11.cpp
@@ -346,6 +346,7 @@ PYBIND11_MODULE(_pygalmesh, m) {
         "_generate_from_inr", &generate_from_inr,
         py::arg("inr_filename"),
         py::arg("outfile"),
+        py::arg("iso_value"),
         py::arg("lloyd") = false,
         py::arg("odt") = false,
         py::arg("perturb") = true,
@@ -360,8 +361,8 @@ PYBIND11_MODULE(_pygalmesh, m) {
         py::arg("exude_time_limit") = 0.0,
         py::arg("exude_sliver_bound") = 0.0,
         py::arg("relative_error_bound") = 1e-3,
-        py::arg("iso_value") = std::nan(""),
-        py::arg("value_outside") = 0.,
+        py::arg("niso_value") = 0,
+        py::arg("value_outside") = 1e-3,
         py::arg("verbose") = true,
         py::arg("seed") = 0
         );