Merge pull request #375 from jo-mueller/replace-pygeodesics

Replace pygeodesics with potpourri3d and gdist

Author: jo-mueller

docs/04_FAQ/geodesics_tools.ipynb

@@ -68,6 +68,7 @@ parts:
       - file: 04_FAQ/validation/validation_stress_paper
       - file: 04_FAQ/validation/analyze_validation
       - file: 04_FAQ/validation/create_validation_data
+    - file: 04_FAQ/geodesics_tools
 
   - caption: API
     chapters:

docs/_toc.yml

@@ -33,7 +33,6 @@ dependencies = [
     "napari-tools-menu>=0.1.15",
     "numpy<1.24.0",
     "pandas",
-    "pygeodesic",
     "scikit-image",
     "scipy>=1.9.0",
     "seaborn",
@@ -42,6 +41,8 @@ dependencies = [
     "vedo>=2023.5.0",
     "vispy",
     "deprecation",
+    "gdist",
+    "potpourri3d",
     # "bokeh >= 3.1.0",
 ]
 

pyproject.toml

@@ -1,16 +1,13 @@
 from typing import List
 
 import numpy as np
-import tqdm
 from napari.types import LayerDataTuple, SurfaceData
 from napari_tools_menu import register_function
 
 from .._utils.frame_by_frame import frame_by_frame
 
 
-def geodesic_distance_matrix(
-    surface: SurfaceData, show_progress: bool = False
-) -> np.ndarray:
+def geodesic_distance_matrix(surface: SurfaceData) -> np.ndarray:
     """
     Calculate a pairwise distance matrix for vertices of a surface.
 
@@ -28,27 +25,16 @@ def geodesic_distance_matrix(
         at `distance_matrix[i, j]`
 
     """
-    from pygeodesic import geodesic
+    import gdist
+    from .._utils import sanitize_faces
 
-    geoalg = geodesic.PyGeodesicAlgorithmExact(surface[0], surface[1])
+    # Reorder the faces of the surface to ensure consistent orientation
+    sanitized_surface = sanitize_faces(surface)
+    vertices = sanitized_surface[0]
+    faces = sanitized_surface[1]
 
-    n_points = len(surface[0])
-    distance_matrix = np.zeros((n_points, n_points))
-    points = surface[0]
-
-    if show_progress:
-        iterator = tqdm.tqdm(
-            enumerate(points), desc="Calculating geodesic distances"
-        )
-    else:
-        iterator = enumerate(points)
-
-    for idx, pt in iterator:
-        distances, _ = geoalg.geodesicDistances(
-            [idx], np.arange(idx + 1, n_points)
-        )
-        distance_matrix[idx, idx + 1 :] = distances
-        distance_matrix[idx + 1 :, idx] = distances
+    distance_matrix = gdist.local_gdist_matrix(
+        vertices, faces, max_distance=1e9).toarray()
 
     return distance_matrix
 
@@ -77,10 +63,15 @@ def geodesic_path(
         coordinates, the second dimension is the vector from point to point.
 
     """
-    from pygeodesic import geodesic
+    import potpourri3d as pp3d
+    from .._utils import sanitize_faces
 
-    geoalg = geodesic.PyGeodesicAlgorithmExact(surface[0], surface[1])
-    distances, path = geoalg.geodesicDistance(index_1, index_2)
+    sanitized_surface = sanitize_faces(surface)
+    vertices = sanitized_surface[0]
+    faces = sanitized_surface[1]
+
+    path_solver = pp3d.EdgeFlipGeodesicSolver(vertices, faces)
+    path = path_solver.find_geodesic_path(index_1, index_2)
 
     # convert points to vectors from point to point
     vectors = []
@@ -240,8 +231,10 @@ def local_extrema_analysis(
             - `min_max_pair_distances`: Distances between all pairs of local minima and maxima.
             - `min_max_pair_anisotropies`: Difference in input value `(vertices, faces, values)` between all pairs of local minima and maxima.
     """
-    triangles = surface[1]
+    from .._utils import sanitize_faces
     feature = surface[2]
+    surface = sanitize_faces(surface)
+    triangles = surface[1]
 
     quad_fit = len(feature)
 

src/napari_stress/_measurements/geodesics.py

@@ -595,7 +595,7 @@ def comprehensive_analysis(
     GDM = None
     if GDM is None:
         GDM = measurements.geodesic_distance_matrix(
-            surface_cell_stress, show_progress=verbose
+            surface_cell_stress
         )
 
     if maximal_distance is None:

src/napari_stress/_measurements/toolbox.py

@@ -1,5 +1,4 @@
 import numpy as np
-from napari.types import PointsData, SurfaceData
 from napari_tools_menu import register_function
 
 from .._utils.frame_by_frame import frame_by_frame
@@ -10,8 +9,8 @@
 )
 @frame_by_frame
 def reconstruct_surface_from_quadrature_points(
-    points: PointsData,
-) -> SurfaceData:
+    points: "napari.types.PointsData",
+) -> "napari.types.SurfaceData":
     """
     Reconstruct the surface for a given set of quadrature points.
 
@@ -22,12 +21,12 @@ def reconstruct_surface_from_quadrature_points(
 
     Returns
     -------
-    Sphere_Triangulation_Array : TYPE
-        DESCRIPTION.
+    tuple
+        Tuple of points and faces
 
     """
     from scipy.spatial import Delaunay
-
+    from .._utils import sanitize_faces
     from .._stress import lebedev_write_SPB as lebedev_write
 
     n_quadrature_points = len(points)
@@ -53,4 +52,4 @@ def reconstruct_surface_from_quadrature_points(
                 delauney_triangles[tri_i, vert_ind] = vertex
                 vert_ind = vert_ind + 1
 
-    return (points, delauney_triangles.astype(int))
+    return sanitize_faces((points, delauney_triangles.astype(int)))

src/napari_stress/_reconstruction/reconstruct_surface.py

@@ -823,9 +823,3 @@ def test_stresses():
         H_i, H0, H_i_ellipsoid, H0_ellipsoid, gamma
     )
     measurements.maximal_tissue_anisotropy(ellipsoid, gamma=gamma)
-
-
-if __name__ == "__main__":
-    import napari
-
-    test_comprehensive_stress_toolbox_4d(napari.Viewer)

src/napari_stress/_tests/test_measurements.py

@@ -1,11 +1,13 @@
 from ._aggregate_measurements import compile_data_from_layers
 from .frame_by_frame import TimelapseConverter, frame_by_frame
 from .import_export_settings import export_settings, import_settings
+from ._utils import sanitize_faces
 
 __all__ = [
     "import_settings",
     "export_settings",
     "compile_data_from_layers",
     "TimelapseConverter",
     "frame_by_frame",
+    "sanitize_faces",
 ]

src/napari_stress/_utils/__init__.py

@@ -0,0 +1,38 @@
+import numpy as np
+
+def sanitize_faces(surface: "napari.types.SurfaceData") -> "napari.types.SurfaceData":
+    """
+    Sanitize the faces of a surface to ensure they are in the correct format.
+
+    This function ensures that the faces of the surface are in the correct
+    format and order. It also ensures that the vertices are in the correct
+    format.
+
+    Parameters
+    ----------
+    surface : 'napari.types.SurfaceData'
+        The surface data to sanitize.
+
+    Returns
+    -------
+    'napari.types.SurfaceData'
+        The sanitized surface data.
+    """
+    # Ensure faces are in the correct format
+    faces = surface[1]
+    vertices = surface[0]
+    center = np.mean(vertices, axis=0)
+
+    for idx, triangle in enumerate(faces):
+        v1 = vertices[int(triangle[0])]
+        v2 = vertices[int(triangle[1])]
+        v3 = vertices[int(triangle[2])]
+
+        normal = np.cross(v2 - v1, v3 - v1)
+        normal = normal / np.linalg.norm(normal)
+        center_triangle = (v1 + v2 + v3) / 3
+
+        if np.dot(normal, center_triangle - center) < 0:
+            faces[idx] = faces[idx][::-1]
+
+    return (vertices.astype(np.float64), faces.astype(np.int32))