diff --git a/scikits/image/filter/__init__.py b/scikits/image/filter/__init__.py new file mode 100644 index 00000000..4bdff030 --- /dev/null +++ b/scikits/image/filter/__init__.py @@ -0,0 +1,2 @@ +from lpi_filter import * + diff --git a/scikits/image/filter/lpi_filter.py b/scikits/image/filter/lpi_filter.py new file mode 100644 index 00000000..57564921 --- /dev/null +++ b/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 diff --git a/scikits/image/filter/tests/data/camera.png b/scikits/image/filter/tests/data/camera.png new file mode 100644 index 00000000..a196b835 Binary files /dev/null and b/scikits/image/filter/tests/data/camera.png differ diff --git a/scikits/image/filter/tests/test_lpi_filter.py b/scikits/image/filter/tests/test_lpi_filter.py new file mode 100644 index 00000000..6cb804f9 --- /dev/null +++ b/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() + diff --git a/scikits/image/io/__init__.py b/scikits/image/io/__init__.py new file mode 100644 index 00000000..0d5394a5 --- /dev/null +++ b/scikits/image/io/__init__.py @@ -0,0 +1 @@ +from pil_imread import * diff --git a/scikits/image/io/pil_imread.py b/scikits/image/io/pil_imread.py new file mode 100644 index 00000000..d57348a8 --- /dev/null +++ b/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)