diff --git a/MIGRATION_GUIDE.md b/MIGRATION_GUIDE.md index 87500a9..de9f8aa 100644 --- a/MIGRATION_GUIDE.md +++ b/MIGRATION_GUIDE.md @@ -1933,6 +1933,47 @@ Important details: See `development/mesh/check_marching_cubes.py` for reference comparisons and `development/mesh/benchmark_marching_cubes.py` for reproducible timings. +### Mesh Smoothing + +Laplacian mesh smoothing is available as `bic.mesh.smooth_mesh`. It was moved +from `bic.utils.smooth_mesh`; the old location is no longer available. The C++ +API likewise moved from `#include "bioimage_cpp/mesh_smoothing.hxx"` and +`bioimage_cpp::smooth_mesh` to `#include "bioimage_cpp/mesh/smoothing.hxx"` +and `bioimage_cpp::mesh::smooth_mesh`. + +### Mesh Simplification + +Topology-preserving triangle-mesh simplification is available as +`bic.mesh.simplify_mesh`. It uses quadric-error edge collapses, preserves +manifold topology and open boundaries, and applies configurable soft penalties +to sharp features: + + vertices, faces, _, values = bic.mesh.marching_cubes( + volume, level=0.5, allow_degenerate=False + ) + vertices, faces, normals, values = bic.mesh.simplify_mesh( + vertices, + faces, + reduction=0.8, + values=values, + feature_angle=45.0, + feature_weight=10.0, + ) + +`reduction` is the approximate fraction of input faces to remove. A mesh can +stop above its target if no topology- and boundary-safe collapse remains. +Normals are recomputed from the final face geometry. Optional scalar values are +interpolated along collapsed edges; the fourth return item is `None` when +values are omitted. Vertices and values preserve float32/float64 input dtypes, +while returned faces are int64. + +The input must be a non-empty, consistently oriented triangle 2-manifold with +no duplicate, zero-area, or unreferenced geometry. Use +`marching_cubes(..., allow_degenerate=False)` before simplification; invalid +meshes are rejected rather than silently repaired. + + + ### Anti-Aliased Resampling `affine_transform` itself never pre-smooths the input; downsampling without diff --git a/README.md b/README.md index 73b6be8..6ec893a 100644 --- a/README.md +++ b/README.md @@ -8,7 +8,7 @@ Image processing and segmentation functionality in C++ with light-weight python bindings through nanobind and minimal dependencies to enable distribution via pip. The package includes dependency-free triangle-mesh extraction from 3D volumes -and segmentation masks under `bioimage_cpp.mesh`. +and segmentation masks, plus Laplacian mesh smoothing, under `bioimage_cpp.mesh`. The `bioimage_cpp` python library can be installed via pip: ```bash diff --git a/development/utils/_mesh_smoothing_reference.py b/development/mesh/_mesh_smoothing_reference.py similarity index 100% rename from development/utils/_mesh_smoothing_reference.py rename to development/mesh/_mesh_smoothing_reference.py diff --git a/development/utils/benchmark_mesh_smoothing.py b/development/mesh/benchmark_mesh_smoothing.py similarity index 93% rename from development/utils/benchmark_mesh_smoothing.py rename to development/mesh/benchmark_mesh_smoothing.py index 3dd0d6a..de39daa 100644 --- a/development/utils/benchmark_mesh_smoothing.py +++ b/development/mesh/benchmark_mesh_smoothing.py @@ -8,8 +8,8 @@ Run:: - python development/utils/benchmark_mesh_smoothing.py --n-vertices 5000 --iterations 5 --repeats 3 - python development/utils/benchmark_mesh_smoothing.py --n-vertices 20000 --threads 1,2,4,0 + python development/mesh/benchmark_mesh_smoothing.py --n-vertices 5000 --iterations 5 --repeats 3 + python development/mesh/benchmark_mesh_smoothing.py --n-vertices 20000 --threads 1,2,4,0 Notes ----- @@ -117,7 +117,7 @@ def main() -> int: label = f"bioimage_cpp[n_threads={n_threads}]" def run(verts=mesh.verts, normals=mesh.normals, faces=mesh.faces, n=n_threads): - bic.utils.smooth_mesh(verts, normals, faces, iterations=args.iterations, n_threads=n) + bic.mesh.smooth_mesh(verts, normals, faces, iterations=args.iterations, n_threads=n) times = time_call(run, args.repeats, args.warmup) rows.append( @@ -141,7 +141,7 @@ def run_ref(): # exact agreement; see module docstring). if reference is not None: ref_v, ref_n = reference(mesh.verts, mesh.normals, mesh.faces, 1) - ours_v, ours_n = bic.utils.smooth_mesh(mesh.verts, mesh.normals, mesh.faces, iterations=1) + ours_v, ours_n = bic.mesh.smooth_mesh(mesh.verts, mesh.normals, mesh.faces, iterations=1) v_max_diff = float(np.max(np.abs(ours_v - ref_v))) n_max_diff = float(np.max(np.abs(ours_n - ref_n))) print(f"Correctness vs reference (iterations=1): " diff --git a/include/bioimage_cpp/mesh/simplification.hxx b/include/bioimage_cpp/mesh/simplification.hxx new file mode 100644 index 0000000..cdb9eb2 --- /dev/null +++ b/include/bioimage_cpp/mesh/simplification.hxx @@ -0,0 +1,869 @@ +#pragma once + +#include "bioimage_cpp/array_view.hxx" +#include "bioimage_cpp/detail/edge_hash.hxx" +#include "bioimage_cpp/detail/indexed_heap.hxx" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace bioimage_cpp::mesh { + +template +struct SimplifyMeshResult { + std::vector vertices; + std::vector faces; + std::vector normals; + std::optional> values; +}; + +namespace detail::simplification { + +using NodeId = std::uint64_t; +using FaceId = std::size_t; +using Edge = bioimage_cpp::detail::Edge; +using EdgeHash = bioimage_cpp::detail::EdgeHash; +using Vec3 = std::array; +using Triangle = std::array; + +inline Vec3 subtract(const Vec3 &a, const Vec3 &b) { + return {a[0] - b[0], a[1] - b[1], a[2] - b[2]}; +} + +inline Vec3 add_scaled(const Vec3 &a, const Vec3 &b, const double t) { + return {a[0] + t * b[0], a[1] + t * b[1], a[2] + t * b[2]}; +} + +inline double dot(const Vec3 &a, const Vec3 &b) { + return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; +} + +inline Vec3 cross(const Vec3 &a, const Vec3 &b) { + return { + a[1] * b[2] - a[2] * b[1], + a[2] * b[0] - a[0] * b[2], + a[0] * b[1] - a[1] * b[0], + }; +} + +inline double norm(const Vec3 &v) { + return std::sqrt(dot(v, v)); +} + +struct Quadric { + std::array matrix{}; + + Quadric &operator+=(const Quadric &other) { + for (std::size_t i = 0; i < matrix.size(); ++i) { + matrix[i] += other.matrix[i]; + } + return *this; + } + + void add_scaled(const Quadric &other, const double scale) { + for (std::size_t i = 0; i < matrix.size(); ++i) { + matrix[i] += scale * other.matrix[i]; + } + } + + [[nodiscard]] double evaluate(const Vec3 &point) const { + const std::array h{point[0], point[1], point[2], 1.0}; + double value = 0.0; + for (std::size_t row = 0; row < 4; ++row) { + for (std::size_t col = 0; col < 4; ++col) { + value += h[row] * matrix[row * 4 + col] * h[col]; + } + } + return value; + } + + [[nodiscard]] double bilinear( + const std::array &a, + const std::array &b + ) const { + double value = 0.0; + for (std::size_t row = 0; row < 4; ++row) { + for (std::size_t col = 0; col < 4; ++col) { + value += a[row] * matrix[row * 4 + col] * b[col]; + } + } + return value; + } +}; + +inline Quadric operator+(Quadric lhs, const Quadric &rhs) { + lhs += rhs; + return lhs; +} + +inline Quadric plane_quadric(const Vec3 &unit_normal, const double d, const double weight) { + const std::array plane{unit_normal[0], unit_normal[1], unit_normal[2], d}; + Quadric quadric; + for (std::size_t row = 0; row < 4; ++row) { + for (std::size_t col = 0; col < 4; ++col) { + quadric.matrix[row * 4 + col] = weight * plane[row] * plane[col]; + } + } + return quadric; +} + +inline Triangle triangle_key(Triangle triangle) { + std::sort(triangle.begin(), triangle.end()); + return triangle; +} + +struct TriangleHash { + std::size_t operator()(const Triangle &triangle) const { + std::size_t seed = 0; + for (const auto value : triangle) { + const auto item = static_cast(value); + seed ^= item + 0x9e3779b97f4a7c15ULL + (seed << 6U) + (seed >> 2U); + } + return seed; + } +}; + +struct Face { + Triangle vertices{}; + bool alive = true; +}; + +struct Vertex { + Vec3 position{}; + Quadric quadric; + double value = 0.0; + bool alive = true; + std::unordered_set faces; + std::unordered_set neighbours; +}; + +struct EdgeRecord { + std::array faces{}; + std::size_t count = 0; +}; + +struct FaceGeometry { + Vec3 normal{}; + double area = 0.0; + Quadric quadric; +}; + +struct Candidate { + Vec3 position{}; + double cost = 0.0; + double interpolation = 0.0; +}; + +inline bool solve_3x3(const Quadric &quadric, Vec3 &solution) { + double augmented[3][4] = { + {quadric.matrix[0], quadric.matrix[1], quadric.matrix[2], -quadric.matrix[3]}, + {quadric.matrix[4], quadric.matrix[5], quadric.matrix[6], -quadric.matrix[7]}, + {quadric.matrix[8], quadric.matrix[9], quadric.matrix[10], -quadric.matrix[11]}, + }; + double max_entry = 0.0; + for (const auto &row : augmented) { + for (std::size_t col = 0; col < 3; ++col) { + max_entry = std::max(max_entry, std::abs(row[col])); + } + } + const double tolerance = 64.0 * std::numeric_limits::epsilon() * max_entry; + if (max_entry == 0.0) { + return false; + } + + for (std::size_t col = 0; col < 3; ++col) { + std::size_t pivot = col; + for (std::size_t row = col + 1; row < 3; ++row) { + if (std::abs(augmented[row][col]) > std::abs(augmented[pivot][col])) { + pivot = row; + } + } + if (std::abs(augmented[pivot][col]) <= tolerance) { + return false; + } + if (pivot != col) { + for (std::size_t entry = col; entry < 4; ++entry) { + std::swap(augmented[col][entry], augmented[pivot][entry]); + } + } + const double divisor = augmented[col][col]; + for (std::size_t entry = col; entry < 4; ++entry) { + augmented[col][entry] /= divisor; + } + for (std::size_t row = 0; row < 3; ++row) { + if (row == col) continue; + const double factor = augmented[row][col]; + for (std::size_t entry = col; entry < 4; ++entry) { + augmented[row][entry] -= factor * augmented[col][entry]; + } + } + } + solution = {augmented[0][3], augmented[1][3], augmented[2][3]}; + return std::isfinite(solution[0]) && std::isfinite(solution[1]) && std::isfinite(solution[2]); +} + +template +class Simplifier { +public: + Simplifier( + const ConstArrayView &vertices, + const ConstArrayView &faces, + const ConstArrayView *values, + const double feature_angle, + const double feature_weight + ) : has_values_(values != nullptr), feature_angle_(feature_angle), feature_weight_(feature_weight) { + initialize(vertices, faces, values); + } + + SimplifyMeshResult run(const std::size_t target_faces) { + build_heap(); + while (alive_faces_ > target_faces && !heap_.empty()) { + const Edge edge = heap_.pop().key; + const auto candidate = make_candidate(edge); + if (!candidate.has_value() || !collapse_is_valid(edge, candidate->position)) { + continue; + } + collapse(edge, *candidate); + } + return finalize(); + } + +private: + using Priority = std::tuple; + using Heap = bioimage_cpp::detail::SparseIndexedHeap< + Edge, Priority, EdgeHash, std::greater + >; + + std::vector vertices_; + std::vector faces_; + std::vector face_geometry_; + std::unordered_map edges_; + std::unordered_map face_lookup_; + Heap heap_; + std::size_t alive_faces_ = 0; + bool has_values_ = false; + double feature_angle_ = 45.0; + double feature_weight_ = 10.0; + double area_tolerance_ = 0.0; + + void initialize( + const ConstArrayView &vertices, + const ConstArrayView &faces, + const ConstArrayView *values + ) { + if (vertices.shape.size() != 2 || vertices.shape[1] != 3) { + throw std::invalid_argument("vertices must have shape (n_vertices, 3)"); + } + if (faces.shape.size() != 2 || faces.shape[1] != 3) { + throw std::invalid_argument("faces must have shape (n_faces, 3)"); + } + const std::size_t n_vertices = static_cast(vertices.shape[0]); + const std::size_t n_faces = static_cast(faces.shape[0]); + if (n_vertices == 0 || n_faces == 0) { + throw std::invalid_argument("vertices and faces must describe a non-empty mesh"); + } + if (values != nullptr + && (values->shape.size() != 1 + || static_cast(values->shape[0]) != n_vertices)) { + throw std::invalid_argument("values must have shape (n_vertices,)"); + } + + vertices_.resize(n_vertices); + Vec3 lower{ + std::numeric_limits::infinity(), + std::numeric_limits::infinity(), + std::numeric_limits::infinity(), + }; + Vec3 upper{ + -std::numeric_limits::infinity(), + -std::numeric_limits::infinity(), + -std::numeric_limits::infinity(), + }; + for (std::size_t vertex = 0; vertex < n_vertices; ++vertex) { + for (std::size_t axis = 0; axis < 3; ++axis) { + const double coordinate = static_cast(vertices.data[vertex * 3 + axis]); + if (!std::isfinite(coordinate)) { + throw std::invalid_argument("vertices must contain only finite values"); + } + vertices_[vertex].position[axis] = coordinate; + lower[axis] = std::min(lower[axis], coordinate); + upper[axis] = std::max(upper[axis], coordinate); + } + if (values != nullptr) { + const double value = static_cast(values->data[vertex]); + if (!std::isfinite(value)) { + throw std::invalid_argument("values must contain only finite values"); + } + vertices_[vertex].value = value; + } + } + const Vec3 diagonal = subtract(upper, lower); + const double scale_squared = dot(diagonal, diagonal); + if (scale_squared == 0.0) { + throw std::invalid_argument("vertices must span a non-zero bounding box"); + } + area_tolerance_ = 64.0 * std::numeric_limits::epsilon() * scale_squared; + + faces_.resize(n_faces); + face_geometry_.resize(n_faces); + std::vector referenced(n_vertices, false); + edges_.reserve(n_faces * 3 / 2 + 1); + face_lookup_.reserve(n_faces); + for (std::size_t face_id = 0; face_id < n_faces; ++face_id) { + Triangle triangle{}; + for (std::size_t corner = 0; corner < 3; ++corner) { + const auto raw = faces.data[face_id * 3 + corner]; + if constexpr (std::is_signed_v) { + if (raw < 0) { + throw std::invalid_argument("faces contains a negative vertex index"); + } + } + const auto vertex = static_cast(raw); + if (vertex >= n_vertices) { + throw std::invalid_argument("faces contains a vertex index outside [0, n_vertices)"); + } + triangle[corner] = vertex; + referenced[static_cast(vertex)] = true; + } + if (triangle[0] == triangle[1] || triangle[0] == triangle[2] + || triangle[1] == triangle[2]) { + throw std::invalid_argument("faces must not contain repeated vertex indices"); + } + const Triangle key = triangle_key(triangle); + if (!face_lookup_.emplace(key, face_id).second) { + throw std::invalid_argument("faces must not contain duplicate triangles"); + } + faces_[face_id].vertices = triangle; + vertices_[static_cast(triangle[0])].faces.insert(face_id); + vertices_[static_cast(triangle[1])].faces.insert(face_id); + vertices_[static_cast(triangle[2])].faces.insert(face_id); + face_geometry_[face_id] = geometry_of(triangle); + if (face_geometry_[face_id].area * 2.0 <= area_tolerance_) { + throw std::invalid_argument("faces must not contain zero-area triangles"); + } + for (const auto vertex : triangle) { + vertices_[static_cast(vertex)].quadric += face_geometry_[face_id].quadric; + } + add_edge(triangle[0], triangle[1], face_id); + add_edge(triangle[1], triangle[2], face_id); + add_edge(triangle[2], triangle[0], face_id); + } + if (std::find(referenced.begin(), referenced.end(), false) != referenced.end()) { + throw std::invalid_argument("every vertex must be referenced by at least one face"); + } + alive_faces_ = n_faces; + validate_edge_orientations(); + validate_vertex_links(); + add_feature_quadrics(); + } + + FaceGeometry geometry_of(const Triangle &triangle) const { + const Vec3 &a = vertices_[static_cast(triangle[0])].position; + const Vec3 &b = vertices_[static_cast(triangle[1])].position; + const Vec3 &c = vertices_[static_cast(triangle[2])].position; + const Vec3 raw_normal = cross(subtract(b, a), subtract(c, a)); + const double length = norm(raw_normal); + FaceGeometry geometry; + geometry.area = 0.5 * length; + if (length > 0.0) { + geometry.normal = { + raw_normal[0] / length, + raw_normal[1] / length, + raw_normal[2] / length, + }; + const double d = -dot(geometry.normal, a); + geometry.quadric = plane_quadric(geometry.normal, d, geometry.area); + } + return geometry; + } + + void add_edge(const NodeId first, const NodeId second, const FaceId face_id) { + const Edge key = bioimage_cpp::detail::edge_key(first, second); + auto &record = edges_[key]; + if (record.count == 2) { + throw std::invalid_argument("faces must form a 2-manifold: an edge has more than two faces"); + } + record.faces[record.count++] = face_id; + vertices_[static_cast(first)].neighbours.insert(second); + vertices_[static_cast(second)].neighbours.insert(first); + } + + void remove_edge(const NodeId first, const NodeId second, const FaceId face_id) { + const Edge key = bioimage_cpp::detail::edge_key(first, second); + auto edge_it = edges_.find(key); + if (edge_it == edges_.end()) return; + auto &record = edge_it->second; + std::size_t position = record.count; + for (std::size_t index = 0; index < record.count; ++index) { + if (record.faces[index] == face_id) position = index; + } + if (position == record.count) return; + record.faces[position] = record.faces[record.count - 1]; + --record.count; + if (record.count == 0) { + edges_.erase(edge_it); + vertices_[static_cast(first)].neighbours.erase(second); + vertices_[static_cast(second)].neighbours.erase(first); + } + } + + bool face_traverses_edge_forward(const Face &face, const Edge &edge) const { + for (std::size_t corner = 0; corner < 3; ++corner) { + const NodeId first = face.vertices[corner]; + const NodeId second = face.vertices[(corner + 1) % 3]; + if (first == edge.first && second == edge.second) return true; + if (first == edge.second && second == edge.first) return false; + } + return false; + } + + void validate_edge_orientations() const { + for (const auto &[edge, record] : edges_) { + if (record.count == 2) { + const bool first = face_traverses_edge_forward(faces_[record.faces[0]], edge); + const bool second = face_traverses_edge_forward(faces_[record.faces[1]], edge); + if (first == second) { + throw std::invalid_argument( + "faces must have consistent winding across shared edges" + ); + } + } + } + } + + bool is_boundary_vertex(const NodeId vertex) const { + for (const NodeId neighbour : vertices_[static_cast(vertex)].neighbours) { + const auto it = edges_.find(bioimage_cpp::detail::edge_key(vertex, neighbour)); + if (it != edges_.end() && it->second.count == 1) return true; + } + return false; + } + + void validate_vertex_links() const { + for (NodeId vertex = 0; vertex < vertices_.size(); ++vertex) { + std::unordered_map> link; + for (const FaceId face_id : vertices_[static_cast(vertex)].faces) { + const auto &triangle = faces_[face_id].vertices; + NodeId first = 0; + NodeId second = 0; + bool found_first = false; + for (const NodeId item : triangle) { + if (item == vertex) continue; + if (!found_first) { + first = item; + found_first = true; + } else { + second = item; + } + } + link[first].push_back(second); + link[second].push_back(first); + } + std::size_t degree_one = 0; + for (const auto &[_, adjacent] : link) { + if (adjacent.size() == 1) ++degree_one; + else if (adjacent.size() != 2) { + throw std::invalid_argument("faces must form a manifold around every vertex"); + } + } + const bool boundary = is_boundary_vertex(vertex); + if ((!boundary && degree_one != 0) || (boundary && degree_one != 2)) { + throw std::invalid_argument("faces must form one manifold fan around every vertex"); + } + std::unordered_set visited; + std::vector pending{link.begin()->first}; + while (!pending.empty()) { + const NodeId current = pending.back(); + pending.pop_back(); + if (!visited.insert(current).second) continue; + for (const NodeId adjacent : link.at(current)) pending.push_back(adjacent); + } + if (visited.size() != link.size()) { + throw std::invalid_argument("faces must form one manifold fan around every vertex"); + } + } + } + + void add_feature_quadrics() { + if (feature_weight_ == 0.0) return; + const double radians = feature_angle_ * std::acos(-1.0) / 180.0; + const double threshold = std::cos(radians); + for (const auto &[edge, record] : edges_) { + if (record.count != 2) continue; + const auto &first = face_geometry_[record.faces[0]]; + const auto &second = face_geometry_[record.faces[1]]; + const double cosine = std::clamp(dot(first.normal, second.normal), -1.0, 1.0); + if (cosine <= threshold) { + Quadric penalty; + penalty.add_scaled(first.quadric, feature_weight_); + penalty.add_scaled(second.quadric, feature_weight_); + vertices_[static_cast(edge.first)].quadric += penalty; + vertices_[static_cast(edge.second)].quadric += penalty; + } + } + } + + bool boundary_policy_allows(const Edge &edge) const { + const auto edge_it = edges_.find(edge); + if (edge_it == edges_.end()) return false; + const bool first_boundary = is_boundary_vertex(edge.first); + const bool second_boundary = is_boundary_vertex(edge.second); + if (first_boundary != second_boundary) return false; + if (first_boundary && edge_it->second.count != 1) return false; + return true; + } + + bool link_condition_holds(const Edge &edge) const { + const auto edge_it = edges_.find(edge); + if (edge_it == edges_.end()) return false; + std::unordered_set expected; + for (std::size_t index = 0; index < edge_it->second.count; ++index) { + const auto &triangle = faces_[edge_it->second.faces[index]].vertices; + for (const NodeId vertex : triangle) { + if (vertex != edge.first && vertex != edge.second) expected.insert(vertex); + } + } + std::unordered_set common; + const auto &first_neighbours = vertices_[static_cast(edge.first)].neighbours; + const auto &second_neighbours = vertices_[static_cast(edge.second)].neighbours; + for (const NodeId neighbour : first_neighbours) { + if (neighbour != edge.second && second_neighbours.contains(neighbour)) { + common.insert(neighbour); + } + } + return common == expected; + } + + Candidate best_of_samples(const Quadric &quadric, const Edge &edge) const { + const Vec3 &first = vertices_[static_cast(edge.first)].position; + const Vec3 &second = vertices_[static_cast(edge.second)].position; + const Vec3 midpoint = add_scaled(first, subtract(second, first), 0.5); + const std::array points{first, second, midpoint}; + const std::array parameters{0.0, 1.0, 0.5}; + std::size_t best = 0; + double best_cost = quadric.evaluate(points[0]); + for (std::size_t index = 1; index < points.size(); ++index) { + const double cost = quadric.evaluate(points[index]); + if (cost < best_cost) { + best = index; + best_cost = cost; + } + } + return {points[best], std::max(0.0, best_cost), parameters[best]}; + } + + std::optional make_candidate(const Edge &edge) const { + const auto edge_it = edges_.find(edge); + if (edge_it == edges_.end() || !boundary_policy_allows(edge) || !link_condition_holds(edge)) { + return std::nullopt; + } + const auto &first_vertex = vertices_[static_cast(edge.first)]; + const auto &second_vertex = vertices_[static_cast(edge.second)]; + const Quadric quadric = first_vertex.quadric + second_vertex.quadric; + Candidate candidate; + if (edge_it->second.count == 1) { + const Vec3 direction = subtract(second_vertex.position, first_vertex.position); + const std::array origin{ + first_vertex.position[0], first_vertex.position[1], first_vertex.position[2], 1.0, + }; + const std::array delta{direction[0], direction[1], direction[2], 0.0}; + const double a = quadric.bilinear(delta, delta); + const double b = 2.0 * quadric.bilinear(origin, delta); + if (a > 64.0 * std::numeric_limits::epsilon()) { + candidate.interpolation = std::clamp(-b / (2.0 * a), 0.0, 1.0); + candidate.position = add_scaled( + first_vertex.position, direction, candidate.interpolation + ); + candidate.cost = std::max(0.0, quadric.evaluate(candidate.position)); + } else { + candidate = best_of_samples(quadric, edge); + } + } else if (!solve_3x3(quadric, candidate.position)) { + candidate = best_of_samples(quadric, edge); + } else { + const Vec3 direction = subtract(second_vertex.position, first_vertex.position); + candidate.interpolation = std::clamp( + dot(subtract(candidate.position, first_vertex.position), direction) / dot(direction, direction), + 0.0, + 1.0 + ); + candidate.cost = std::max(0.0, quadric.evaluate(candidate.position)); + } + if (!std::isfinite(candidate.cost)) return std::nullopt; + return candidate; + } + + std::vector affected_faces(const Edge &edge) const { + std::vector result; + const auto &first = vertices_[static_cast(edge.first)].faces; + const auto &second = vertices_[static_cast(edge.second)].faces; + result.reserve(first.size() + second.size()); + result.insert(result.end(), first.begin(), first.end()); + result.insert(result.end(), second.begin(), second.end()); + std::sort(result.begin(), result.end()); + result.erase(std::unique(result.begin(), result.end()), result.end()); + return result; + } + + bool collapse_is_valid(const Edge &edge, const Vec3 &position) const { + if (!boundary_policy_allows(edge) || !link_condition_holds(edge)) return false; + const auto affected = affected_faces(edge); + const std::unordered_set affected_set(affected.begin(), affected.end()); + std::unordered_set new_faces; + for (const FaceId face_id : affected) { + const auto &face = faces_[face_id]; + const bool has_first = std::find( + face.vertices.begin(), face.vertices.end(), edge.first + ) != face.vertices.end(); + const bool has_second = std::find( + face.vertices.begin(), face.vertices.end(), edge.second + ) != face.vertices.end(); + if (has_first && has_second) continue; + + Triangle updated = face.vertices; + for (NodeId &vertex : updated) { + if (vertex == edge.second) vertex = edge.first; + } + Vec3 positions[3]; + for (std::size_t corner = 0; corner < 3; ++corner) { + positions[corner] = updated[corner] == edge.first + ? position + : vertices_[static_cast(updated[corner])].position; + } + const Vec3 new_normal = cross( + subtract(positions[1], positions[0]), subtract(positions[2], positions[0]) + ); + const auto &old_normal = face_geometry_[face_id].normal; + if (norm(new_normal) <= area_tolerance_ || dot(new_normal, old_normal) <= 0.0) { + return false; + } + const Triangle key = triangle_key(updated); + if (!new_faces.insert(key).second) return false; + const auto existing = face_lookup_.find(key); + if (existing != face_lookup_.end() && !affected_set.contains(existing->second)) { + return false; + } + } + return true; + } + + void build_heap() { + heap_.reserve(edges_.size()); + std::vector entries; + entries.reserve(edges_.size()); + for (const auto &[edge, _] : edges_) { + const auto candidate = make_candidate(edge); + if (candidate.has_value() && collapse_is_valid(edge, candidate->position)) { + entries.push_back({edge, {candidate->cost, edge.first, edge.second}}); + } + } + heap_.build_heap(std::move(entries)); + } + + void erase_face_adjacency(const FaceId face_id) { + const auto triangle = faces_[face_id].vertices; + face_lookup_.erase(triangle_key(triangle)); + for (const NodeId vertex : triangle) { + vertices_[static_cast(vertex)].faces.erase(face_id); + } + remove_edge(triangle[0], triangle[1], face_id); + remove_edge(triangle[1], triangle[2], face_id); + remove_edge(triangle[2], triangle[0], face_id); + } + + void add_face_adjacency(const FaceId face_id) { + const auto triangle = faces_[face_id].vertices; + face_lookup_.emplace(triangle_key(triangle), face_id); + for (const NodeId vertex : triangle) { + vertices_[static_cast(vertex)].faces.insert(face_id); + } + add_edge(triangle[0], triangle[1], face_id); + add_edge(triangle[1], triangle[2], face_id); + add_edge(triangle[2], triangle[0], face_id); + } + + void collapse(const Edge &edge, const Candidate &candidate) { + std::unordered_set dirty_vertices; + dirty_vertices.insert(edge.first); + dirty_vertices.insert(edge.second); + for (const NodeId neighbour : vertices_[static_cast(edge.first)].neighbours) { + dirty_vertices.insert(neighbour); + } + for (const NodeId neighbour : vertices_[static_cast(edge.second)].neighbours) { + dirty_vertices.insert(neighbour); + } + for (const NodeId vertex : dirty_vertices) { + if (!vertices_[static_cast(vertex)].alive) continue; + for (const NodeId neighbour : vertices_[static_cast(vertex)].neighbours) { + heap_.erase(bioimage_cpp::detail::edge_key(vertex, neighbour)); + } + } + + const auto affected = affected_faces(edge); + for (const FaceId face_id : affected) erase_face_adjacency(face_id); + + auto &kept = vertices_[static_cast(edge.first)]; + auto &removed = vertices_[static_cast(edge.second)]; + kept.position = candidate.position; + kept.quadric += removed.quadric; + if (has_values_) { + kept.value = (1.0 - candidate.interpolation) * kept.value + + candidate.interpolation * removed.value; + } + removed.alive = false; + removed.faces.clear(); + removed.neighbours.clear(); + + for (const FaceId face_id : affected) { + auto &face = faces_[face_id]; + bool has_first = false; + bool has_second = false; + for (const NodeId vertex : face.vertices) { + has_first = has_first || vertex == edge.first; + has_second = has_second || vertex == edge.second; + } + if (has_first && has_second) { + face.alive = false; + --alive_faces_; + continue; + } + for (NodeId &vertex : face.vertices) { + if (vertex == edge.second) vertex = edge.first; + } + face_geometry_[face_id] = geometry_of(face.vertices); + add_face_adjacency(face_id); + } + + dirty_vertices.erase(edge.second); + for (const NodeId neighbour : kept.neighbours) dirty_vertices.insert(neighbour); + for (const NodeId vertex : dirty_vertices) { + if (!vertices_[static_cast(vertex)].alive) continue; + std::vector neighbours( + vertices_[static_cast(vertex)].neighbours.begin(), + vertices_[static_cast(vertex)].neighbours.end() + ); + std::sort(neighbours.begin(), neighbours.end()); + for (const NodeId neighbour : neighbours) { + const Edge candidate_edge = bioimage_cpp::detail::edge_key(vertex, neighbour); + const auto next = make_candidate(candidate_edge); + if (next.has_value() && collapse_is_valid(candidate_edge, next->position)) { + heap_.push_or_change( + candidate_edge, + {next->cost, candidate_edge.first, candidate_edge.second} + ); + } else { + heap_.erase(candidate_edge); + } + } + } + } + + SimplifyMeshResult finalize() const { + SimplifyMeshResult result; + std::vector compact(vertices_.size(), -1); + std::size_t output_vertices = 0; + for (std::size_t vertex = 0; vertex < vertices_.size(); ++vertex) { + if (vertices_[vertex].alive && !vertices_[vertex].faces.empty()) { + compact[vertex] = static_cast(output_vertices++); + for (const double coordinate : vertices_[vertex].position) { + result.vertices.push_back(static_cast(coordinate)); + } + if (has_values_) { + if (!result.values.has_value()) result.values.emplace(); + result.values->push_back(static_cast(vertices_[vertex].value)); + } + } + } + result.faces.reserve(alive_faces_ * 3); + for (const auto &face : faces_) { + if (!face.alive) continue; + for (const NodeId vertex : face.vertices) { + result.faces.push_back(compact[static_cast(vertex)]); + } + } + + std::vector accumulated_normals(output_vertices, Vec3{0.0, 0.0, 0.0}); + for (std::size_t face = 0; face < result.faces.size() / 3; ++face) { + Vec3 positions[3]; + for (std::size_t corner = 0; corner < 3; ++corner) { + const auto vertex = static_cast(result.faces[face * 3 + corner]); + positions[corner] = { + static_cast(result.vertices[vertex * 3]), + static_cast(result.vertices[vertex * 3 + 1]), + static_cast(result.vertices[vertex * 3 + 2]), + }; + } + const Vec3 normal = cross( + subtract(positions[1], positions[0]), subtract(positions[2], positions[0]) + ); + for (std::size_t corner = 0; corner < 3; ++corner) { + const auto vertex = static_cast(result.faces[face * 3 + corner]); + for (std::size_t axis = 0; axis < 3; ++axis) { + accumulated_normals[vertex][axis] += normal[axis]; + } + } + } + result.normals.resize(output_vertices * 3); + for (std::size_t vertex = 0; vertex < output_vertices; ++vertex) { + const Vec3 &normal = accumulated_normals[vertex]; + const double length = norm(normal); + for (std::size_t axis = 0; axis < 3; ++axis) { + result.normals[vertex * 3 + axis] = static_cast(normal[axis] / length); + } + } + return result; + } +}; + +} // namespace detail::simplification + +// Simplify an oriented, manifold triangle mesh with constrained quadric-error +// edge collapses. Topology and mesh boundaries are preserved; sharp features +// receive a configurable soft penalty. Optional scalar values are interpolated +// along each collapsed edge and vertex normals are recomputed from final faces. +template +SimplifyMeshResult simplify_mesh( + const ConstArrayView &vertices, + const ConstArrayView &faces, + const double reduction, + const ConstArrayView *values = nullptr, + const double feature_angle = 45.0, + const double feature_weight = 10.0 +) { + if (!std::isfinite(reduction) || reduction < 0.0 || reduction >= 1.0) { + throw std::invalid_argument("reduction must be finite and in [0, 1)"); + } + if (!std::isfinite(feature_angle) || feature_angle < 0.0 || feature_angle > 180.0) { + throw std::invalid_argument("feature_angle must be finite and in [0, 180]"); + } + if (!std::isfinite(feature_weight) || feature_weight < 0.0) { + throw std::invalid_argument("feature_weight must be finite and non-negative"); + } + if (faces.shape.size() != 2) { + throw std::invalid_argument("faces must have shape (n_faces, 3)"); + } + const auto n_faces = static_cast(faces.shape[0]); + const auto target = static_cast( + std::ceil((1.0 - reduction) * static_cast(n_faces)) + ); + detail::simplification::Simplifier simplifier( + vertices, faces, values, feature_angle, feature_weight + ); + return simplifier.run(target); +} + +} // namespace bioimage_cpp::mesh diff --git a/include/bioimage_cpp/mesh_smoothing.hxx b/include/bioimage_cpp/mesh/smoothing.hxx similarity index 95% rename from include/bioimage_cpp/mesh_smoothing.hxx rename to include/bioimage_cpp/mesh/smoothing.hxx index 9637e5a..eb53186 100644 --- a/include/bioimage_cpp/mesh_smoothing.hxx +++ b/include/bioimage_cpp/mesh/smoothing.hxx @@ -12,9 +12,9 @@ #include #include -namespace bioimage_cpp { +namespace bioimage_cpp::mesh { -namespace detail::mesh_smoothing { +namespace detail::smoothing { template struct Adjacency { @@ -74,7 +74,7 @@ Adjacency build_adjacency(const ConstArrayView &faces, std::ptrdiff_t n_ve return Adjacency{std::move(offsets), std::move(neighbours)}; } -} // namespace detail::mesh_smoothing +} // namespace detail::smoothing // Laplacian smoothing of a triangular mesh: each vertex (and corresponding // normal) is replaced by the mean of itself and its 1-ring neighbours, @@ -120,7 +120,7 @@ void smooth_mesh( return; } - const auto adjacency = detail::mesh_smoothing::build_adjacency(faces, n_verts); + const auto adjacency = detail::smoothing::build_adjacency(faces, n_verts); std::vector scratch_verts(n_total); std::vector scratch_normals(n_total); @@ -140,13 +140,13 @@ void smooth_mesh( const auto &offsets = adjacency.offsets; const auto &neighbours = adjacency.neighbours; - const std::size_t threads = detail::normalize_thread_count( + const std::size_t threads = bioimage_cpp::detail::normalize_thread_count( n_threads < 0 ? 0 : static_cast(n_threads), static_cast(n_verts) ); auto smooth_pass = [&](const V *src_verts, const V *src_normals, V *dst_verts, V *dst_normals) { - detail::parallel_for_chunks( + bioimage_cpp::detail::parallel_for_chunks( threads, static_cast(n_verts), [&](std::size_t, std::size_t begin, std::size_t end) { @@ -184,4 +184,4 @@ void smooth_mesh( } } -} // namespace bioimage_cpp +} // namespace bioimage_cpp::mesh diff --git a/src/bindings/mesh.cxx b/src/bindings/mesh.cxx index 7511ebb..4833847 100644 --- a/src/bindings/mesh.cxx +++ b/src/bindings/mesh.cxx @@ -3,9 +3,12 @@ #include "bioimage_cpp/array_view.hxx" #include "bioimage_cpp/detail/profile.hxx" #include "bioimage_cpp/mesh/marching_cubes.hxx" +#include "bioimage_cpp/mesh/simplification.hxx" +#include "bioimage_cpp/mesh/smoothing.hxx" #include #include +#include #include #include @@ -25,6 +28,13 @@ using FloatInput = nb::ndarray; using UInt8Input = nb::ndarray; using FloatOutput = nb::ndarray; using Int32Output = nb::ndarray; +using Int64Output = nb::ndarray; + +template +using InputArray = nb::ndarray; + +template +using OutputArray = nb::ndarray; template Output output_array(std::vector &&values, const std::vector &shape) { @@ -58,6 +68,154 @@ std::vector shape_of(const FloatInput &array) { return shape; } +using FacesArray = nb::ndarray; + +template +std::pair, OutputArray> smooth_mesh_t( + InputArray verts, + InputArray normals, + FacesArray faces, + std::size_t iterations, + int n_threads +) { + if (verts.ndim() != 2) { + throw std::invalid_argument( + "verts must have ndim=2, got ndim=" + std::to_string(verts.ndim()) + ); + } + if (normals.ndim() != 2 || normals.shape(0) != verts.shape(0) + || normals.shape(1) != verts.shape(1)) { + throw std::invalid_argument("normals must have the same shape as verts"); + } + if (faces.ndim() != 2 || faces.shape(1) != 3) { + throw std::invalid_argument("faces must have shape (n_faces, 3)"); + } + + const std::size_t n_verts_size = verts.shape(0); + const std::size_t dim_size = verts.shape(1); + const std::size_t n_faces_size = faces.shape(0); + const std::size_t n_total = n_verts_size * dim_size; + + std::vector verts_ndarray_shape{n_verts_size, dim_size}; + std::vector verts_view_shape{ + static_cast(n_verts_size), + static_cast(dim_size), + }; + std::vector faces_view_shape{ + static_cast(n_faces_size), + std::ptrdiff_t{3}, + }; + + auto *verts_data = new V[n_total](); + nb::capsule verts_owner(verts_data, [](void *p) noexcept { delete[] static_cast(p); }); + auto *normals_data = new V[n_total](); + nb::capsule normals_owner(normals_data, [](void *p) noexcept { delete[] static_cast(p); }); + + ConstArrayView verts_view{verts.data(), verts_view_shape, {}}; + ConstArrayView normals_view{normals.data(), verts_view_shape, {}}; + ConstArrayView faces_view{faces.data(), faces_view_shape, {}}; + ArrayView out_verts_view{verts_data, verts_view_shape, {}}; + ArrayView out_normals_view{normals_data, verts_view_shape, {}}; + + { + nb::gil_scoped_release release; + mesh::smooth_mesh( + verts_view, + normals_view, + faces_view, + iterations, + n_threads, + out_verts_view, + out_normals_view + ); + } + + OutputArray out_verts( + verts_data, verts_ndarray_shape.size(), verts_ndarray_shape.data(), verts_owner + ); + OutputArray out_normals( + normals_data, verts_ndarray_shape.size(), verts_ndarray_shape.data(), normals_owner + ); + return {std::move(out_verts), std::move(out_normals)}; +} + +template +nb::tuple simplify_mesh_t( + InputArray vertices, + FacesArray faces, + double reduction, + std::optional> values, + double feature_angle, + double feature_weight +) { + if (vertices.ndim() != 2 || vertices.shape(1) != 3) { + throw std::invalid_argument("vertices must have shape (n_vertices, 3)"); + } + if (faces.ndim() != 2 || faces.shape(1) != 3) { + throw std::invalid_argument("faces must have shape (n_faces, 3)"); + } + if (values.has_value() + && (values->ndim() != 1 || values->shape(0) != vertices.shape(0))) { + throw std::invalid_argument("values must have shape (n_vertices,)"); + } + + const std::vector vertex_shape{ + static_cast(vertices.shape(0)), std::ptrdiff_t{3} + }; + const std::vector face_shape{ + static_cast(faces.shape(0)), std::ptrdiff_t{3} + }; + ConstArrayView vertex_view{vertices.data(), vertex_shape, {}}; + ConstArrayView face_view{faces.data(), face_shape, {}}; + ConstArrayView value_view; + const ConstArrayView *value_ptr = nullptr; + if (values.has_value()) { + value_view = ConstArrayView{ + values->data(), {static_cast(values->shape(0))}, {} + }; + value_ptr = &value_view; + } + + mesh::SimplifyMeshResult result; + { + nb::gil_scoped_release release; + result = mesh::simplify_mesh( + vertex_view, + face_view, + reduction, + value_ptr, + feature_angle, + feature_weight + ); + } + + const std::size_t n_vertices = result.vertices.size() / 3; + const std::size_t n_faces = result.faces.size() / 3; + auto out_vertices = output_array>( + std::move(result.vertices), {n_vertices, 3} + ); + auto out_faces = output_array( + std::move(result.faces), {n_faces, 3} + ); + auto out_normals = output_array>( + std::move(result.normals), {n_vertices, 3} + ); + if (!result.values.has_value()) { + return nb::make_tuple( + std::move(out_vertices), std::move(out_faces), std::move(out_normals), nb::none() + ); + } + auto out_values = output_array>( + std::move(*result.values), {n_vertices} + ); + return nb::make_tuple( + std::move(out_vertices), + std::move(out_faces), + std::move(out_normals), + std::move(out_values) + ); +} + nb::tuple marching_cubes_float32( FloatInput volume, const double level, @@ -150,6 +308,70 @@ nb::tuple marching_cubes_float32( } // namespace void bind_mesh(nb::module_ &m) { + m.def( + "_simplify_mesh_float32_float32", + &simplify_mesh_t, + nb::arg("vertices"), + nb::arg("faces"), + nb::arg("reduction"), + nb::arg("values") = nb::none(), + nb::arg("feature_angle") = 45.0, + nb::arg("feature_weight") = 10.0, + "Constrained QEM simplification with float32 vertices and values." + ); + m.def( + "_simplify_mesh_float32_float64", + &simplify_mesh_t, + nb::arg("vertices"), + nb::arg("faces"), + nb::arg("reduction"), + nb::arg("values") = nb::none(), + nb::arg("feature_angle") = 45.0, + nb::arg("feature_weight") = 10.0, + "Constrained QEM simplification with float32 vertices and float64 values." + ); + m.def( + "_simplify_mesh_float64_float32", + &simplify_mesh_t, + nb::arg("vertices"), + nb::arg("faces"), + nb::arg("reduction"), + nb::arg("values") = nb::none(), + nb::arg("feature_angle") = 45.0, + nb::arg("feature_weight") = 10.0, + "Constrained QEM simplification with float64 vertices and float32 values." + ); + m.def( + "_simplify_mesh_float64_float64", + &simplify_mesh_t, + nb::arg("vertices"), + nb::arg("faces"), + nb::arg("reduction"), + nb::arg("values") = nb::none(), + nb::arg("feature_angle") = 45.0, + nb::arg("feature_weight") = 10.0, + "Constrained QEM simplification with float64 vertices and values." + ); + m.def( + "_smooth_mesh_float32", + &smooth_mesh_t, + nb::arg("verts"), + nb::arg("normals"), + nb::arg("faces"), + nb::arg("iterations"), + nb::arg("n_threads"), + "Laplacian smoothing of a triangular mesh with float32 vertices and normals." + ); + m.def( + "_smooth_mesh_float64", + &smooth_mesh_t, + nb::arg("verts"), + nb::arg("normals"), + nb::arg("faces"), + nb::arg("iterations"), + nb::arg("n_threads"), + "Laplacian smoothing of a triangular mesh with float64 vertices and normals." + ); m.def( "_marching_cubes_float32", &marching_cubes_float32, diff --git a/src/bindings/utils.cxx b/src/bindings/utils.cxx index 38a84dd..f2be2c4 100644 --- a/src/bindings/utils.cxx +++ b/src/bindings/utils.cxx @@ -1,12 +1,10 @@ #include "utils.hxx" #include "bioimage_cpp/array_view.hxx" -#include "bioimage_cpp/mesh_smoothing.hxx" #include "bioimage_cpp/run_length.hxx" #include "bioimage_cpp/take_dict.hxx" #include -#include #include #include @@ -82,76 +80,6 @@ OutputArray take_dict_t(const nb::dict &relabeling, InputArray to_relabel) return OutputArray(data, ndarray_shape.size(), ndarray_shape.data(), owner); } -using FacesArray = nb::ndarray; - -template -std::pair, OutputArray> smooth_mesh_t( - InputArray verts, - InputArray normals, - FacesArray faces, - std::size_t iterations, - int n_threads -) { - if (verts.ndim() != 2) { - throw std::invalid_argument( - "verts must have ndim=2, got ndim=" + std::to_string(verts.ndim()) - ); - } - if (normals.ndim() != 2 || normals.shape(0) != verts.shape(0) - || normals.shape(1) != verts.shape(1)) { - throw std::invalid_argument("normals must have the same shape as verts"); - } - if (faces.ndim() != 2 || faces.shape(1) != 3) { - throw std::invalid_argument("faces must have shape (n_faces, 3)"); - } - - const std::size_t n_verts_size = verts.shape(0); - const std::size_t dim_size = verts.shape(1); - const std::size_t n_faces_size = faces.shape(0); - const std::size_t n_total = n_verts_size * dim_size; - - std::vector verts_ndarray_shape{n_verts_size, dim_size}; - std::vector verts_view_shape{ - static_cast(n_verts_size), - static_cast(dim_size), - }; - std::vector faces_view_shape{ - static_cast(n_faces_size), - std::ptrdiff_t{3}, - }; - - auto *verts_data = new V[n_total](); - nb::capsule verts_owner(verts_data, [](void *p) noexcept { delete[] static_cast(p); }); - auto *normals_data = new V[n_total](); - nb::capsule normals_owner(normals_data, [](void *p) noexcept { delete[] static_cast(p); }); - - ConstArrayView verts_view{verts.data(), verts_view_shape, {}}; - ConstArrayView normals_view{normals.data(), verts_view_shape, {}}; - ConstArrayView faces_view{faces.data(), faces_view_shape, {}}; - ArrayView out_verts_view{verts_data, verts_view_shape, {}}; - ArrayView out_normals_view{normals_data, verts_view_shape, {}}; - - { - nb::gil_scoped_release release; - smooth_mesh( - verts_view, - normals_view, - faces_view, - iterations, - n_threads, - out_verts_view, - out_normals_view - ); - } - - OutputArray out_verts( - verts_data, verts_ndarray_shape.size(), verts_ndarray_shape.data(), verts_owner - ); - OutputArray out_normals( - normals_data, verts_ndarray_shape.size(), verts_ndarray_shape.data(), normals_owner - ); - return {std::move(out_verts), std::move(out_normals)}; -} template OutputArray compute_rle_t(InputArray mask) { @@ -213,26 +141,6 @@ void bind_utils(nb::module_ &m) { nb::arg("to_relabel"), "Map a contiguous int64 array through an integer dictionary." ); - m.def( - "_smooth_mesh_float32", - &smooth_mesh_t, - nb::arg("verts"), - nb::arg("normals"), - nb::arg("faces"), - nb::arg("iterations"), - nb::arg("n_threads"), - "Laplacian smoothing of a triangular mesh with float32 vertices and normals." - ); - m.def( - "_smooth_mesh_float64", - &smooth_mesh_t, - nb::arg("verts"), - nb::arg("normals"), - nb::arg("faces"), - nb::arg("iterations"), - nb::arg("n_threads"), - "Laplacian smoothing of a triangular mesh with float64 vertices and normals." - ); m.def("_compute_rle_bool", &compute_rle_t, nb::arg("mask"), "COCO-style binary run-length encoding of a contiguous bool array."); m.def("_compute_rle_uint8", &compute_rle_t, nb::arg("mask"), diff --git a/src/bioimage_cpp/__init__.py b/src/bioimage_cpp/__init__.py index e5ce86f..b466c50 100644 --- a/src/bioimage_cpp/__init__.py +++ b/src/bioimage_cpp/__init__.py @@ -56,7 +56,7 @@ - `distance`: distance transform functionality. - `filters`: efficient implementation of convolutional image filters. - `graph`: graph creation and graph (partitioning) algorithms. -- `mesh`: triangle-mesh extraction from 3D volumes and segmentation masks. +- `mesh`: triangle-mesh extraction and processing. - `segmentation`: image segmentation functionality. - `transformation`: affine transformations. - `utils`: misc utility functionality. diff --git a/src/bioimage_cpp/mesh/__init__.py b/src/bioimage_cpp/mesh/__init__.py index edb9dc3..57a3627 100644 --- a/src/bioimage_cpp/mesh/__init__.py +++ b/src/bioimage_cpp/mesh/__init__.py @@ -1,4 +1,4 @@ -"""Triangle-mesh extraction from 3-D scalar volumes.""" +"""Triangle-mesh extraction and processing.""" from __future__ import annotations @@ -9,6 +9,18 @@ from .. import _core +_SMOOTH_MESH_BY_DTYPE = { + np.dtype("float32"): _core._smooth_mesh_float32, + np.dtype("float64"): _core._smooth_mesh_float64, +} + +_SIMPLIFY_MESH_BY_DTYPE = { + (np.dtype("float32"), np.dtype("float32")): _core._simplify_mesh_float32_float32, + (np.dtype("float32"), np.dtype("float64")): _core._simplify_mesh_float32_float64, + (np.dtype("float64"), np.dtype("float32")): _core._simplify_mesh_float64_float32, + (np.dtype("float64"), np.dtype("float64")): _core._simplify_mesh_float64_float64, +} + def _as_spacing(spacing: float | Sequence[float]) -> np.ndarray: if np.isscalar(spacing): @@ -176,4 +188,248 @@ def marching_cubes( return vertices, faces, normals, values -__all__ = ["marching_cubes"] +def smooth_mesh( + verts: np.ndarray, + normals: np.ndarray, + faces: np.ndarray, + iterations: int, + *, + n_threads: int = 0, +) -> tuple[np.ndarray, np.ndarray]: + """Laplacian smoothing of a triangular mesh. + + Each vertex (and the corresponding normal) is iteratively replaced by the + mean of itself and its 1-ring neighbours in the mesh. Updates are applied + in Jacobi (double-buffered) style: every vertex in an iteration is updated + from the same input state, independent of vertex order. The adjacency is + built once from ``faces``; each interior edge is treated as a single + undirected edge regardless of how many faces share it. + + Parameters + ---------- + verts: + 2D ``float32`` or ``float64`` array of shape ``(n_verts, dim)``. + normals: + Array with the same shape and dtype as ``verts``. + faces: + 2D integer array of shape ``(n_faces, 3)`` with values in + ``[0, n_verts)``. Non-``int64`` inputs are converted internally. + iterations: + Number of smoothing passes. ``0`` returns copies of the inputs. + n_threads: + Number of threads for the inner smoothing loop. ``0`` (the default) + picks the hardware concurrency; ``1`` forces serial execution. + + Returns + ------- + (np.ndarray, np.ndarray) + Smoothed vertices and normals with the same shape and dtype as the + inputs. + """ + verts_array = np.asarray(verts) + normals_array = np.asarray(normals) + faces_array = np.asarray(faces) + + if verts_array.ndim != 2: + raise ValueError( + f"verts must have ndim=2, got ndim={verts_array.ndim}" + ) + if normals_array.shape != verts_array.shape: + raise ValueError( + "normals must have the same shape as verts, got " + f"verts shape={verts_array.shape}, normals shape={normals_array.shape}" + ) + if verts_array.dtype != normals_array.dtype: + raise TypeError( + "verts and normals must have the same dtype, got " + f"verts dtype={verts_array.dtype}, normals dtype={normals_array.dtype}" + ) + + try: + smooth = _SMOOTH_MESH_BY_DTYPE[verts_array.dtype] + except KeyError as error: + supported = ", ".join(str(dtype) for dtype in _SMOOTH_MESH_BY_DTYPE) + raise TypeError( + f"verts must have one of dtypes ({supported}), got dtype={verts_array.dtype}" + ) from error + + if faces_array.ndim != 2 or faces_array.shape[1] != 3: + raise ValueError( + f"faces must have shape (n_faces, 3), got shape={faces_array.shape}" + ) + if not np.issubdtype(faces_array.dtype, np.integer): + raise TypeError( + f"faces must have an integer dtype, got dtype={faces_array.dtype}" + ) + + n_verts = verts_array.shape[0] + if faces_array.size > 0: + face_min = int(faces_array.min()) + face_max = int(faces_array.max()) + if face_min < 0 or face_max >= n_verts: + raise ValueError( + "faces must contain indices in [0, n_verts), got values in " + f"[{face_min}, {face_max}] with n_verts={n_verts}" + ) + + iterations_int = int(iterations) + if iterations_int < 0: + raise ValueError(f"iterations must be non-negative, got {iterations_int}") + + n_threads_int = int(n_threads) + if n_threads_int < 0: + raise ValueError(f"n_threads must be non-negative, got {n_threads_int}") + + verts_contiguous = np.ascontiguousarray(verts_array) + normals_contiguous = np.ascontiguousarray(normals_array) + faces_contiguous = np.ascontiguousarray(faces_array, dtype=np.int64) + + return smooth( + verts_contiguous, + normals_contiguous, + faces_contiguous, + iterations_int, + n_threads_int, + ) + + +def simplify_mesh( + vertices: np.ndarray, + faces: np.ndarray, + reduction: float, + *, + values: np.ndarray | None = None, + feature_angle: float = 45.0, + feature_weight: float = 10.0, +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray | None]: + """Simplify an oriented manifold triangle mesh with QEM edge collapses. + + The requested reduction is the approximate fraction of input faces to + remove. Edge collapses preserve topology and mesh boundaries. Interior + edges whose dihedral angle is at least feature_angle receive a soft + quadric penalty scaled by feature_weight. A constrained mesh may stop + above the requested target when no safe collapse remains. + + Parameters + ---------- + vertices: + Float32 or float64 array of shape (n_vertices, 3). + faces: + Integer array of shape (n_faces, 3) describing a consistently + oriented triangle 2-manifold. Degenerate and non-manifold inputs are + rejected. + reduction: + Finite face-reduction fraction in [0, 1). + values: + Optional float32 or float64 scalar array of shape (n_vertices,). + Values are linearly interpolated along collapsed edges and preserve + their input dtype. + feature_angle: + Sharp-feature threshold in degrees, in [0, 180]. + feature_weight: + Non-negative soft-feature penalty. Zero disables feature penalties. + + Returns + ------- + vertices, faces, normals, values: + Compact simplified vertices, int64 triangle indices, recomputed + area-weighted unit vertex normals, and interpolated values. The final + item is None when values was not supplied. + + Notes + ----- + For meshes from marching_cubes, use allow_degenerate=False so the + simplifier receives valid triangles. + """ + vertices_array = np.asarray(vertices) + faces_array = np.asarray(faces) + if vertices_array.ndim != 2 or vertices_array.shape[1] != 3: + raise ValueError( + "vertices must have shape (n_vertices, 3), got " + f"shape={vertices_array.shape}" + ) + if vertices_array.dtype not in (np.dtype("float32"), np.dtype("float64")): + raise TypeError( + "vertices must have dtype float32 or float64, got " + f"dtype={vertices_array.dtype}" + ) + if faces_array.ndim != 2 or faces_array.shape[1] != 3: + raise ValueError( + f"faces must have shape (n_faces, 3), got shape={faces_array.shape}" + ) + if not np.issubdtype(faces_array.dtype, np.integer): + raise TypeError(f"faces must have an integer dtype, got dtype={faces_array.dtype}") + if len(vertices_array) == 0 or len(faces_array) == 0: + raise ValueError("vertices and faces must describe a non-empty mesh") + if not np.all(np.isfinite(vertices_array)): + raise ValueError("vertices must contain only finite values") + face_min = int(faces_array.min()) + face_max = int(faces_array.max()) + if face_min < 0 or face_max >= len(vertices_array): + raise ValueError( + "faces must contain indices in [0, n_vertices), got values in " + f"[{face_min}, {face_max}] with n_vertices={len(vertices_array)}" + ) + + try: + reduction_float = float(reduction) + except (TypeError, ValueError) as error: + raise TypeError("reduction must be a real number") from error + if not np.isfinite(reduction_float) or not 0.0 <= reduction_float < 1.0: + raise ValueError( + f"reduction must be finite and in [0, 1), got {reduction_float}" + ) + try: + feature_angle_float = float(feature_angle) + except (TypeError, ValueError) as error: + raise TypeError("feature_angle must be a real number") from error + if ( + not np.isfinite(feature_angle_float) + or feature_angle_float < 0.0 + or feature_angle_float > 180.0 + ): + raise ValueError( + "feature_angle must be finite and in [0, 180], got " + f"{feature_angle_float}" + ) + try: + feature_weight_float = float(feature_weight) + except (TypeError, ValueError) as error: + raise TypeError("feature_weight must be a real number") from error + if not np.isfinite(feature_weight_float) or feature_weight_float < 0.0: + raise ValueError( + "feature_weight must be finite and non-negative, got " + f"{feature_weight_float}" + ) + + values_array: np.ndarray | None = None + value_dtype = vertices_array.dtype + if values is not None: + values_array = np.asarray(values) + if values_array.shape != (len(vertices_array),): + raise ValueError( + "values must have shape (n_vertices,), got " + f"shape={values_array.shape}" + ) + if values_array.dtype not in (np.dtype("float32"), np.dtype("float64")): + raise TypeError( + "values must have dtype float32 or float64, got " + f"dtype={values_array.dtype}" + ) + if not np.all(np.isfinite(values_array)): + raise ValueError("values must contain only finite values") + value_dtype = values_array.dtype + values_array = np.ascontiguousarray(values_array) + + simplify = _SIMPLIFY_MESH_BY_DTYPE[(vertices_array.dtype, value_dtype)] + return simplify( + np.ascontiguousarray(vertices_array), + np.ascontiguousarray(faces_array, dtype=np.int64), + reduction_float, + values_array, + feature_angle_float, + feature_weight_float, + ) + + +__all__ = ["marching_cubes", "simplify_mesh", "smooth_mesh"] diff --git a/src/bioimage_cpp/utils.py b/src/bioimage_cpp/utils.py index ff99fd1..52c6402 100644 --- a/src/bioimage_cpp/utils.py +++ b/src/bioimage_cpp/utils.py @@ -19,11 +19,6 @@ np.dtype("int64"): _core._take_dict_int64, } -_SMOOTH_MESH_BY_DTYPE = { - np.dtype("float32"): _core._smooth_mesh_float32, - np.dtype("float64"): _core._smooth_mesh_float64, -} - _COMPUTE_RLE_BY_DTYPE = { np.dtype("bool"): _core._compute_rle_bool, np.dtype("uint8"): _core._compute_rle_uint8, @@ -358,110 +353,6 @@ def compute_rle(mask: np.ndarray) -> np.ndarray: return run(np.ascontiguousarray(array)) -def smooth_mesh( - verts: np.ndarray, - normals: np.ndarray, - faces: np.ndarray, - iterations: int, - *, - n_threads: int = 0, -) -> tuple[np.ndarray, np.ndarray]: - """Laplacian smoothing of a triangular mesh. - - Each vertex (and the corresponding normal) is iteratively replaced by the - mean of itself and its 1-ring neighbours in the mesh. Updates are applied - in Jacobi (double-buffered) style: every vertex in an iteration is updated - from the same input state, independent of vertex order. The adjacency is - built once from ``faces``; each interior edge is treated as a single - undirected edge regardless of how many faces share it. - - Parameters - ---------- - verts: - 2D ``float32`` or ``float64`` array of shape ``(n_verts, dim)``. - normals: - Array with the same shape and dtype as ``verts``. - faces: - 2D integer array of shape ``(n_faces, 3)`` with values in - ``[0, n_verts)``. Non-``int64`` inputs are converted internally. - iterations: - Number of smoothing passes. ``0`` returns copies of the inputs. - n_threads: - Number of threads for the inner smoothing loop. ``0`` (the default) - picks the hardware concurrency; ``1`` forces serial execution. - - Returns - ------- - (np.ndarray, np.ndarray) - Smoothed vertices and normals with the same shape and dtype as the - inputs. - """ - verts_array = np.asarray(verts) - normals_array = np.asarray(normals) - faces_array = np.asarray(faces) - - if verts_array.ndim != 2: - raise ValueError( - f"verts must have ndim=2, got ndim={verts_array.ndim}" - ) - if normals_array.shape != verts_array.shape: - raise ValueError( - "normals must have the same shape as verts, got " - f"verts shape={verts_array.shape}, normals shape={normals_array.shape}" - ) - if verts_array.dtype != normals_array.dtype: - raise TypeError( - "verts and normals must have the same dtype, got " - f"verts dtype={verts_array.dtype}, normals dtype={normals_array.dtype}" - ) - - try: - smooth = _SMOOTH_MESH_BY_DTYPE[verts_array.dtype] - except KeyError as error: - supported = ", ".join(str(dtype) for dtype in _SMOOTH_MESH_BY_DTYPE) - raise TypeError( - f"verts must have one of dtypes ({supported}), got dtype={verts_array.dtype}" - ) from error - - if faces_array.ndim != 2 or faces_array.shape[1] != 3: - raise ValueError( - f"faces must have shape (n_faces, 3), got shape={faces_array.shape}" - ) - if not np.issubdtype(faces_array.dtype, np.integer): - raise TypeError( - f"faces must have an integer dtype, got dtype={faces_array.dtype}" - ) - - n_verts = verts_array.shape[0] - if faces_array.size > 0: - face_min = int(faces_array.min()) - face_max = int(faces_array.max()) - if face_min < 0 or face_max >= n_verts: - raise ValueError( - "faces must contain indices in [0, n_verts), got values in " - f"[{face_min}, {face_max}] with n_verts={n_verts}" - ) - - iterations_int = int(iterations) - if iterations_int < 0: - raise ValueError(f"iterations must be non-negative, got {iterations_int}") - - n_threads_int = int(n_threads) - if n_threads_int < 0: - raise ValueError(f"n_threads must be non-negative, got {n_threads_int}") - - verts_contiguous = np.ascontiguousarray(verts_array) - normals_contiguous = np.ascontiguousarray(normals_array) - faces_contiguous = np.ascontiguousarray(faces_array, dtype=np.int64) - - return smooth( - verts_contiguous, - normals_contiguous, - faces_contiguous, - iterations_int, - n_threads_int, - ) - def _normalize_labels(labels: np.ndarray, name: str) -> np.ndarray: array = np.asarray(labels) @@ -526,6 +417,5 @@ def _validate_normalize_by(normalize_by: str) -> None: "UnionFind", "compute_rle", "segmentation_overlap", - "smooth_mesh", "take_dict", ] diff --git a/tests/mesh/test_simplification.py b/tests/mesh/test_simplification.py new file mode 100644 index 0000000..04efd79 --- /dev/null +++ b/tests/mesh/test_simplification.py @@ -0,0 +1,247 @@ +import numpy as np +import pytest + +import bioimage_cpp as bic + + +def _sphere_mesh(dtype=np.float32): + n = 22 + center = (n - 1) / 2.0 + z, y, x = np.ogrid[:n, :n, :n] + volume = ( + (z - center) ** 2 + (y - center) ** 2 + (x - center) ** 2 <= 7.0**2 + ).astype(np.uint8) + vertices, faces, _, values = bic.mesh.marching_cubes( + volume, 0.5, allow_degenerate=False + ) + return vertices.astype(dtype), faces, values + + +def _edge_counts(faces): + edges = np.concatenate( + [faces[:, [0, 1]], faces[:, [1, 2]], faces[:, [2, 0]]], axis=0 + ) + edges.sort(axis=1) + _, counts = np.unique(edges, axis=0, return_counts=True) + return counts + + +def _euler_characteristic(vertices, faces): + edges = np.concatenate( + [faces[:, [0, 1]], faces[:, [1, 2]], faces[:, [2, 0]]], axis=0 + ) + edges.sort(axis=1) + return len(vertices) - len(np.unique(edges, axis=0)) + len(faces) + + +def _grid_mesh(size=6, dtype=np.float64): + y, x = np.mgrid[:size, :size] + vertices = np.stack([np.zeros_like(x), y, x], axis=-1).reshape(-1, 3).astype(dtype) + faces = [] + for row in range(size - 1): + for col in range(size - 1): + first = row * size + col + faces.append((first, first + size, first + 1)) + faces.append((first + 1, first + size, first + size + 1)) + return vertices, np.asarray(faces, dtype=np.int64) + + +def _tetrahedron(dtype=np.float64): + vertices = np.asarray( + [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]], dtype=dtype + ) + faces = np.asarray( + [[0, 2, 1], [0, 1, 3], [0, 3, 2], [1, 2, 3]], dtype=np.int64 + ) + return vertices, faces + + +def _reference_normals(vertices, faces): + triangles = vertices[faces].astype(np.float64) + face_normals = np.cross( + triangles[:, 1] - triangles[:, 0], + triangles[:, 2] - triangles[:, 0], + ) + normals = np.zeros(vertices.shape, dtype=np.float64) + for corner in range(3): + np.add.at(normals, faces[:, corner], face_normals) + normals /= np.linalg.norm(normals, axis=1, keepdims=True) + return normals + + +def test_simplifies_closed_marching_cubes_mesh_and_preserves_topology(): + vertices, faces, values = _sphere_mesh() + target = int(np.ceil(0.5 * len(faces))) + + out_vertices, out_faces, out_normals, out_values = bic.mesh.simplify_mesh( + vertices, faces, 0.5, values=values + ) + + assert len(out_faces) <= target + assert len(out_vertices) < len(vertices) + assert out_vertices.dtype == out_normals.dtype == np.float32 + assert out_faces.dtype == np.int64 + assert out_values.dtype == values.dtype + assert out_values.shape == (len(out_vertices),) + assert set(_edge_counts(out_faces)) == {2} + assert _euler_characteristic(out_vertices, out_faces) == 2 + np.testing.assert_allclose( + np.linalg.norm(out_normals, axis=1), 1.0, rtol=2e-6, atol=2e-6 + ) + np.testing.assert_allclose( + out_normals, _reference_normals(out_vertices, out_faces), rtol=2e-6, atol=2e-6 + ) + assert out_values.min() >= values.min() + assert out_values.max() <= values.max() + + +def test_open_boundary_remains_a_single_manifold_loop(): + vertices, faces = _grid_mesh() + out_vertices, out_faces, _, out_values = bic.mesh.simplify_mesh( + vertices, faces, 0.4 + ) + + counts = _edge_counts(out_faces) + assert set(counts) == {1, 2} + assert _euler_characteristic(out_vertices, out_faces) == 1 + assert out_values is None + + edges = np.concatenate( + [out_faces[:, [0, 1]], out_faces[:, [1, 2]], out_faces[:, [2, 0]]], + axis=0, + ) + edges.sort(axis=1) + unique_edges, edge_counts = np.unique(edges, axis=0, return_counts=True) + boundary = unique_edges[edge_counts == 1] + degrees = np.bincount(boundary.ravel(), minlength=len(out_vertices)) + np.testing.assert_array_equal(degrees[degrees > 0], 2) + + +def test_tetrahedron_stops_before_topology_would_change(): + vertices, faces = _tetrahedron() + output = bic.mesh.simplify_mesh(vertices, faces, 0.9) + np.testing.assert_array_equal(output[0], vertices) + np.testing.assert_array_equal(output[1], faces) + + +@pytest.mark.parametrize("vertex_dtype", [np.float32, np.float64]) +@pytest.mark.parametrize("value_dtype", [np.float32, np.float64]) +def test_vertex_and_value_dtypes_are_preserved(vertex_dtype, value_dtype): + vertices, faces = _grid_mesh(dtype=vertex_dtype) + values = np.linspace(0.0, 1.0, len(vertices), dtype=value_dtype) + out_vertices, out_faces, out_normals, out_values = bic.mesh.simplify_mesh( + vertices[::1], faces.astype(np.int32), 0.2, values=values + ) + assert out_vertices.dtype == out_normals.dtype == vertex_dtype + assert out_faces.dtype == np.int64 + assert out_values.dtype == value_dtype + + +def test_noncontiguous_inputs_are_normalized_at_python_boundary(): + vertices, faces = _grid_mesh() + values = np.linspace(0.0, 1.0, len(vertices)) + large_vertices = np.repeat(vertices, 2, axis=0) + large_values = np.repeat(values, 2) + large_faces = np.empty((len(faces), 6), dtype=np.int64) + large_faces[:, ::2] = faces + vertex_view = large_vertices[::2] + face_view = large_faces[:, ::2] + value_view = large_values[::2] + assert not vertex_view.flags.c_contiguous + assert not face_view.flags.c_contiguous + assert not value_view.flags.c_contiguous + + expected = bic.mesh.simplify_mesh(vertices, faces, 0.3, values=values) + actual = bic.mesh.simplify_mesh( + vertex_view, face_view, 0.3, values=value_view + ) + for got, want in zip(actual, expected, strict=True): + np.testing.assert_array_equal(got, want) + + +def test_reduction_zero_is_a_geometry_copy_and_recomputes_normals(): + vertices, faces = _tetrahedron(np.float32) + out_vertices, out_faces, out_normals, out_values = bic.mesh.simplify_mesh( + vertices, faces, 0.0 + ) + np.testing.assert_array_equal(out_vertices, vertices) + np.testing.assert_array_equal(out_faces, faces) + np.testing.assert_allclose( + out_normals, _reference_normals(vertices, faces), rtol=2e-6, atol=2e-6 + ) + assert out_values is None + + +def test_output_is_deterministic(): + vertices, faces, values = _sphere_mesh(np.float64) + first = bic.mesh.simplify_mesh( + vertices, faces, 0.35, values=values, feature_angle=30.0, feature_weight=20.0 + ) + second = bic.mesh.simplify_mesh( + vertices, faces, 0.35, values=values, feature_angle=30.0, feature_weight=20.0 + ) + for actual, expected in zip(first, second, strict=True): + np.testing.assert_array_equal(actual, expected) + + +def test_soft_feature_weight_improves_fold_preservation(): + size = 9 + y, x = np.mgrid[:size, :size] + vertices = np.stack([np.abs(x - 4), y, x], axis=-1).reshape(-1, 3).astype(np.float64) + faces = [] + for row in range(size - 1): + for col in range(size - 1): + first = row * size + col + faces.append((first, first + size, first + 1)) + faces.append((first + 1, first + size, first + size + 1)) + faces = np.asarray(faces, dtype=np.int64) + + unprotected = bic.mesh.simplify_mesh( + vertices, faces, 0.65, feature_angle=20.0, feature_weight=0.0 + ) + protected = bic.mesh.simplify_mesh( + vertices, faces, 0.65, feature_angle=20.0, feature_weight=100.0 + ) + def fold_error(points): + return np.max(np.abs(points[:, 0] - np.abs(points[:, 2] - 4.0))) + + assert not np.array_equal(unprotected[0], protected[0]) + assert fold_error(protected[0]) <= fold_error(unprotected[0]) + + +@pytest.mark.parametrize("reduction", [-0.1, 1.0, np.nan]) +def test_rejects_invalid_reduction(reduction): + vertices, faces = _tetrahedron() + with pytest.raises(ValueError, match="reduction"): + bic.mesh.simplify_mesh(vertices, faces, reduction) + + +def test_rejects_bad_values_and_feature_options(): + vertices, faces = _tetrahedron() + with pytest.raises(ValueError, match="values must have shape"): + bic.mesh.simplify_mesh(vertices, faces, 0.2, values=np.zeros(3)) + with pytest.raises(TypeError, match="values must have dtype"): + bic.mesh.simplify_mesh(vertices, faces, 0.2, values=np.zeros(4, np.int32)) + with pytest.raises(ValueError, match="feature_angle"): + bic.mesh.simplify_mesh(vertices, faces, 0.2, feature_angle=181) + with pytest.raises(ValueError, match="feature_weight"): + bic.mesh.simplify_mesh(vertices, faces, 0.2, feature_weight=-1) + + +def test_rejects_degenerate_inconsistently_oriented_and_nonmanifold_meshes(): + vertices, faces = _tetrahedron() + + degenerate = faces.copy() + degenerate[0] = [0, 0, 1] + with pytest.raises(ValueError, match="repeated vertex"): + bic.mesh.simplify_mesh(vertices, degenerate, 0.2) + + inconsistent = faces.copy() + inconsistent[0] = inconsistent[0, ::-1] + with pytest.raises(ValueError, match="consistent winding"): + bic.mesh.simplify_mesh(vertices, inconsistent, 0.2) + + nonmanifold_vertices = np.vstack([vertices, [[0.0, -1.0, 0.0]]]) + nonmanifold_faces = np.vstack([faces, [[0, 1, 4]]]) + with pytest.raises(ValueError, match="more than two faces"): + bic.mesh.simplify_mesh(nonmanifold_vertices, nonmanifold_faces, 0.2) diff --git a/tests/test_utils_mesh_smoothing.py b/tests/mesh/test_smoothing.py similarity index 77% rename from tests/test_utils_mesh_smoothing.py rename to tests/mesh/test_smoothing.py index 7f3b46e..00e080b 100644 --- a/tests/test_utils_mesh_smoothing.py +++ b/tests/mesh/test_smoothing.py @@ -7,6 +7,11 @@ import bioimage_cpp as bic +def test_public_namespace(): + assert hasattr(bic.mesh, "smooth_mesh") + assert not hasattr(bic.utils, "smooth_mesh") + + def _tetrahedron(dtype): verts = np.array( [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]], @@ -50,7 +55,7 @@ def _octahedron(): @pytest.mark.parametrize("dtype", [np.float32, np.float64]) def test_tetrahedron_one_iteration_collapses_to_centroid(dtype): verts, normals, faces = _tetrahedron(dtype) - out_verts, out_normals = bic.utils.smooth_mesh(verts, normals, faces, iterations=1) + out_verts, out_normals = bic.mesh.smooth_mesh(verts, normals, faces, iterations=1) assert out_verts.dtype == verts.dtype assert out_normals.dtype == normals.dtype @@ -66,7 +71,7 @@ def test_tetrahedron_one_iteration_collapses_to_centroid(dtype): def test_iterations_zero_returns_copy_of_inputs(): verts, normals, faces = _tetrahedron(np.float64) - out_verts, out_normals = bic.utils.smooth_mesh(verts, normals, faces, iterations=0) + out_verts, out_normals = bic.mesh.smooth_mesh(verts, normals, faces, iterations=0) np.testing.assert_array_equal(out_verts, verts) np.testing.assert_array_equal(out_normals, normals) @@ -78,7 +83,7 @@ def test_iterations_zero_returns_copy_of_inputs(): def test_octahedron_converges_to_centroid(): verts, faces = _octahedron() normals = verts.copy() - out_verts, _ = bic.utils.smooth_mesh(verts, normals, faces, iterations=100) + out_verts, _ = bic.mesh.smooth_mesh(verts, normals, faces, iterations=100) centroid = verts.mean(axis=0) # Every vertex should be near the centroid after enough iterations. np.testing.assert_allclose(out_verts, np.broadcast_to(centroid, verts.shape), atol=1e-6) @@ -94,7 +99,7 @@ def test_non_3d_feature_vectors(dim): [[0, 1, 2], [1, 2, 3], [4, 5, 6], [5, 6, 7], [0, 4, 7], [2, 3, 5]], dtype=np.int64, ) - out_verts, out_normals = bic.utils.smooth_mesh(verts, normals, faces, iterations=3) + out_verts, out_normals = bic.mesh.smooth_mesh(verts, normals, faces, iterations=3) assert out_verts.shape == (n_verts, dim) assert out_normals.shape == (n_verts, dim) @@ -104,8 +109,8 @@ def test_parity_with_iterations_alternation(): verts, normals, faces = _tetrahedron(np.float64) cur_v, cur_n = verts, normals for _ in range(5): - cur_v, cur_n = bic.utils.smooth_mesh(cur_v, cur_n, faces, iterations=1) - multi_v, multi_n = bic.utils.smooth_mesh(verts, normals, faces, iterations=5) + cur_v, cur_n = bic.mesh.smooth_mesh(cur_v, cur_n, faces, iterations=1) + multi_v, multi_n = bic.mesh.smooth_mesh(verts, normals, faces, iterations=5) np.testing.assert_allclose(multi_v, cur_v, rtol=1e-12) np.testing.assert_allclose(multi_n, cur_n, rtol=1e-12) @@ -121,9 +126,9 @@ def test_threading_matches_serial(): axis=1, ).astype(np.int64) - serial = bic.utils.smooth_mesh(verts, normals, faces, iterations=8, n_threads=1) - parallel = bic.utils.smooth_mesh(verts, normals, faces, iterations=8, n_threads=4) - auto = bic.utils.smooth_mesh(verts, normals, faces, iterations=8, n_threads=0) + serial = bic.mesh.smooth_mesh(verts, normals, faces, iterations=8, n_threads=1) + parallel = bic.mesh.smooth_mesh(verts, normals, faces, iterations=8, n_threads=4) + auto = bic.mesh.smooth_mesh(verts, normals, faces, iterations=8, n_threads=0) np.testing.assert_array_equal(serial[0], parallel[0]) np.testing.assert_array_equal(serial[1], parallel[1]) np.testing.assert_array_equal(serial[0], auto[0]) @@ -139,17 +144,17 @@ def test_non_contiguous_inputs_are_accepted(): normals_view = big_normals[::2] assert not verts_view.flags.c_contiguous - out_verts, out_normals = bic.utils.smooth_mesh(verts_view, normals_view, faces, iterations=1) - expected_v, expected_n = bic.utils.smooth_mesh(verts, normals, faces, iterations=1) + out_verts, out_normals = bic.mesh.smooth_mesh(verts_view, normals_view, faces, iterations=1) + expected_v, expected_n = bic.mesh.smooth_mesh(verts, normals, faces, iterations=1) np.testing.assert_array_equal(out_verts, expected_v) np.testing.assert_array_equal(out_normals, expected_n) def test_faces_dtype_is_normalised_to_int64(): verts, normals, faces = _tetrahedron(np.float64) - out_int64, _ = bic.utils.smooth_mesh(verts, normals, faces.astype(np.int64), iterations=1) - out_int32, _ = bic.utils.smooth_mesh(verts, normals, faces.astype(np.int32), iterations=1) - out_uint32, _ = bic.utils.smooth_mesh(verts, normals, faces.astype(np.uint32), iterations=1) + out_int64, _ = bic.mesh.smooth_mesh(verts, normals, faces.astype(np.int64), iterations=1) + out_int32, _ = bic.mesh.smooth_mesh(verts, normals, faces.astype(np.int32), iterations=1) + out_uint32, _ = bic.mesh.smooth_mesh(verts, normals, faces.astype(np.uint32), iterations=1) np.testing.assert_array_equal(out_int64, out_int32) np.testing.assert_array_equal(out_int64, out_uint32) @@ -159,7 +164,7 @@ def test_rejects_mismatched_shapes(): normals_bad = np.zeros((4, 2), dtype=np.float64) faces = np.array([[0, 1, 2]], dtype=np.int64) with pytest.raises(ValueError, match="same shape as verts"): - bic.utils.smooth_mesh(verts, normals_bad, faces, iterations=1) + bic.mesh.smooth_mesh(verts, normals_bad, faces, iterations=1) def test_rejects_1d_verts(): @@ -167,33 +172,33 @@ def test_rejects_1d_verts(): normals = verts.copy() faces = np.array([[0, 1, 2]], dtype=np.int64) with pytest.raises(ValueError, match="verts must have ndim=2"): - bic.utils.smooth_mesh(verts, normals, faces, iterations=1) + bic.mesh.smooth_mesh(verts, normals, faces, iterations=1) def test_rejects_bad_faces_shape(): verts, normals, _ = _tetrahedron(np.float64) bad_faces = np.array([[0, 1, 2, 3]], dtype=np.int64) with pytest.raises(ValueError, match=r"faces must have shape"): - bic.utils.smooth_mesh(verts, normals, bad_faces, iterations=1) + bic.mesh.smooth_mesh(verts, normals, bad_faces, iterations=1) def test_rejects_out_of_range_faces(): verts, normals, _ = _tetrahedron(np.float64) bad_faces = np.array([[0, 1, 10]], dtype=np.int64) with pytest.raises(ValueError, match=r"\[0, n_verts\)"): - bic.utils.smooth_mesh(verts, normals, bad_faces, iterations=1) + bic.mesh.smooth_mesh(verts, normals, bad_faces, iterations=1) def test_rejects_negative_iterations(): verts, normals, faces = _tetrahedron(np.float64) with pytest.raises(ValueError, match="iterations must be non-negative"): - bic.utils.smooth_mesh(verts, normals, faces, iterations=-1) + bic.mesh.smooth_mesh(verts, normals, faces, iterations=-1) def test_rejects_negative_n_threads(): verts, normals, faces = _tetrahedron(np.float64) with pytest.raises(ValueError, match="n_threads must be non-negative"): - bic.utils.smooth_mesh(verts, normals, faces, iterations=1, n_threads=-1) + bic.mesh.smooth_mesh(verts, normals, faces, iterations=1, n_threads=-1) def test_rejects_unsupported_verts_dtype(): @@ -201,7 +206,7 @@ def test_rejects_unsupported_verts_dtype(): normals = verts.copy() faces = np.array([[0, 1, 2]], dtype=np.int64) with pytest.raises(TypeError, match="verts must have one of dtypes"): - bic.utils.smooth_mesh(verts, normals, faces, iterations=1) + bic.mesh.smooth_mesh(verts, normals, faces, iterations=1) def test_rejects_mismatched_verts_normals_dtype(): @@ -209,7 +214,7 @@ def test_rejects_mismatched_verts_normals_dtype(): normals = np.zeros((4, 3), dtype=np.float64) faces = np.array([[0, 1, 2]], dtype=np.int64) with pytest.raises(TypeError, match="same dtype"): - bic.utils.smooth_mesh(verts, normals, faces, iterations=1) + bic.mesh.smooth_mesh(verts, normals, faces, iterations=1) def test_parity_with_python_reference(): @@ -223,8 +228,8 @@ def test_parity_with_python_reference(): """ pytest.importorskip("nifty") - repo_root = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) - dev_path = os.path.join(repo_root, "development", "utils") + repo_root = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) + dev_path = os.path.join(repo_root, "development", "mesh") sys.path.insert(0, dev_path) try: from _mesh_smoothing_reference import smooth_mesh as smooth_mesh_reference @@ -241,7 +246,7 @@ def test_parity_with_python_reference(): faces = hull.simplices.astype(np.int64) normals = verts.copy() - out_v, out_n = bic.utils.smooth_mesh(verts, normals, faces, iterations=1) + out_v, out_n = bic.mesh.smooth_mesh(verts, normals, faces, iterations=1) ref_v, ref_n = smooth_mesh_reference(verts, normals, faces, 1) np.testing.assert_allclose(out_v, ref_v, rtol=1e-10, atol=1e-12) np.testing.assert_allclose(out_n, ref_n, rtol=1e-10, atol=1e-12)