Pull in LPI filters and add image reading (via PIL for now).

scikits/image/filter/__init__.py

@@ -0,0 +1,2 @@
+from lpi_filter import *
+

scikits/image/filter/lpi_filter.py

@@ -0,0 +1,158 @@
+"""
+:author: Stefan van der Walt, 2008
+:license: modified BSD
+"""
+
+__all__ = ['LPIFilter2D']
+__docformat__ = 'restructuredtext en'
+
+import numpy as np
+from scipy.fftpack import fftshift, ifftshift
+
+eps = np.finfo(float).eps
+
+class LPIFilter2D(object):
+    """Linear Position-Invariant Filter (2-dimensional)
+
+    """
+    def __init__(self,impulse_response,**filter_params):
+        """
+        *Parameters*:
+            impulse_response : callable f(r,c,**filter_params)
+                Function that yields the impulse response.  `r` and
+                `c` are 1-dimensional vectors that represent row and
+                column positions, in other words coordinates are
+                (r[0],c[0]),(r[0],c[1]) etc.  `**filter_params` are
+                passed through.
+
+                In other words, example would be called like this:
+
+                r = [0,0,0,1,1,1,2,2,2]
+                c = [0,1,2,0,1,2,0,1,2]
+                impulse_response(r,c,**filter_params)
+
+        *Example*:
+
+           Gaussian filter:
+
+           >>> def filt_func(r,c):
+                   return np.exp(-np.hypot(r,c)/1)
+
+           >>> filter = LPIFilter2D(filt_func)
+
+
+        """
+        self.impulse_response = impulse_response
+        self.filter_params = filter_params
+        self._cache = None
+
+    def _pad(self,data,shape):
+        """Pad the data to the given shape with zeros.
+
+        *Parameters*:
+            data : 2-d ndarray
+                Input data
+            shape : (2,) tuple
+
+        """
+        out = np.zeros(shape)
+        out[[slice(0,n) for n in data.shape]] = data
+        return out
+
+    def _prepare(self,data):
+        """Calculate filter and data FFT in preparation for filtering.
+
+        """
+        dshape = np.array(data.shape)
+        dshape += (dshape %2 == 0) # all filter dimensions must be uneven
+        oshape = np.array(data.shape)*2-1
+
+        if self._cache is None or np.any(self._cache.shape != oshape):
+            coords = np.mgrid[[slice(0,float(n)) for n in dshape]]
+            # this steps over two sets of coordinates,
+            # not over the coordinates individually
+            for k,coord in enumerate(coords):
+                coord -= (dshape[k]-1)/2.
+            coords = coords.reshape(2,-1).T # coordinate pairs (r,c)
+
+            f = self.impulse_response(coords[:,0],coords[:,1],
+                                      **self.filter_params).reshape(dshape)
+
+            f = self._pad(f,oshape)
+            F = np.dual.fftn(f)
+            self._cache = F
+        else:
+            F = self._cache
+
+        data = self._pad(data,oshape)
+        G = np.dual.fftn(data)
+
+        return F,G
+
+    def _min_limit(self,x,val=eps):
+        mask = np.abs(x) < eps
+        x[mask] = np.sign(x[mask])*eps
+
+    def _centre(self,x,oshape):
+        """Return an array of oshape from the centre of x.
+
+        """
+        start = (np.array(x.shape) - np.array(oshape))/2.+1
+        out = x[[slice(s,s+n) for s,n in zip(start,oshape)]]
+        return out
+
+    def __call__(self,data):
+        """Apply the filter to the given data.
+
+        *Parameters*:
+            data : (M,N) ndarray
+
+        """
+        F,G = self._prepare(data)
+        out = np.dual.ifftn(F*G)
+        out = np.abs(self._centre(out,data.shape))
+        return out
+
+    def inverse(self,data,max_gain=2):
+        """Apply the filter in reverse to the given data.
+
+        *Parameters*:
+            data : (M,N) ndarray
+                Input data.
+            max_gain : float
+                Limit the filter gain.  Often, the filter contains
+                zeros, which would cause the inverse filter to have
+                infinite gain.  High gain causes amplification of
+                artefacts, so a conservative limit is recommended.
+
+        """
+        F,G = self._prepare(data)
+        self._min_limit(F)
+
+        F = 1/F
+        mask = np.abs(F) > max_gain
+        F[mask] = np.sign(F[mask])*max_gain
+
+        return self._centre(np.abs(ifftshift(np.dual.ifftn(G*F))),data.shape)
+
+    def wiener(self,data,K=0.25):
+        """Minimum Mean Square Error (Wiener) inverse filter.
+
+        *Parameters*:
+            data : (M,N) ndarray
+                Input data.
+            K : float or (M,N) ndarray
+                Ratio between power spectrum of noise and undegraded
+                image.
+
+        """
+        F,G = self._prepare(data)
+        self._min_limit(F)
+
+        H_mag_sqr = np.abs(F)**2
+        F = 1/F * H_mag_sqr / (H_mag_sqr + K)
+
+        return self._centre(np.abs(ifftshift(np.dual.ifftn(G*F))),data.shape)
+
+    def constrained_least_squares(self,data,lam):
+        pass

