Fix LPI filter for data with even dimensions

skimage/filters/__init__.pyi

@@ -13,6 +13,7 @@ __all__ = [
     "farid_h",
     "farid_v",
     "filter_inverse",
+    "filter_forward",
     "frangi",
     "gabor",
     "gabor_kernel",
@@ -84,6 +85,7 @@ from .edges import (
 from .lpi_filter import (
     LPIFilter2D,
     filter_inverse,
+    filter_forward,
     inverse,
     wiener,
 )

skimage/filters/lpi_filter.py

@@ -18,7 +18,7 @@ def _center(x, oshape):
     """Return an array of shape ``oshape`` from the center of array ``x``.
 
     """
-    start = (np.array(x.shape) - np.array(oshape)) // 2 + 1
+    start = (np.array(x.shape) - np.array(oshape)) // 2
     out = x[tuple(slice(s, s + n) for s, n in zip(start, oshape))]
     return out
 
@@ -66,11 +66,10 @@ def __init__(self, impulse_response, **filter_params):
 
         Examples
         --------
-        Gaussian filter: Use a 1-D gaussian in each direction without
-        normalization coefficients.
+        Gaussian filter without normalization of coefficients:
 
-        >>> def filt_func(r, c, sigma = 1):
-        ...     return np.exp(-np.hypot(r, c)/sigma)
+        >>> def filt_func(r, c, sigma=1):
+        ...     return np.exp(-(r**2 + c**2)/(2 * sigma**2))
         >>> filter = LPIFilter2D(filt_func)
 
         """
@@ -86,14 +85,18 @@ def _prepare(self, data):
 
         """
         dshape = np.array(data.shape)
-        dshape += (dshape % 2 == 0)  # all filter dimensions must be uneven
+        even_offset = (dshape % 2 == 0).astype(int)
+        dshape += even_offset  # all filter dimensions must be uneven
         oshape = np.array(data.shape) * 2 - 1
 
         float_dtype = _supported_float_type(data.dtype)
         data = data.astype(float_dtype, copy=False)
 
         if self._cache is None or np.any(self._cache.shape != oshape):
-            coords = np.mgrid[[slice(0, float(n)) for n in dshape]]
+            coords = np.mgrid[
+                [slice(0 + offset, float(n + offset))
+                 for (n, offset) in zip(dshape, even_offset)]
+            ]
             # this steps over two sets of coordinates,
             # not over the coordinates individually
             for k, coord in enumerate(coords):
@@ -152,10 +155,10 @@ def filter_forward(data, impulse_response=None, filter_params=None,
     Examples
     --------
 
-    Gaussian filter:
+    Gaussian filter without normalization:
 
-    >>> def filt_func(r, c):
-    ...     return np.exp(-np.hypot(r, c)/1)
+    >>> def filt_func(r, c, sigma=1):
+    ...     return np.exp(-(r**2 + c**2)/(2 * sigma**2))
     >>>
     >>> from skimage import data
     >>> filtered = filter_forward(data.coins(), filt_func)

skimage/filters/tests/test_lpi_filter.py

@@ -1,18 +1,46 @@
 import numpy as np
 import pytest
-from numpy.testing import assert_, assert_equal
+from numpy.testing import assert_, assert_equal, assert_array_almost_equal
 
 from skimage._shared.utils import _supported_float_type
-from skimage.data import camera
-from skimage.filters.lpi_filter import LPIFilter2D, filter_inverse, wiener
+from skimage.data import camera, coins
+from skimage.filters import (
+    LPIFilter2D, filter_inverse, filter_forward, wiener, gaussian
+)
+
+
+def test_filter_forward():
+    def filt_func(r, c, sigma=2):
+        return (1 / (2 * np.pi * sigma**2)) *  np.exp(-(r**2 + c**2) / (2 * sigma ** 2))
+
+    gaussian_args = {
+        'sigma': 2,
+        'preserve_range': True,
+        'mode': 'constant',
+        'truncate': 20  # LPI filtering is more precise than the truncated
+                        # Gaussian, so don't truncate at the default of 4 sigma
+    }
+
+    # Odd image size
+    image = coins()[:303, :383]
+    filtered = filter_forward(image, filt_func)
+    filtered_gaussian = gaussian(image, **gaussian_args)
+    assert_array_almost_equal(filtered, filtered_gaussian)
+
+    # Even image size
+    image = coins()
+    filtered = filter_forward(image, filt_func)
+    filtered_gaussian = gaussian(image, **gaussian_args)
+    assert_array_almost_equal(filtered, filtered_gaussian)
+
 
 
 class TestLPIFilter2D:
 
     img = camera()[:50, :50]
 
     def filt_func(self, r, c):
-        return np.exp(-np.hypot(r, c) / 1)
+        return np.exp(-(r**2 + c**2) / 1)
 
     def setup_method(self):
         self.f = LPIFilter2D(self.filt_func)