Reduce ridge filters memory footprints (#6509)

* Reduce ridge filters memory footprints by not storing intermediate result and taking max between iterations.

* Clean PEP8 whitespace issue.

* Add memory benchmarks for ridge filters

* Reformat indent in math expression

Co-authored-by: Riadh Fezzani <rfezzani@gmail.com>

* Add CPU benchmarks for ridge filters

* Remake memory footprint fix.

* Use zeros_like where possible

Co-authored-by: Riadh Fezzani <rfezzani@gmail.com>

Co-authored-by: Lars Grüter <lagru@mailbox.org>
Co-authored-by: Lars Grüter <lagru+github@mailbox.org>
Co-authored-by: Riadh Fezzani <rfezzani@gmail.com>

benchmarks/benchmark_filters.py

@@ -2,7 +2,7 @@
 # https://asv.readthedocs.io/en/latest/writing_benchmarks.html
 import numpy as np
 
-from skimage import data, filters
+from skimage import data, filters, color
 from skimage.filters.thresholding import threshold_li
 
 
@@ -103,3 +103,50 @@ def time_integer_image(self):
 
     def time_float32_image(self):
         result1 = threshold_li(self.image_float32)
+
+
+class RidgeFilters:
+    """Benchmark ridge filters in scikit-image."""
+
+    def setup(self):
+        # Ensure memory footprint of lazy import is included in reference
+        self._ = filters.meijering, filters.sato, filters.frangi, filters.hessian
+        self.image = color.rgb2gray(data.retina())
+
+    def peakmem_setup(self):
+        """peakmem includes the memory used by setup.
+        Peakmem benchmarks measure the maximum amount of RAM used by a
+        function. However, this maximum also includes the memory used
+        by ``setup`` (as of asv 0.2.1; see [1]_)
+        Measuring an empty peakmem function might allow us to disambiguate
+        between the memory used by setup and the memory used by slic (see
+        ``peakmem_slic_basic``, below).
+        References
+        ----------
+        .. [1]: https://asv.readthedocs.io/en/stable/writing_benchmarks.html#peak-memory
+        """
+        pass
+
+    def time_meijering(self):
+        filters.meijering(self.image)
+
+    def peakmem_meijering(self):
+        filters.meijering(self.image)
+
+    def time_sato(self):
+        filters.sato(self.image)
+
+    def peakmem_sato(self):
+        filters.sato(self.image)
+
+    def time_frangi(self):
+        filters.frangi(self.image)
+
+    def peakmem_frangi(self):
+        filters.frangi(self.image)
+
+    def time_hessian(self):
+        filters.hessian(self.image)
+
+    def peakmem_hessian(self):
+        filters.hessian(self.image)

skimage/filters/ridges.py

@@ -149,7 +149,9 @@ def meijering(image, sigmas=range(1, 10, 2), alpha=None,
     mtx = linalg.circulant(
         [1, *[alpha] * (image.ndim - 1)]).astype(image.dtype)
 
-    filtered = []
+    # Generate empty array for storing maximum value
+    # from different (sigma) scales
+    filtered_max = np.zeros_like(image)
     for sigma in sigmas:  # Filter for all sigmas.
         eigvals = hessian_matrix_eigvals(hessian_matrix(
             image, sigma, mode=mode, cval=cval, use_gaussian_derivatives=True))
@@ -164,8 +166,9 @@ def meijering(image, sigmas=range(1, 10, 2), alpha=None,
         max_val = vals.max()
         if max_val > 0:
             vals /= max_val
-        filtered.append(vals)
-    return np.max(filtered, axis=0)  # Return pixel-wise max over all sigmas.
+        filtered_max = np.maximum(filtered_max, vals)
+
+    return filtered_max  # Return pixel-wise max over all sigmas.
 
 
 def sato(image, sigmas=range(1, 10, 2), black_ridges=True,
@@ -221,7 +224,9 @@ def sato(image, sigmas=range(1, 10, 2), black_ridges=True,
     if not black_ridges:  # Normalize to black ridges.
         image = -image
 
-    filtered = []
+    # Generate empty array for storing maximum value
+    # from different (sigma) scales
+    filtered_max = np.zeros_like(image)
     for sigma in sigmas:  # Filter for all sigmas.
         eigvals = hessian_matrix_eigvals(hessian_matrix(
             image, sigma, mode=mode, cval=cval, use_gaussian_derivatives=True))
@@ -232,8 +237,8 @@ def sato(image, sigmas=range(1, 10, 2), black_ridges=True,
         eigvals = eigvals[:-1]
         vals = (sigma ** 2
                 * np.prod(np.maximum(eigvals, 0), 0) ** (1 / len(eigvals)))
-        filtered.append(vals)
-    return np.max(filtered, axis=0)  # Return pixel-wise max over all sigmas.
+        filtered_max = np.maximum(filtered_max, vals)
+    return filtered_max  # Return pixel-wise max over all sigmas.
 
 
 def frangi(image, sigmas=range(1, 10, 2), scale_range=None,
@@ -318,7 +323,9 @@ def frangi(image, sigmas=range(1, 10, 2), scale_range=None,
     if not black_ridges:  # Normalize to black ridges.
         image = -image
 
-    filtered = []
+    # Generate empty array for storing maximum value
+    # from different (sigma) scales
+    filtered_max = np.zeros_like(image)
     for sigma in sigmas:  # Filter for all sigmas.
         eigvals = hessian_matrix_eigvals(hessian_matrix(
             image, sigma, mode=mode, cval=cval, use_gaussian_derivatives=True))
@@ -345,8 +352,8 @@ def frangi(image, sigmas=range(1, 10, 2), scale_range=None,
         vals = 1.0 - np.exp(-r_a**2 / (2 * alpha**2))  # plate sensitivity
         vals *= np.exp(-r_b**2 / (2 * beta**2))  # blobness
         vals *= 1.0 - np.exp(-s**2 / (2 * gamma**2))  # structuredness
-        filtered.append(vals)
-    return np.max(filtered, axis=0)  # Return pixel-wise max over all sigmas.
+        filtered_max = np.maximum(filtered_max, vals)
+    return filtered_max  # Return pixel-wise max over all sigmas.
 
 
 def hessian(image, sigmas=range(1, 10, 2), scale_range=None, scale_step=None,