Refactor mesh smoothing function to the mesh header and module

Author: constantinpape

MIGRATION_GUIDE.md

@@ -1933,6 +1933,16 @@ 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`.
+
+
+
 ### Anti-Aliased Resampling
 
 `affine_transform` itself never pre-smooths the input; downsampling without

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

…pment/utils/_mesh_smoothing_reference.py …opment/mesh/_mesh_smoothing_reference.pydevelopment/utils/_mesh_smoothing_reference.py renamed to development/mesh/_mesh_smoothing_reference.py development/utils/_mesh_smoothing_reference.py renamed to development/mesh/_mesh_smoothing_reference.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): "

…opment/utils/benchmark_mesh_smoothing.py …lopment/mesh/benchmark_mesh_smoothing.pydevelopment/utils/benchmark_mesh_smoothing.py renamed to development/mesh/benchmark_mesh_smoothing.py development/utils/benchmark_mesh_smoothing.py renamed to development/mesh/benchmark_mesh_smoothing.py

@@ -12,9 +12,9 @@
 #include <utility>
 #include <vector>
 
-namespace bioimage_cpp {
+namespace bioimage_cpp::mesh {
 
-namespace detail::mesh_smoothing {
+namespace detail::smoothing {
 
 template <class I>
 struct Adjacency {
@@ -74,7 +74,7 @@ Adjacency<I> build_adjacency(const ConstArrayView<I> &faces, std::ptrdiff_t n_ve
     return Adjacency<I>{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<I>(faces, n_verts);
+    const auto adjacency = detail::smoothing::build_adjacency<I>(faces, n_verts);
 
     std::vector<V> scratch_verts(n_total);
     std::vector<V> 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<std::size_t>(n_threads),
         static_cast<std::size_t>(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<std::size_t>(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

include/bioimage_cpp/mesh_smoothing.hxx include/bioimage_cpp/mesh/smoothing.hxxinclude/bioimage_cpp/mesh_smoothing.hxx renamed to include/bioimage_cpp/mesh/smoothing.hxx

@@ -3,9 +3,11 @@
 #include "bioimage_cpp/array_view.hxx"
 #include "bioimage_cpp/detail/profile.hxx"
 #include "bioimage_cpp/mesh/marching_cubes.hxx"
+#include "bioimage_cpp/mesh/smoothing.hxx"
 
 #include <nanobind/ndarray.h>
 #include <nanobind/stl/optional.h>
+#include <nanobind/stl/pair.h>
 
 #include <cstddef>
 #include <cstdint>
@@ -26,6 +28,12 @@ using UInt8Input = nb::ndarray<nb::numpy, const std::uint8_t, nb::c_contig>;
 using FloatOutput = nb::ndarray<nb::numpy, float, nb::c_contig>;
 using Int32Output = nb::ndarray<nb::numpy, std::int32_t, nb::c_contig>;
 
+template <class T>
+using InputArray = nb::ndarray<nb::numpy, const T, nb::c_contig>;
+
+template <class T>
+using OutputArray = nb::ndarray<nb::numpy, T, nb::c_contig>;
+
 template <class T, class Output>
 Output output_array(std::vector<T> &&values, const std::vector<std::size_t> &shape) {
     if (values.empty()) {
@@ -58,6 +66,77 @@ std::vector<std::ptrdiff_t> shape_of(const FloatInput &array) {
     return shape;
 }
 
+using FacesArray = nb::ndarray<nb::numpy, const std::int64_t, nb::c_contig>;
+
+template <class V>
+std::pair<OutputArray<V>, OutputArray<V>> smooth_mesh_t(
+    InputArray<V> verts,
+    InputArray<V> 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<std::size_t> verts_ndarray_shape{n_verts_size, dim_size};
+    std::vector<std::ptrdiff_t> verts_view_shape{
+        static_cast<std::ptrdiff_t>(n_verts_size),
+        static_cast<std::ptrdiff_t>(dim_size),
+    };
+    std::vector<std::ptrdiff_t> faces_view_shape{
+        static_cast<std::ptrdiff_t>(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<V *>(p); });
+    auto *normals_data = new V[n_total]();
+    nb::capsule normals_owner(normals_data, [](void *p) noexcept { delete[] static_cast<V *>(p); });
+
+    ConstArrayView<V> verts_view{verts.data(), verts_view_shape, {}};
+    ConstArrayView<V> normals_view{normals.data(), verts_view_shape, {}};
+    ConstArrayView<std::int64_t> faces_view{faces.data(), faces_view_shape, {}};
+    ArrayView<V> out_verts_view{verts_data, verts_view_shape, {}};
+    ArrayView<V> out_normals_view{normals_data, verts_view_shape, {}};
+
+    {
+        nb::gil_scoped_release release;
+        mesh::smooth_mesh<V, std::int64_t>(
+            verts_view,
+            normals_view,
+            faces_view,
+            iterations,
+            n_threads,
+            out_verts_view,
+            out_normals_view
+        );
+    }
+
+    OutputArray<V> out_verts(
+        verts_data, verts_ndarray_shape.size(), verts_ndarray_shape.data(), verts_owner
+    );
+    OutputArray<V> out_normals(
+        normals_data, verts_ndarray_shape.size(), verts_ndarray_shape.data(), normals_owner
+    );
+    return {std::move(out_verts), std::move(out_normals)};
+}
+
 nb::tuple marching_cubes_float32(
     FloatInput volume,
     const double level,
@@ -150,6 +229,26 @@ nb::tuple marching_cubes_float32(
 } // namespace
 
 void bind_mesh(nb::module_ &m) {
+    m.def(
+        "_smooth_mesh_float32",
+        &smooth_mesh_t<float>,
+        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<double>,
+        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,

src/bindings/mesh.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 <nanobind/ndarray.h>
-#include <nanobind/stl/pair.h>
 #include <nanobind/stl/string.h>
 
 #include <algorithm>
@@ -82,76 +80,6 @@ OutputArray<T> take_dict_t(const nb::dict &relabeling, InputArray<T> to_relabel)
     return OutputArray<T>(data, ndarray_shape.size(), ndarray_shape.data(), owner);
 }
 
-using FacesArray = nb::ndarray<nb::numpy, const std::int64_t, nb::c_contig>;
-
-template <class V>
-std::pair<OutputArray<V>, OutputArray<V>> smooth_mesh_t(
-    InputArray<V> verts,
-    InputArray<V> 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<std::size_t> verts_ndarray_shape{n_verts_size, dim_size};
-    std::vector<std::ptrdiff_t> verts_view_shape{
-        static_cast<std::ptrdiff_t>(n_verts_size),
-        static_cast<std::ptrdiff_t>(dim_size),
-    };
-    std::vector<std::ptrdiff_t> faces_view_shape{
-        static_cast<std::ptrdiff_t>(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<V *>(p); });
-    auto *normals_data = new V[n_total]();
-    nb::capsule normals_owner(normals_data, [](void *p) noexcept { delete[] static_cast<V *>(p); });
-
-    ConstArrayView<V> verts_view{verts.data(), verts_view_shape, {}};
-    ConstArrayView<V> normals_view{normals.data(), verts_view_shape, {}};
-    ConstArrayView<std::int64_t> faces_view{faces.data(), faces_view_shape, {}};
-    ArrayView<V> out_verts_view{verts_data, verts_view_shape, {}};
-    ArrayView<V> out_normals_view{normals_data, verts_view_shape, {}};
-
-    {
-        nb::gil_scoped_release release;
-        smooth_mesh<V, std::int64_t>(
-            verts_view,
-            normals_view,
-            faces_view,
-            iterations,
-            n_threads,
-            out_verts_view,
-            out_normals_view
-        );
-    }
-
-    OutputArray<V> out_verts(
-        verts_data, verts_ndarray_shape.size(), verts_ndarray_shape.data(), verts_owner
-    );
-    OutputArray<V> out_normals(
-        normals_data, verts_ndarray_shape.size(), verts_ndarray_shape.data(), normals_owner
-    );
-    return {std::move(out_verts), std::move(out_normals)};
-}
 
 template <class T>
 OutputArray<std::int64_t> compute_rle_t(InputArray<T> 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<float>,
-        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<double>,
-        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<bool>, nb::arg("mask"),
           "COCO-style binary run-length encoding of a contiguous bool array.");
     m.def("_compute_rle_uint8", &compute_rle_t<std::uint8_t>, nb::arg("mask"),

src/bindings/utils.cxx

@@ -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.