dont recompute mesh but reuse it. if not provided return NaN in compute cristae junctions

Author: lufre1

synapse_net/cristae_analysis.py

@@ -608,24 +608,25 @@ def compute_junction_distances(
     the nearest vertex of a triangle mesh, and pairwise surface geodesics are computed with
     ``bioimage_cpp.distance.geodesic_distances_mesh``. The mesh is the **eroded-mito (lumen) surface**
     passed in as ``mesh_vertices``/``mesh_faces`` by :func:`_single_mito_row` (a clean, single-wall
-    surface at the membrane's inner edge); if none is supplied a mesh is built from ``membrane_mask``
-    as a convenience. When no usable surface mesh exists (empty membrane / degenerate mesh), the
-    junction distances are NaN.
+    surface at the membrane's inner edge). If no mesh is supplied — or no usable surface mesh exists
+    (empty membrane / degenerate mesh) — the junction distances are NaN. (There is no membrane-band
+    fallback mesh: the metric is defined on the lumen surface, and meshing the thick membrane band
+    would give a different, capped double-wall surface.)
 
     A Clark-Evans nearest-neighbour index summarises whether the junctions are clustered.
 
     Args:
         contact_labels: Integer junction label array (0 = background, 1..n = junctions),
             e.g. the first return value of :func:`detect_contact_sites`.
-        membrane_mask: Binary mitochondrial membrane mask the junctions sit on (used to build a
-            fallback mesh when ``mesh_vertices``/``mesh_faces`` are not supplied).
+        membrane_mask: Binary mitochondrial membrane mask the junctions sit on. Only used for the
+            empty-membrane early-out (no membrane → NaN); it is not meshed.
         voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys.
         surface_area_nm2: Membrane/mito surface area used as the reference area for the
             Clark-Evans expectation. If None or non-positive, the clustering index is NaN.
         n_jobs: 1 = serial, -1 = all cores (forwarded to the mesh solver's thread count).
         mesh_vertices: Optional (n_vertices, 3) mesh vertices in nm (unpadded mask frame; see
-            :func:`_surface_mesh`) — the eroded-mito (lumen) surface. If omitted, a mesh is built
-            from ``membrane_mask``.
+            :func:`_surface_mesh`) — the eroded-mito (lumen) surface. If omitted, the junction
+            distances are NaN.
         mesh_faces: Optional (n_faces, 3) triangle indices matching ``mesh_vertices``.
 
     Returns:
@@ -659,12 +660,9 @@ def compute_junction_distances(
     )
 
     if mesh_vertices is not None and mesh_faces is not None and len(mesh_faces) > 0:
-        mesh = (mesh_vertices, mesh_faces)
+        distance_matrix = _junction_matrix_mesh(centroids, sampling, mesh_vertices, mesh_faces, n_jobs)
     else:
-        mesh = _surface_mesh(membrane, sampling)
-    distance_matrix = (
-        _junction_matrix_mesh(centroids, sampling, mesh[0], mesh[1], n_jobs) if mesh is not None else None
-    )
+        distance_matrix = None
 
     if distance_matrix is None:
         summary = dict(_JUNCTION_DISTANCE_NAN)

test/test_cristae_analysis.py

@@ -570,12 +570,13 @@ def _flat_membrane_with_junctions(positions, shape=(12, 40, 40), z=6):
             labels[z, y, x] = i
         return labels, membrane
 
-    # These exercise compute_junction_distances directly; with no mesh supplied it meshes the given
-    # membrane and takes the surface geodesic, so they require the bioimage-cpp geodesic API.
+    # These exercise compute_junction_distances directly. The metric is defined on the lumen surface,
+    # so the mesh must be supplied explicitly (there is no membrane-band fallback); here the membrane's
+    # own surface is meshed as the test fixture. They require the bioimage-cpp geodesic API.
     def test_geodesic_follows_bent_membrane(self):
         # An L-shaped membrane: the geodesic around the bend is longer than the straight line
         # between the two seed voxels.
-        from synapse_net.cristae_analysis import compute_junction_distances
+        from synapse_net.cristae_analysis import compute_junction_distances, _surface_mesh
         shape = (5, 40, 40)
         membrane = np.zeros(shape, dtype=bool)
         z = 2
@@ -584,7 +585,10 @@ def test_geodesic_follows_bent_membrane(self):
         labels = np.zeros(shape, dtype=np.int32)
         labels[z, 5, 6] = 1               # near the far end of the horizontal arm
         labels[z, 33, 34] = 2             # near the far end of the vertical arm
-        dist, _ = compute_junction_distances(labels, membrane, voxel_size=1.0)
+        verts, faces = _surface_mesh(membrane, np.ones(3))
+        dist, _ = compute_junction_distances(
+            labels, membrane, voxel_size=1.0, mesh_vertices=verts, mesh_faces=faces
+        )
         straight = np.sqrt((33 - 5) ** 2 + (34 - 6) ** 2)
         self.assertGreater(dist[0, 1], straight * 1.2)
 
@@ -600,16 +604,21 @@ def test_fewer_than_two_junctions_is_nan(self):
     def test_clustered_index_lower_than_dispersed(self):
         # Same membrane/area and junction count, but tightly grouped vs evenly spread:
         # the clustered arrangement must give a smaller Clark-Evans index.
-        from synapse_net.cristae_analysis import compute_junction_distances
+        from synapse_net.cristae_analysis import compute_junction_distances, _surface_mesh
         area = 40.0 * 40.0
         clustered_pos = [(18, 18), (18, 20), (20, 18), (20, 20)]
         dispersed_pos = [(8, 8), (8, 30), (30, 8), (30, 30)]
-        _, clustered = compute_junction_distances(
-            *self._flat_membrane_with_junctions(clustered_pos), voxel_size=1.0, surface_area_nm2=area
-        )
-        _, dispersed = compute_junction_distances(
-            *self._flat_membrane_with_junctions(dispersed_pos), voxel_size=1.0, surface_area_nm2=area
-        )
+
+        def _cjd(positions):
+            labels, membrane = self._flat_membrane_with_junctions(positions)
+            verts, faces = _surface_mesh(membrane, np.ones(3))
+            return compute_junction_distances(
+                labels, membrane, voxel_size=1.0, surface_area_nm2=area,
+                mesh_vertices=verts, mesh_faces=faces,
+            )
+
+        _, clustered = _cjd(clustered_pos)
+        _, dispersed = _cjd(dispersed_pos)
         self.assertLess(clustered["junction_clustering_index"], dispersed["junction_clustering_index"])
 
 
@@ -936,18 +945,18 @@ def test_clipped_mito_open_mesh_finite(self):
         area_closed = _surface_area(mito_binary, np.ones(ndim))
         self.assertLess(area_open, area_closed)  # the fabricated z-caps are no longer counted
 
-    def test_primitive_mesh_on_flat_membrane(self):
-        # With no mesh supplied the primitive meshes the given membrane surface — finite, positive,
-        # and of the right order (junctions ~17-20 apart). Absolute accuracy on a 1-voxel sheet is not
-        # asserted (that degenerate mesh is only the fallback; the pipeline meshes the thicker lumen).
+    def test_no_mesh_supplied_is_nan(self):
+        # With no lumen mesh supplied the junction distances are NaN — there is no membrane-band
+        # fallback mesh (the metric is defined on the eroded-mito lumen surface).
         from synapse_net.cristae_analysis import compute_junction_distances
         labels, membrane = TestJunctionDistances._flat_membrane_with_junctions(
             [(20, 8), (20, 28), (8, 20), (32, 20)]
         )
         _, summary = compute_junction_distances(labels, membrane, 1.0, surface_area_nm2=1000.0)
-        mean_nn = summary["mean_nn_junction_distance_nm"]
-        self.assertTrue(np.isfinite(mean_nn) and mean_nn > 0)
-        self.assertLess(mean_nn, 100.0)  # sane order of magnitude, not a runaway path
+        self.assertEqual(summary["junction_count"], 4)
+        self.assertTrue(np.isnan(summary["mean_nn_junction_distance_nm"]))
+        self.assertTrue(np.isnan(summary["median_nn_junction_distance_nm"]))
+        self.assertTrue(np.isnan(summary["junction_clustering_index"]))
 
 
 _EXPECTED_COLUMNS = [