Add function to directly compute eigenvalues of the Hessian matrix

diptorch/_modidx.py

@@ -11,7 +11,8 @@
                                   'diptorch.filters._hessian_3d': ('filters.html#_hessian_3d', 'diptorch/filters.py'),
                                   'diptorch.filters._hessian_as_matrix': ('filters.html#_hessian_as_matrix', 'diptorch/filters.py'),
                                   'diptorch.filters.gaussian_filter': ('filters.html#gaussian_filter', 'diptorch/filters.py'),
-                                  'diptorch.filters.hessian': ('filters.html#hessian', 'diptorch/filters.py')},
+                                  'diptorch.filters.hessian': ('filters.html#hessian', 'diptorch/filters.py'),
+                                  'diptorch.filters.hessian_eigenvalues': ('filters.html#hessian_eigenvalues', 'diptorch/filters.py')},
             'diptorch.linalg': { 'diptorch.linalg._is_hermitian': ('linalg.html#_is_hermitian', 'diptorch/linalg.py'),
                                  'diptorch.linalg._is_square': ('linalg.html#_is_square', 'diptorch/linalg.py'),
                                  'diptorch.linalg.deth3': ('linalg.html#deth3', 'diptorch/linalg.py'),

diptorch/filters.py

@@ -1,7 +1,7 @@
 # AUTOGENERATED! DO NOT EDIT! File to edit: ../notebooks/00_filters.ipynb.
 
 # %% auto 0
-__all__ = ['gaussian_filter', 'hessian']
+__all__ = ['gaussian_filter', 'hessian', 'hessian_eigenvalues']
 
 # %% ../notebooks/00_filters.ipynb 3
 from math import ceil, sqrt
@@ -94,6 +94,9 @@ def _conv(
     return convfn(x, weight=kernel.view(*view))
 
 # %% ../notebooks/00_filters.ipynb 8
+from .linalg import eigvalsh2, eigvalsh3
+
+
 def hessian(
     img: torch.Tensor, sigma: float, as_matrix: bool = False, **kwargs
 ) -> torch.Tensor:
@@ -110,6 +113,16 @@ def hessian(
     else:
         return hessian
 
+
+def hessian_eigenvalues(img: torch.Tensor, sigma: float, **kwargs):
+    H = hessian(img, sigma, **kwargs)
+    if len(H) == 3:
+        return eigvalsh2(*H)
+    elif len(H) == 6:
+        return eigvalsh3(*H)
+    else:
+        raise ValueError(f"Unrecognized number of upper triangular elements: {len(H)}")
+
 # %% ../notebooks/00_filters.ipynb 9
 def _hessian_2d(img: torch.Tensor, sigma: float, **kwargs):
     xx = gaussian_filter(img, sigma, order=[0, 2], **kwargs)

notebooks/00_filters.ipynb

@@ -152,7 +152,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Hessian matrix of an image"
+    "## Hessian matrix of an image\n",
+    "\n",
+    "Compute a symmetric matrix of all second-order partial Gaussian derivatives of an image."
    ]
   },
   {
@@ -162,7 +164,12 @@
    "outputs": [],
    "source": [
     "#| export\n",
-    "def hessian(img: torch.Tensor, sigma: float, as_matrix: bool = False, **kwargs) -> torch.Tensor:\n",
+    "from diptorch.linalg import eigvalsh2, eigvalsh3\n",
+    "\n",
+    "\n",
+    "def hessian(\n",
+    "    img: torch.Tensor, sigma: float, as_matrix: bool = False, **kwargs\n",
+    ") -> torch.Tensor:\n",
     "    \"\"\"Compute the Hessian of a 2D or 3D image.\"\"\"\n",
     "    if img.ndim == 4:\n",
     "        hessian = _hessian_2d(img, sigma, **kwargs)\n",
@@ -174,7 +181,17 @@
     "    if as_matrix:\n",
     "        return _hessian_as_matrix(*hessian)\n",
     "    else:\n",
-    "        return hessian"
+    "        return hessian\n",
+    "\n",
+    "\n",
+    "def hessian_eigenvalues(img: torch.Tensor, sigma: float, **kwargs):\n",
+    "    H = hessian(img, sigma, **kwargs)\n",
+    "    if len(H) == 3:\n",
+    "        return eigvalsh2(*H)\n",
+    "    elif len(H) == 6:\n",
+    "        return eigvalsh3(*H)\n",
+    "    else:\n",
+    "        raise ValueError(f\"Unrecognized number of upper triangular elements: {len(H)}\")"
    ]
   },
   {