diff --git a/pygalmesh/main.py b/pygalmesh/main.py index 50b3777..06d116d 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, ) @@ -295,6 +297,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, @@ -303,20 +306,29 @@ 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, + 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, 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, @@ -325,6 +337,8 @@ 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, + niso_value=len(iso_value), verbose=verbose, seed=seed, ) @@ -353,6 +367,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, ) @@ -446,33 +461,41 @@ 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, 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, + iso_value: float | tuple = 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, - max_edge_size_at_feature_edges, - min_facet_angle, - max_radius_surface_delaunay_ball, - max_facet_distance, - max_circumradius_edge_ratio, - max_cell_circumradius, - 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.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 1d63fe6..953dd06 100644 --- a/src/generate_from_inr.cpp +++ b/src/generate_from_inr.cpp @@ -14,12 +14,16 @@ #include #include #include +#include +#include +#include namespace pygalmesh { typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef CGAL::Labeled_mesh_domain_3 Mesh_domain; +typedef CGAL::Mesh_domain_with_polyline_features_3 Mesh_domain_with_features; // Triangulation typedef CGAL::Mesh_triangulation_3::type Tr; @@ -30,6 +34,22 @@ typedef CGAL::Mesh_criteria_3 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 + ImageValueToSubdomainIdx(const boost::container::flat_set& 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 mySet; // The only member of the class +}; + + typedef CGAL::Mesh_constant_domain_field_3 Sizing_field_cell; @@ -37,10 +57,12 @@ void generate_from_inr( const std::string & inr_filename, const std::string & outfile, + const std::vector & iso_value, const bool lloyd, 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, @@ -49,6 +71,9 @@ generate_from_inr( const double max_cell_circumradius, const double exude_time_limit, const double exude_sliver_bound, + const double relative_error_bound, + const int niso_value, + const double value_outside, const bool verbose, const int seed ) @@ -60,7 +85,59 @@ generate_from_inr( if (!success) { throw "Could not read image file"; } - Mesh_domain cgal_domain = Mesh_domain::create_labeled_image_mesh_domain(image); + + Mesh_domain_with_features *cgal_domain; + if (niso_value==0) { + 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(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 Flat_set; + Flat_set iso_values_set; + for(int i=0;i(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). + 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::image_values_to_subdomain_indices(f). + value_outside(value_outside))); + } + } Mesh_criteria criteria( CGAL::parameters::edge_size=max_edge_size_at_feature_edges, @@ -68,7 +145,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 @@ -77,20 +155,20 @@ generate_from_inr( std::cerr.setstate(std::ios_base::failbit); } C3t3 c3t3 = CGAL::make_mesh_3( - 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_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() - ); + CGAL::parameters::no_exude() + ); if (!verbose) { - std::cerr.clear(); + std::cerr.clear(); } // Output @@ -119,6 +197,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 ) @@ -130,7 +209,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 79b04b7..455174d 100644 --- a/src/generate_from_inr.hpp +++ b/src/generate_from_inr.hpp @@ -3,16 +3,19 @@ #include #include +#include namespace pygalmesh { void generate_from_inr( const std::string & inr_filename, const std::string & outfile, + const std::vector & iso_value, const bool lloyd = false, 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::max(), const double min_facet_angle = 0.0, const double max_radius_surface_delaunay_ball = 0.0, @@ -21,6 +24,9 @@ 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 int niso_value = 0, + const double value_outside = 1e-3, const bool verbose = true, const int seed = 0 ); @@ -43,6 +49,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 9c8afc1..00c1cc7 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 ); @@ -345,10 +346,12 @@ 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, 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, @@ -357,6 +360,9 @@ 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("niso_value") = 0, + py::arg("value_outside") = 1e-3, py::arg("verbose") = true, py::arg("seed") = 0 ); @@ -378,6 +384,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 );