Merge pull request #390 from campaslab/revert-389-use-sorted-semi-axe…

…s-for-curvature-calculation

Revert "Use sorted semi axes for curvature calculation"

src/napari_stress/_approximation/expansion.py

@@ -554,11 +554,11 @@ def properties(self):
             The maximum and minimum curvatures :math:`H_{max}` and :math:`H_{min}` are calculated as follows:
 
             .. math::
-                H_{max} = a / (2 * c^2) + a / (2 * b^2)
+                H_{max} = 1 / (2 * a^2) + 1 / (2 * b^2)
 
-                H_{min} = c / (2 * b^2) + c / (2 * a^2)
+                H_{min} = 1 / (2 * b^2) + 1 / (2 * a^2)
 
-            where a, b and c are the lengths of the ellipsoid axes along the three spatial dimensions.
+            where a, b are the largest and smallest axes of the ellipsoid, respectively.
         """
         return self._properties
 
@@ -571,14 +571,21 @@ def _measure_max_min_curvatures(self):
         None
 
         """
-        semi_axis_sorted = np.sort(self._axes)
-        a = semi_axis_sorted[2]
-        b = semi_axis_sorted[1]
-        c = semi_axis_sorted[0]
+        # get and remove the largest, smallest and medial axis
+        axes = list(self._axes)
+        largest_axis = max(axes)
+        axes.remove(largest_axis)
+        smallest_axis = min(axes)
+        axes.remove(smallest_axis)
+        medial_axis = axes[0]
 
         # accoording to paper (https://www.biorxiv.org/content/10.1101/2021.03.26.437148v1.full)
-        maximum_mean_curvature = a / (2 * c**2) + a / (2 * b**2)
-        minimum_mean_curvature = c / (2 * b**2) + c / (2 * a**2)
+        maximum_mean_curvature = largest_axis / (
+            2 * smallest_axis**2
+        ) + largest_axis / (2 * medial_axis**2)
+        minimum_mean_curvature = smallest_axis / (
+            2 * medial_axis**2
+        ) + smallest_axis / (2 * largest_axis**2)
 
         self._properties["maximum_mean_curvature"] = maximum_mean_curvature
         self._properties["minimum_mean_curvature"] = minimum_mean_curvature

src/napari_stress/_tests/test_approximation.py

@@ -124,23 +124,6 @@ def test_lsq_ellipsoid0(n_tests=10):
         max_mean_curvatures.append(expander.properties["maximum_mean_curvature"])
         min_mean_curvatures.append(expander.properties["minimum_mean_curvature"])
 
-        # The mean curvature of an ellipsoid is given by the formula
-        # H = a / (2 * c^2) + a / (2 * b^2) for the maximum mean curvature
-        # and H = c / (2 * b^2) + c / (2 * a^2) for the minimum mean curvature
-        # where a, b, c are the semi-axes of the ellipsoid
-        axes_sorted = np.sort(expander.axes_)
-        a = axes_sorted[2]
-        b = axes_sorted[1]
-        c = axes_sorted[0]
-        h_max_theory = a / (2 * c**2) + a / (2 * b**2)
-        h_min_theory = c / (2 * b**2) + c / (2 * a**2)
-        assert expander.properties["maximum_mean_curvature"] == h_max_theory
-        assert expander.properties["minimum_mean_curvature"] == h_min_theory
-        assert (
-            expander.properties["maximum_mean_curvature"]
-            > expander.properties["minimum_mean_curvature"]
-        )
-
         assert np.allclose(center, x0)
         assert np.allclose(a2, axes_lengths[2])
         assert np.allclose(a1, axes_lengths[1])