scikits/image/filter/tests/test_lpi_filter.py

@@ -0,0 +1,56 @@
+import os.path
+
+import numpy as np
+from numpy.testing import *
+
+from scikits.image.io import *
+from scikits.image.filter import *
+
+data_dir = os.path.join(os.path.dirname(__file__), './data/')
+
+class TestLPIFilter2D():
+    img = imread(os.path.join(data_dir + 'camera.png'),
+                 flatten=True)[:-101,:-100]
+
+    def filt_func(self,r,c):
+        return np.exp(-np.hypot(r,c)/1)
+
+    def setUp(self):
+        self.f = LPIFilter2D(self.filt_func)
+
+    def tst_shape(self, x):
+        X = self.f(x)
+        assert_equal(X.shape,x.shape)
+
+    def test_ip_shape(self):
+        rows,columns = self.img.shape[:2]
+
+        for c_slice in [slice(0,columns),slice(0,columns-5),
+                        slice(0,columns-100)]:
+            yield (self.tst_shape,self.img[:,c_slice])
+
+    def test_inverse(self):
+        F = self.f(self.img)
+        g = self.f.inverse(F)
+        assert_equal(g.shape,self.img.shape)
+
+        g1 = self.f.inverse(F[::-1,::-1])
+        assert ((g-g1[::-1,::-1]).sum() < 55)
+
+        # test cache
+        g1 = self.f.inverse(F[::-1,::-1])
+        assert ((g-g1[::-1,::-1]).sum() < 55)
+
+
+    def test_wiener(self):
+        F = self.f(self.img)
+        g = self.f.wiener(F)
+        assert_equal(g.shape,self.img.shape)
+
+        g1 = self.f.wiener(F[::-1,::-1])
+        assert ((g-g1[::-1,::-1]).sum() < 1)
+
+
+if __name__ == "__main__":
+    run_module_suite()
+

scikits/image/io/__init__.py

@@ -0,0 +1 @@
+from pil_imread import *

scikits/image/io/pil_imread.py

@@ -0,0 +1,34 @@
+__all__ = ['imread']
+
+import numpy as np
+
+def imread(fname, flatten=False):
+    """Load an image from file.
+
+    Parameters
+    ----------
+    fname : string
+        Image file name, e.g. ``test.jpg``.
+    flatten : bool
+        If true, convert the output to grey-scale.
+
+    Returns
+    -------
+    img_array : ndarray
+        The different colour bands/channels are stored in the
+        third dimension, such that a grey-image is MxN, an
+        RGB-image MxNx3 and an RGBA-image MxNx4.
+
+    """
+    try:
+        from PIL import Image
+    except ImportError:
+        raise ImportError("Could not import the Python Imaging Library (PIL)"
+                          " required to load image files.  Please refer to"
+                          " http://pypi.python.org/pypi/PIL/ for installation"
+                          " instructions.")
+
+    im = Image.open(fname)
+    if flatten:
+        im = im.convert('F')
+    return np.array(im)