generate_from_inr_feature_input function included, reads polylines file

pygalmesh/main.py

@@ -12,6 +12,7 @@
     _generate_2d,
     _generate_from_inr,
     _generate_from_inr_with_subdomain_sizing,
+    _generate_from_inr_feature_input,
     _generate_from_off,
     _generate_mesh,
     _generate_periodic_mesh,
@@ -291,6 +292,7 @@ def generate_volume_mesh_from_surface_mesh(
 
 def generate_from_inr(
     inr_filename: str,
+    feature_filename: str = None,
     lloyd: bool = False,
     odt: bool = False,
     perturb: bool = True,
@@ -309,59 +311,84 @@ def generate_from_inr(
     fh, outfile = tempfile.mkstemp(suffix=".mesh")
     os.close(fh)
 
-    if isinstance(max_cell_circumradius, float):
-        _generate_from_inr(
-            inr_filename,
-            outfile,
-            lloyd=lloyd,
-            odt=odt,
-            perturb=perturb,
-            exude=exude,
-            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,
-            exude_time_limit=exude_time_limit,
-            exude_sliver_bound=exude_sliver_bound,
-            verbose=verbose,
-            seed=seed,
-        )
-    else:
-        assert isinstance(max_cell_circumradius, dict)
-        if "default" in max_cell_circumradius.keys():
-            default_max_cell_circumradius = max_cell_circumradius.pop("default")
+    if feature_filename is None:
+        if isinstance(max_cell_circumradius, float):
+            _generate_from_inr(
+                inr_filename,
+                outfile,
+                lloyd=lloyd,
+                odt=odt,
+                perturb=perturb,
+                exude=exude,
+                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,
+                exude_time_limit=exude_time_limit,
+                exude_sliver_bound=exude_sliver_bound,
+                verbose=verbose,
+                seed=seed,
+            )
         else:
-            default_max_cell_circumradius = 0.0
-
-        max_cell_circumradiuss = list(max_cell_circumradius.values())
-        subdomain_labels = list(max_cell_circumradius.keys())
-
-        _generate_from_inr_with_subdomain_sizing(
-            inr_filename,
-            outfile,
-            default_max_cell_circumradius,
-            max_cell_circumradiuss,
-            subdomain_labels,
-            lloyd=lloyd,
-            odt=odt,
-            perturb=perturb,
-            exude=exude,
-            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,
-            verbose=verbose,
-            seed=seed,
-        )
+            assert isinstance(max_cell_circumradius, dict)
+            if "default" in max_cell_circumradius.keys():
+                default_max_cell_circumradius = max_cell_circumradius.pop("default")
+            else:
+                default_max_cell_circumradius = 0.0
+
+            max_cell_circumradiuss = list(max_cell_circumradius.values())
+            subdomain_labels = list(max_cell_circumradius.keys())
+
+            _generate_from_inr_with_subdomain_sizing(
+                inr_filename,
+                outfile,
+                default_max_cell_circumradius,
+                max_cell_circumradiuss,
+                subdomain_labels,
+                lloyd=lloyd,
+                odt=odt,
+                perturb=perturb,
+                exude=exude,
+                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,
+                verbose=verbose,
+                seed=seed,
+            )
 
+    else:
+        if isinstance(max_cell_circumradius, float):
+            _generate_from_inr_feature_input(
+                inr_filename,
+                feature_filename=feature_filename,
+                outfile=outfile,
+                lloyd=lloyd,
+                odt=odt,
+                perturb=perturb,
+                exude=exude,
+                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,
+                exude_time_limit=exude_time_limit,
+                exude_sliver_bound=exude_sliver_bound,
+                verbose=verbose,
+                seed=seed,
+            )
+
+        else:
+            raise ValueError("max_cell_circumradius must be a float")
+
     mesh = meshio.read(outfile)
     os.remove(outfile)
     return mesh
 
-
 def remesh_surface(
     filename: str,
     max_edge_size_at_feature_edges: float = 0.0,
@@ -442,6 +469,7 @@ def save_inr(vol, voxel_size: tuple[float, float, float], fname: str):
 def generate_from_array(
     vol,
     voxel_size: tuple[float, float, float],
+    feature_filename: str = None,
     lloyd: bool = False,
     odt: bool = False,
     perturb: bool = True,
@@ -461,6 +489,7 @@ def generate_from_array(
     save_inr(vol, voxel_size, inr_filename)
     mesh = generate_from_inr(
         inr_filename,
+        feature_filename,
         lloyd,
         odt,
         perturb,

src/generate_from_inr.cpp

@@ -1,6 +1,7 @@
 #define CGAL_MESH_3_VERBOSE 1
 
 #include "generate_from_inr.hpp"
+#include "read_polylines.h"
 
 #include <cassert>
 
@@ -19,7 +20,8 @@ namespace pygalmesh {
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
 
-typedef CGAL::Labeled_mesh_domain_3<K> Mesh_domain;
+typedef CGAL::Labeled_mesh_domain_3<K> Image_domain;
+typedef CGAL::Mesh_domain_with_polyline_features_3<Image_domain> Mesh_domain;
 
 // Triangulation
 typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr;
@@ -175,4 +177,94 @@ generate_from_inr_with_subdomain_sizing(
   return;
 }
 
+void
+generate_from_inr_feature_input(
+    const std::string & inr_filename,
+    const std::string & feature_filename,
+    const std::string & outfile,
+    const bool lloyd,
+    const bool odt,
+    const bool perturb,
+    const bool exude,
+    const double max_edge_size_at_feature_edges,
+    const double min_facet_angle,
+    const double max_radius_surface_delaunay_ball,
+    const double max_facet_distance,
+    const double max_circumradius_edge_ratio,
+    const double max_cell_circumradius,
+    const double exude_time_limit,
+    const double exude_sliver_bound,
+    const bool verbose,
+    const int seed
+    )
+{
+  CGAL::get_default_random() = CGAL::Random(seed);
+
+  CGAL::Image_3 image;
+  const bool success = image.read(inr_filename.c_str());
+  if (!success) {
+    throw "Could not read image file";
+  }
+
+  const std::string lines_fname = CGAL::data_file_path(feature_filename);
+  using Point_3 = K::Point_3;
+  std::vector<std::vector<Point_3>> features_inside;
+
+  if (!read_polylines(lines_fname, features_inside)) // see file "read_polylines.h"
+  {
+    std::cerr << "Error: Cannot read file " << lines_fname << std::endl;
+  }
+
+  Mesh_domain cgal_domain =
+    Mesh_domain::create_labeled_image_mesh_domain(image, CGAL::parameters::input_features = std::cref(features_inside));
+
+  CGAL::Bbox_3 bbox = cgal_domain.bbox();
+  double diag = CGAL::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +
+                           CGAL::square(bbox.ymax() - bbox.ymin()) +
+                           CGAL::square(bbox.zmax() - bbox.zmin()));
+
+  double sizing;
+  if (max_edge_size_at_feature_edges < 0.001) {
+    sizing = diag * 0.05;
+  } else {
+    sizing = max_edge_size_at_feature_edges;
+  }
+
+  Mesh_criteria criteria(
+      CGAL::parameters::edge_size=sizing,
+      CGAL::parameters::facet_angle=min_facet_angle,
+      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
+      );
+
+  // Mesh generation
+  if (!verbose) {
+    // 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::time_limit(0)) : CGAL::parameters::no_lloyd(),
+      odt ? CGAL::parameters::odt(CGAL::parameters::time_limit(0)) : 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;
+}
 } // namespace pygalmesh

src/generate_from_inr.hpp

@@ -47,6 +47,27 @@ generate_from_inr_with_subdomain_sizing(
     const int seed = 0
     );
 
+void
+generate_from_inr_feature_input(
+    const std::string & inr_filename,
+    const std::string & feature_filename,
+    const std::string & outfile,
+    const bool lloyd = false,
+    const bool odt = false,
+    const bool perturb = true,
+    const bool exude = true,
+    const double max_edge_size_at_feature_edges = 0.0,
+    const double min_facet_angle = 0.0,
+    const double max_radius_surface_delaunay_ball = 0.0,
+    const double max_facet_distance = 0.0,
+    const double max_circumradius_edge_ratio = 0.0,
+    const double max_cell_circumradius = 0.0,
+    const double exude_time_limit = 0.0,
+    const double exude_sliver_bound = 0.0,
+    const bool verbose = true,
+    const int seed = 0
+    );
+
 } // namespace pygalmesh
 
 #endif // GENERATE_FROM_INR_HPP

src/pybind11.cpp

@@ -381,6 +381,26 @@ PYBIND11_MODULE(_pygalmesh, m) {
         py::arg("verbose") = true,
         py::arg("seed") = 0
         );
+    m.def(
+        "_generate_from_inr_feature_input", &generate_from_inr_feature_input,
+        py::arg("inr_filename"),
+        py::arg("feature_filename"),
+        py::arg("outfile"),
+        py::arg("lloyd") = false,
+        py::arg("odt") = false,
+        py::arg("perturb") = true,
+        py::arg("exude") = true,
+        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,
+        py::arg("max_facet_distance") = 0.0,
+        py::arg("max_circumradius_edge_ratio") = 0.0,
+        py::arg("max_cell_circumradius") = 0.0,
+        py::arg("exude_time_limit") = 0.0,
+        py::arg("exude_sliver_bound") = 0.0,
+        py::arg("verbose") = true,
+        py::arg("seed") = 0
+        );
     m.def(
         "_remesh_surface", &remesh_surface,
         py::arg("infile"),

src/read_polylines.h

@@ -0,0 +1,29 @@
+#ifndef READ_POLYLINES_H
+#define READ_POLYLINES_H
+
+#include <cstddef>
+#include <vector>
+#include <fstream>
+
+template <typename Point_3>
+bool read_polylines(const std::string fname,
+                    std::vector<std::vector<Point_3> >& polylines)
+{
+  std::ifstream ifs(fname);
+  if(ifs.bad()) return false;
+  std::size_t n;
+  while(ifs >> n) {
+    polylines.resize(polylines.size()+1);
+    std::vector<Point_3>& polyline = polylines.back();
+    while(n-- != 0) {
+      Point_3 p;
+      ifs >> p;
+      if(ifs.fail()) return false;
+      polyline.push_back(p);
+    }
+  }
+  if(ifs.bad()) return false;
+  else return ifs.eof();
+}
+
+#endif // READ_POLYLINES_H