Add example scripts for geodesic distances

Author: constantinpape

examples/distance/complex-shape.tif

@@ -0,0 +1,21 @@
+import imageio.v3 as imageio
+import napari
+import numpy as np
+
+from bioimage_cpp.distance import geodesic_distance_field, distance_transform
+
+mask = (imageio.imread("./complex-shape.tif") == 2).astype("uint8")
+# compute the most central point --- distance transform maximum
+max_point = np.argmax(distance_transform(mask).ravel())
+max_point = np.unravel_index(max_point, mask.shape)
+
+sources = np.array([max_point])
+distance_field, gradients = geodesic_distance_field(mask, sources, return_gradient=True)
+gradients = gradients.transpose((2, 0, 1))
+
+v = napari.Viewer()
+v.add_image(distance_field)
+v.add_image(gradients)
+v.add_labels(mask)
+v.add_points([max_point])
+napari.run()

examples/distance/geodesic_2d.py

@@ -0,0 +1,20 @@
+import napari
+import numpy as np
+
+from bioimage_cpp.distance import geodesic_distance_field, distance_transform
+from shape import make_shape
+
+mask = make_shape()
+max_point = np.argmax(distance_transform(mask).ravel())
+max_point = np.unravel_index(max_point, mask.shape)
+
+sources = np.array([max_point])
+distance_field, gradients = geodesic_distance_field(mask, sources, return_gradient=True)
+gradients = gradients.transpose((3, 0, 1, 2))
+
+v = napari.Viewer()
+v.add_image(distance_field)
+v.add_image(gradients)
+v.add_labels(mask)
+v.add_points([max_point])
+napari.run()

examples/distance/geodesic_3d.py

@@ -0,0 +1,36 @@
+import napari
+import numpy as np
+from scipy.spatial import cKDTree
+from skimage.measure import marching_cubes
+
+from bioimage_cpp.distance import geodesic_distance_field_mesh
+from shape import make_shape
+
+# Create the mesh based on our sample shape.
+mask = make_shape()
+verts, faces, _, _ = marching_cubes(mask, level=0.5)
+
+# Create some points in the volume as sources for the distances.
+c = mask.shape[0] / 2.0
+angles = np.linspace(0.0, 2.0 * np.pi, num=5, endpoint=False)
+points = np.array(
+    [[c, c + 33.0 * np.sin(a), c + 33.0 * np.cos(a)] for a in angles]
+)
+
+# Snap the points to the closest mesh vertex.
+_, idx = cKDTree(verts).query(points)
+sources = np.unique(idx)
+
+# Compute the distance field on the mesh.
+field = geodesic_distance_field_mesh(verts, faces, sources)
+
+# Alternatively, you can also compute the pairwise geodesic distance
+# between all the points like this:
+# pairwise_distances = geodesic_distances_mesh(verts, faces, sources, number_of_threads=1)
+
+# Display the mesh and points in napari.
+# The mesh gets the distance field as values.
+v = napari.Viewer()
+v.add_surface((verts, faces, field))
+v.add_points(points)
+napari.run()

examples/distance/geodesic_mesh.py

@@ -0,0 +1,10 @@
+import numpy as np
+
+
+def make_shape(size=80):
+    R, r = 24.0, 9.0
+    c = size / 2.0
+    zz, yy, xx = np.ogrid[:size, :size, :size]
+    q = np.sqrt((yy - c) ** 2 + (xx - c) ** 2) - R
+    mask = (q ** 2 + (zz - c) ** 2 <= r ** 2).astype("uint8")
+    return mask