improve efficiency of filters.frangi and filters.hessian by directly …

…obtaining the eigenvalues sorted by magnitude

python/cucim/src/cucim/skimage/feature/corner.py

@@ -634,14 +634,18 @@ def _image_orthogonal_matrix33_eigvals(
     return eigs
 
 
-def _symmetric_compute_eigenvalues(S_elems, fast_3d=True):
+def _symmetric_compute_eigenvalues(S_elems, sort='descending', abs_sort=False):
     """Compute eigenvalues from the upperdiagonal entries of a symmetric matrix
 
     Parameters
     ----------
     S_elems : list of ndarray
         The upper-diagonal elements of the matrix, as returned by
         `hessian_matrix` or `structure_tensor`.
+    sort : {"ascending", "descending"}, optional
+        Eigenvalues should be sorted in the specified order.
+    abs_sort : boolean, optional
+        If ``True``, sort based on the absolute values.
 
     Returns
     -------
@@ -653,19 +657,24 @@ def _symmetric_compute_eigenvalues(S_elems, fast_3d=True):
 
     if len(S_elems) == 3:  # Use fast analytical kernel for 2D
         eigs = _image_orthogonal_matrix22_eigvals(
-            *S_elems, sort='descending', abs_sort=False
+            *S_elems, sort=sort, abs_sort=abs_sort
         )
-    elif fast_3d and len(S_elems) == 6:  # Use fast analytical kernel for 3D
+    elif len(S_elems) == 6:  # Use fast analytical kernel for 3D
         eigs = _image_orthogonal_matrix33_eigvals(
-            *S_elems, sort='descending', abs_sort=False
+            *S_elems, sort=sort, abs_sort=abs_sort
         )
     else:
         # n-dimensional case. warning: extremely memory inefficient!
         matrices = _symmetric_image(S_elems)
         # eigvalsh returns eigenvalues in increasing order. We want decreasing
-        eigs = cp.linalg.eigvalsh(matrices)[..., ::-1]
+        eigs = cp.linalg.eigvalsh(matrices)
         leading_axes = tuple(range(eigs.ndim - 1))
         eigs = cp.transpose(eigs, (eigs.ndim - 1,) + leading_axes)
+        if abs_sort:
+            # (sort by magnitude)
+            eigs = cp.take_along_axis(eigs, cp.abs(eigs).argsort(0), 0)
+        if sort == 'descending':
+            eigs = eigs[::-1, ...]
     return eigs
 
 

python/cucim/src/cucim/skimage/filters/ridges.py

@@ -15,7 +15,8 @@ class of ridge filters relies on the eigenvalues of the Hessian matrix of
 import numpy as np
 
 from .._shared.utils import _supported_float_type, check_nD, deprecated
-from ..feature.corner import hessian_matrix, hessian_matrix_eigvals
+from ..feature.corner import (_symmetric_compute_eigenvalues, hessian_matrix,
+                              hessian_matrix_eigvals)
 from ..util import img_as_float
 
 
@@ -479,12 +480,12 @@ def frangi(image, sigmas=range(1, 10, 2), scale_range=None,
     for i, sigma in enumerate(sigmas):  # Filter for all sigmas.
         H = hessian_matrix(image, sigma, mode=mode, cval=cval,
                            use_gaussian_derivatives=True)
-        eigvals = hessian_matrix_eigvals(H)
 
-        # Sort eigenvalues by ascending magnitude
-        # (hessian_matrix_eigvals are sorted in descending order, but not by
-        #  magnitude)
-        eigvals = cp.take_along_axis(eigvals, cp.abs(eigvals).argsort(0), 0)
+        # using _symmetric_compute_eigenvalues rather than hessian_matrix_eigvals
+        # so we can directly sort by ascending magnitude
+        eigvals = _symmetric_compute_eigenvalues(
+            H, sort='ascending', abs_sort=True
+        )
 
         # compute squared sum of the eigenvalues
         if ndim == 2: