diff --git a/doc/examples/plot_deconvolution.py b/doc/examples/plot_deconvolution.py new file mode 100644 index 00000000..fd94e77f --- /dev/null +++ b/doc/examples/plot_deconvolution.py @@ -0,0 +1,58 @@ +# -*- coding: utf-8 -*- +""" +===================== +Deconvolution of Lena +===================== + +In this example, we deconvolve a noisy version of Lena using wiener +and unsupervised wiener algorithms. This algorithms are known to not +have the best respect of sharp edge in the image. + +Wiener filter +------------- + +This is simply the inverse filter based on the PSF, the prior +regularisation (penalisation of high frequency) and the tradeoff +between the data and prior adequacy. The regularization parameter must +be hand tuned. + +Unsupervised wiener +------------------- + +This algorithm has a self tuned regularisation parameters based on +data learning. This is not common and based on the following publication + +.. [1] François Orieux, Jean-François Giovannelli, and Thomas + Rodet, "Bayesian estimation of regularization and point + spread function parameters for Wiener-Hunt deconvolution", + J. Opt. Soc. Am. A 27, 1593-1607 (2010) +""" +import numpy as np +import matplotlib.pyplot as plt + +from skimage import color, data, deconvolution + +lena = color.rgb2gray(data.lena()) +from scipy.signal import convolve2d as conv2 +psf = np.ones((5, 5)) +lena = conv2(lena, psf, 'same') +lena += 0.1 * lena.std() * np.random.standard_normal(lena.shape) + +deconvolued, _ = deconvolution.unsupervised_wiener(lena, psf) + +fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8, 5)) + +plt.gray() + +ax[0].imshow(lena) +ax[0].axis('off') +ax[0].set_title('Data') + +ax[1].imshow(deconvolued) +ax[1].axis('off') +ax[1].set_title('Deconvolution') + +fig.subplots_adjust(wspace=0.02, hspace=0.2, + top=0.9, bottom=0.05, left=0, right=1) + +plt.show() diff --git a/skimage/deconvolution/__init__.py b/skimage/deconvolution/__init__.py new file mode 100644 index 00000000..508f9d1c --- /dev/null +++ b/skimage/deconvolution/__init__.py @@ -0,0 +1,27 @@ +"""Skimage module for image deconvolution + +This module implement various algorithm of the literarture for image +deconvolution. + +References +---------- +.. [1] François Orieux, Jean-François Giovannelli, and Thomas + Rodet, "Bayesian estimation of regularization and point + spread function parameters for Wiener-Hunt deconvolution", + J. Opt. Soc. Am. A 27, 1593-1607 (2010) + + http://www.opticsinfobase.org/josaa/abstract.cfm?URI=josaa-27-7-1593 + +.. [2] Richardson, William Hadley, "Bayesian-Based Iterative Method of + Image Restoration". JOSA 62 (1): 55–59. doi:10.1364/JOSA.62.000055, 1972 + +.. [3] B. R. Hunt "A matrix theory proof of the discrete + convolution theorem", IEEE Trans. on Audio and + Electroacoustics, vol. au-19, no. 4, pp. 285-288, dec. 1971 +""" + +from .wiener import wiener, unsupervised_wiener, richardson_lucy + +__all__ = ['wiener', + "unsupervised_wiener", + "richardson_lucy"] diff --git a/skimage/deconvolution/test/camera_unsup.npy b/skimage/deconvolution/test/camera_unsup.npy new file mode 100644 index 00000000..a10fbd76 Binary files /dev/null and b/skimage/deconvolution/test/camera_unsup.npy differ diff --git a/skimage/deconvolution/test/camera_wiener.npy b/skimage/deconvolution/test/camera_wiener.npy new file mode 100644 index 00000000..f095ff8c Binary files /dev/null and b/skimage/deconvolution/test/camera_wiener.npy differ diff --git a/skimage/deconvolution/test/test_deconvolution.py b/skimage/deconvolution/test/test_deconvolution.py new file mode 100644 index 00000000..b542db6a --- /dev/null +++ b/skimage/deconvolution/test/test_deconvolution.py @@ -0,0 +1,38 @@ +import warnings + +import numpy as np +import numpy.testing.assert_array_almost_equal + +from scipy.signal import convolve2d as conv2 +from skimage import data, deconvolution + +# Test deconvolution +# =========================== + +test_img = data.camera().astype(np.float) + + +def test_wiener(): + psf = np.ones((5, 5)) + data = conv2(test_img, psf, 'same') + np.random.seed(0) + data += 0.1 * data.std() * np.random.standard_normal(data.shape) + deconvolued = deconvolution.wiener(data, psf, 25) + + numpy.testing.assert_array_almost_equal(deconvolued, + np.load("./camera_wiener.npy")) + + +def test_unsupervised_wiener(): + psf = np.ones((5, 5)) + data = conv2(test_img, psf, 'same') + np.random.seed(0) + data += 0.1 * data.std() * np.random.standard_normal(data.shape) + deconvolued, _ = deconvolution.unsupervised_wiener(data, psf) + + numpy.testing.assert_array_almost_equal(deconvolued, + np.load("./camera_unsup.npy")) + + +def test_rychardson_lucy(): + return True diff --git a/skimage/deconvolution/uft.py b/skimage/deconvolution/uft.py new file mode 100644 index 00000000..fedbcb3f --- /dev/null +++ b/skimage/deconvolution/uft.py @@ -0,0 +1,473 @@ +# -*- coding: utf-8 -*- +# uft.py --- Unitary fourier transform + +# Copyright (c) 2011, 2012, 2013 François Orieux + +# Permission is hereby granted, free of charge, to any person +# obtaining a copy of this software and associated documentation files +# (the "Software"), to deal in the Software without restriction, +# including without limitation the rights to use, copy, modify, merge, +# publish, distribute, sublicense, and/or sell copies of the Software, +# and to permit persons to whom the Software is furnished to do so, +# subject to the following conditions: + +# The above copyright notice and this permission notice shall be +# included in all copies or substantial portions of the Software. + +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS +# BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN +# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN +# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +# Commentary: + +"""Function of unitary fourier transform and utilities + +This module implement unitary fourier transform, that is ortho-normal +transform. They are specially usefull for convolution [1]: they +respect the parseval equality, the value of the null frequency is +equal to + +.. math:: \frac{1}{\sqrt{n}} \sum_i x_i. + +If the anfft module is present, his function are used. anfft wrap fftw +C library. Otherwise, numpy fft functions are used. + +You must keep in mind that the transform are applied from the last +axes. this is a fftw convention for performance reason (c order +array). If you want more sofisticated use, you must use directly the +numpy.fft module. + +References +---------- +.. [1] B. R. Hunt "A matrix theory proof of the discrete convolution + theorem", IEEE Trans. on Audio and Electroacoustics, + vol. au-19, no. 4, pp. 285-288, dec. 1971 + +""" + +# code: + +import logging + +import numpy as np +try: + import anfft + ANFFTMOD = True +except ImportError: + logging.info("Installation of the anfft package improve preformance" + " by using fftw library.") + ANFFTMOD = False + +__author__ = "François Orieux" +__copyright__ = "Copyright (C) 2011, 2012, 2013 F. Orieux " +__credits__ = ["François Orieux"] +__license__ = "mit" +__version__ = "0.1.0" +__maintainer__ = "François Orieux" +__email__ = "orieux@iap.fr" +__status__ = "development" +__url__ = "" +__keywords__ = "fft" + + +def _circshift(inarray, shifts): + """Shift array circularly. + + Circularly shifts the values in the array `a` by `s` + elements. Return a copy. + + Parameters + ---------- + a : ndarray + The array to shift. + + s : tuple of int + A tuple of integer scalars where the N-th element specifies the + shift amount for the N-th dimension of array `a`. If an element + is positive, the values of `a` are shifted down (or to the + right). If it is negative, the values of `a` are shifted up (or + to the left). + + Returns + ------- + y : ndarray + The shifted array (elements are copied) + + Examples + -------- + >>> circshift(np.arange(10), 2) + array([8, 9, 0, 1, 2, 3, 4, 5, 6, 7]) + + """ + # Initialize array of indices + idx = [] + + # Loop through each dimension of the input matrix to calculate + # shifted indices + for dim in range(inarray.ndim): + length = inarray.shape[dim] + try: + shift = shifts[dim] + except IndexError: + shift = 0 # no shift if not specify + + # Lets start for fancy indexing. First we build the shifted + # index for dim k. It will be broadcasted to other dim so + # ndmin is specified + index = np.mod(np.array(range(length), + ndmin=inarray.ndim) - shift, + length) + # Shape adaptation + shape = np.ones(inarray.ndim) + shape[dim] = inarray.shape[dim] + index = np.reshape(index, shape) + + idx.append(index.astype(int)) + + # Perform the actual conversion by indexing into the input matrix + return inarray[idx] + + +def ufftn(inarray, dim=None): + """N-dim unitary Fourier transform + + Parameters + ---------- + inarray : ndarray + The array to transform. + + dim : int, optional + The `dim` last axis along wich to compute the transform. All + axes by default. + + Returns + ------- + outarray : array-like (same shape than inarray) + """ + if not dim: + dim = inarray.ndim + + if ANFFTMOD: + outarray = anfft.fftn(inarray, k=dim) + else: + outarray = np.fft.fftn(inarray, axes=range(-dim, 0)) + + return outarray / np.sqrt(np.prod(inarray.shape[-dim:])) + + +def uifftn(inarray, dim=None): + """N-dim unitary inverse Fourier transform + + Parameters + ---------- + inarray : ndarray + The array to transform. + + dim : int, optional + The `dim` last axis along wich to compute the transform. All + axes by default. + + Returns + ------- + outarray : array-like (same shape than inarray) + """ + if not dim: + dim = inarray.ndim + + if ANFFTMOD: + outarray = anfft.ifftn(inarray, k=dim) + else: + outarray = np.fft.ifftn(inarray, axes=range(-dim, 0)) + + return outarray * np.sqrt(np.prod(inarray.shape[-dim:])) + + +def urfftn(inarray, dim=None): + """N-dim real unitary Fourier transform + + This transform consider the Hermitian property of the transform on + real input + + Parameters + ---------- + inarray : ndarray + The array to transform. + + dim : int, optional + The `dim` last axis along wich to compute the transform. All + axes by default. + + Returns + ------- + outarray : array-like (the last dim as N / 2 + 1 lenght) + """ + if not dim: + dim = inarray.ndim + + if ANFFTMOD: + outarray = anfft.rfftn(inarray, k=dim) + else: + outarray = np.fft.rfftn(inarray, axes=range(-dim, 0)) + + return outarray / np.sqrt(np.prod(inarray.shape[-dim:])) + + +def uirfftn(inarray, dim=None): + """N-dim real unitary Fourier transform + + This transform consider the Hermitian property of the transform + from complex to real real input. + + Parameters + ---------- + inarray : ndarray + The array to transform. + + dim : int, optional + The `dim` last axis along wich to compute the transform. All + axes by default. + + Returns + ------- + outarray : array-like (the last dim as (N - 1) *2 lenght) + """ + if not dim: + dim = inarray.ndim + + if ANFFTMOD: + outarray = anfft.irfftn(inarray, k=dim) + else: + outarray = np.fft.irfftn(inarray, axes=range(-dim, 0)) + + return outarray * np.sqrt(np.prod(inarray.shape[-dim:-1]) * + (inarray.shape[-1] - 1) * 2) + + +def ufft2(inarray): + """2-dim unitary Fourier transform + + Compute the Fourier transform on the last 2 axes. + + Parameters + ---------- + inarray : ndarray + The array to transform. + + Returns + ------- + outarray : array-like (same shape than inarray) + + See Also + -------- + uifft2, ufftn, urfftn + """ + return ufftn(inarray, 2) + + +def uifft2(inarray): + """2-dim inverse unitary Fourier transform + + Compute the inverse Fourier transform on the last 2 axes. + + Parameters + ---------- + inarray : ndarray + The array to transform. + + Returns + ------- + outarray : array-like (same shape than inarray) + + See Also + -------- + uifft2, uifftn, uirfftn + """ + return uifftn(inarray, 2) + + +def urfft2(inarray): + """2-dim real unitary Fourier transform + + Compute the real Fourier transform on the last 2 axes. This + transform consider the Hermitian property of the transform from + complex to real real input. + + Parameters + ---------- + inarray : ndarray + The array to transform. + + Returns + ------- + outarray : array-like (the last dim as (N - 1) *2 lenght) + + See Also + -------- + ufft2, ufftn, urfftn + """ + return urfftn(inarray, 2) + + +def uirfft2(inarray): + """2-dim real unitary Fourier transform + + Compute the real inverse Fourier transform on the last 2 axes. + This transform consider the Hermitian property of the transform + from complex to real real input. + + Parameters + ---------- + inarray : ndarray + The array to transform. + + Returns + ------- + outarray : array-like (the last dim as (N - 1) *2 lenght) + + See Also + -------- + urfft2, uifftn, uirfftn + """ + return uirfftn(inarray, 2) + + +def image_quad_norm(inarray): + """Return quadratic norm of images in Fourier space + + This function detect if the image suppose the hermitian property. + + Parameters + ---------- + inarray : array-like + The images are supposed to be in the last two axes + + Returns + ------- + norm : float + + """ + # If there is an hermitian symmetry + if inarray.shape[-1] != inarray.shape[-2]: + return 2 * np.sum(np.sum(np.abs(inarray)**2, axis=-1), axis=-1) - \ + np.sum(np.abs(inarray[..., 0])**2, axis=-1) + else: + return np.sum(np.sum(np.abs(inarray)**2, axis=-1), axis=-1) + + +def crandn(shape): + """white complex gaussian noise + + Generate directly the unitary Fourier transform of white gaussian + noise noise field (with given shape) of zero mean and variance + unity (ie N(0,1)). + """ + return np.sqrt(0.5) * (np.random.standard_normal(shape) + + 1j * np.random.standard_normal(shape)) + + +def ir2tf(imp_resp, shape, dim=None, real=True): + """Compute the transfer function of IR + + This function make the necessary correct zero-padding, zero + convention, correct fft2 etc... to compute the transfer function + of IR. To use with unitary Fourier transform for the signal (ufftn + or equivalent). + + Parameters + ---------- + imp_resp : ndarray + The impulsionnal responses. + + shape : tuple of int + A tuple of integer corresponding to the target shape of the + tranfert function. + + dim : int, optional + The `dim` last axis along wich to compute the transform. All + axes by default. + + real : boolean (optionnal, default True) + If True, imp_resp is supposed real and the hermissian property + is used with rfftn Fourier transform. + + Returns + ------- + y : complex ndarray + The tranfert function of shape `shape`. + + See Also + -------- + ufftn, uifftn, urfftn, uirfftn + + Notes + ----- + The input array can be composed of multiple dimentionnal IR with + an arbitraru number of IR. The individual IR must be accesed + through first axes. The last `dim` axes of space definition. The + `dim` parameter must be specified to compute the transform only + along these last axes. + """ + if not dim: + dim = imp_resp.ndim + + # Zero padding and fill + irpadded = np.zeros(shape) + irpadded[tuple([slice(0, s) for s in imp_resp.shape])] = imp_resp + # Circshift fo zero convention of the fft to avoid the phase + # problem. Work with odd and even size. + irpadded = _circshift(irpadded, + [-np.floor(s / 2) + if i >= imp_resp.ndim - dim + else 0 + for i, s in enumerate(imp_resp.shape)]) + + if real: + if anfft: + return anfft.rfftn(irpadded, k=dim) + else: + return np.fft.rfftn(irpadded, axes=range(-dim, 0)) + else: + if anfft: + return anfft.fftn(irpadded, k=dim) + else: + return np.fft.fftn(irpadded, axes=range(-dim, 0)) + + +def laplacian(ndim, shape): + """Return the transfert function of the laplacian + + Laplacian is the second order difference, on line and column. + + Parameters + ---------- + ndim : int + The dimension of the laplacian + + shape : tuple, shape + The support on which to compute the transfert function + + Returns + ------- + tf : array_like, complex + The transfert function + + impr : array_like, real + The laplacian + """ + impr = np.zeros([3] * ndim) + for dim in range(ndim): + idx = tuple([slice(1, 2)] * dim + + [slice(None)] + + [slice(1, 2)] * (ndim - dim - 1)) + impr[idx] = np.array([-1.0, + 0.0, + -1.0]).reshape([-1 if i == dim else 1 + for i in range(ndim)]) + impr[([slice(1, 2)] * ndim)] = 2.0 * ndim + + return ir2tf(impr, shape), impr diff --git a/skimage/deconvolution/wiener.py b/skimage/deconvolution/wiener.py new file mode 100644 index 00000000..7914679a --- /dev/null +++ b/skimage/deconvolution/wiener.py @@ -0,0 +1,325 @@ +# -*- coding: utf-8 -*- + +# Copyright (c) 2013 François Orieux + +# Permission is hereby granted, free of charge, to any person +# obtaining a copy of this software and associated documentation files +# (the "Software"), to deal in the Software without restriction, +# including without limitation the rights to use, copy, modify, merge, +# publish, distribute, sublicense, and/or sell copies of the Software, +# and to permit persons to whom the Software is furnished to do so, +# subject to the following conditions: + +# The above copyright notice and this permission notice shall be +# included in all copies or substantial portions of the Software. + +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS +# BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN +# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN +# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +# Commentary: + +"""Implementations deconvolution functions""" + +# code: + +from __future__ import division + +import numpy as np +import numpy.random as npr +from scipy.signal import convolve2d as conv2 + +import uft + +__author__ = "François Orieux" +__copyright__ = "Copyright (C) 2013 F. Orieux " +__credits__ = ["François Orieux"] +__license__ = "mit" +__version__ = "0.1.0" +__maintainer__ = "François Orieux" +__email__ = "orieux@iap.fr" +__status__ = "stable" +__url__ = "http://research.orieux.fr" +__keywords__ = "deconvolution, image" + + +def wiener(data, psf, reg_val, reg=None, real=True): + """Wiener-Hunt deconvolution + + return the deconvolution with a wiener-hunt approach (ie with + Fourier diagonalisation). + + Parameters + ---------- + data : (M, N) ndarray + The data + + psf : ndarray + The impulsionnal response in real space or the transfer + function. Differentiation is done with the dtype where + transfer function is supposed complex. + + reg_val : float + The regularisation parameter value. + + reg : ndarray, optional + The regularisation operator. The laplacian by + default. Otherwise, the same constraints that for `psf` + apply. + + real : boolean, optional + True by default. Specify if `psf` or `reg` are provided + with hermitian hypothesis or not. See uft module. + + Returns + ------- + im_deconv : (M, N) ndarray + The deconvolued data + + Exemples + -------- + >>> from skimage import color, data, deconvolution + >>> lena = color.rgb2gray(data.lena()) + >>> from scipy.signal import convolve2d as conv2 + >>> psf = np.ones((5, 5)) + >>> lena = conv2(lena, psf, 'same') + >>> lena += 0.1 * lena.std() * np.random.standard_normal(lena.shape) + >>> deconvolued_lena = deconvolution.wiener(lena, psf, 1100) + + References + ---------- + .. [1] François Orieux, Jean-François Giovannelli, and Thomas + Rodet, "Bayesian estimation of regularization and point + spread function parameters for Wiener-Hunt deconvolution", + J. Opt. Soc. Am. A 27, 1593-1607 (2010) + + http://www.opticsinfobase.org/josaa/abstract.cfm?URI=josaa-27-7-1593 + + .. [2] B. R. Hunt "A matrix theory proof of the discrete + convolution theorem", IEEE Trans. on Audio and + Electroacoustics, vol. au-19, no. 4, pp. 285-288, dec. 1971 + """ + if not reg: + reg, _ = uft.laplacian(data.ndim, data.shape) + if reg.dtype != np.complex: + reg = uft.ir2tf(reg, data.shape) + + if psf.shape != reg.shape: + trans_func = uft.ir2tf(psf, data.shape) + else: + trans_func = psf + + wiener_filter = np.conj(trans_func) / (np.abs(trans_func)**2 + + reg_val * np.abs(reg)**2) + if real: + return uft.uirfft2(wiener_filter * uft.urfft2(data)) + else: + return uft.uifft2(wiener_filter * uft.ufft2(data)) + + +def unsupervised_wiener(data, psf, reg=None, user_params=None): + """Unsupervised Wiener-Hunt deconvolution + + return the deconvolution with a wiener-hunt approach, where the + hyperparameters are estimated (or automatically tuned from a + practical point of view). The algorithm is a stochastic iterative + process (Gibbs sampler). + + If you use this work, please add a citation to the reference below. + + Parameters + ---------- + image : (M, N) ndarray + The data + + psf : ndarray + The impulsionnal response in real space or the transfer + function. Differentiation is done with the dtype where + transfer function is supposed complex. + + reg : ndarray, optional + The regularisation operator. The laplacian by + default. Otherwise, the same constraints that for `psf` + apply + + user_params : dict + dictionary of gibbs parameters. See below. + + Returns + ------- + x_postmean : (M, N) ndarray + The deconvolued data (the posterior mean) + + chains : dict + The keys 'noise' and 'prior' contains the chain list of noise and + prior precision respectively + + Other parameters + ---------------- + The key of user_params are + + threshold : float + The stopping criterion: the norm of the difference between to + successive approximated solution (empirical mean of object + sample). 1e-4 by default. + + burnin : int + The number of sample to ignore to start computation of the + mean. 100 by default. + + min_iter : int + The minimum number of iteration. 30 by default. + + max_iter : int + The maximum number of iteration if `threshold` is not + satisfied. 150 by default. + + callback : None + A user provided callable to which is passed, if the function + exists, the current image sample. This function can be used to + store the sample, or compute other moments than the mean. + + Exemples + -------- + >>> from skimage import color, data, deconvolution + >>> lena = color.rgb2gray(data.lena()) + >>> from scipy.signal import convolve2d as conv2 + >>> psf = np.ones((5, 5)) + >>> lena = conv2(lena, psf, 'same') + >>> lena += 0.1 * lena.std() * np.random.standard_normal(lena.shape) + >>> deconvolued_lena = deconvolution.unsupervised_wiener(lena, psf) + + References + ---------- + .. [1] François Orieux, Jean-François Giovannelli, and Thomas + Rodet, "Bayesian estimation of regularization and point + spread function parameters for Wiener-Hunt deconvolution", + J. Opt. Soc. Am. A 27, 1593-1607 (2010) + + http://www.opticsinfobase.org/josaa/abstract.cfm?URI=josaa-27-7-1593 + """ + params = {'threshold': 1e-4, 'max_iter': 200, + 'min_iter': 30, 'burnin': 15, 'callback': None} + params.update(user_params if user_params else {}) + + if not reg: + reg, _ = uft.laplacian(data.ndim, data.shape) + if reg.dtype != np.complex: + reg = uft.ir2tf(reg, data.shape) + + if psf.shape != reg.shape: + trans_fct = uft.ir2tf(psf, data.shape) + else: + trans_fct = psf + + # The mean of the object + x_postmean = np.zeros(trans_fct.shape) + # The previous computed mean in the iterative loop + prev_x_postmean = np.zeros(trans_fct.shape) + + # Difference between two successive mean + delta = np.NAN + + # Initial state of the chain + gn_chain, gx_chain = [1], [1] + + # The correlation of the object in Fourier space (if size is big, + # this can reduce computation time in the loop) + areg2 = np.abs(reg)**2 + atf2 = np.abs(trans_fct)**2 + + data = uft.urfft2(data.astype(np.float)) + + # Gibbs sampling + for iteration in range(params['max_iter']): + # Sample of Eq. 27 p(circX^k | gn^k-1, gx^k-1, y). + + # weighing (correlation in direct space) + precision = gn_chain[-1] * atf2 + gx_chain[-1] * areg2 # Eq. 29 + excursion = uft.crandn(data.shape) / np.sqrt(precision) + + # mean Eq. 30 (RLS for fixed gn, gamma0 and gamma1 ...) + wiener_filter = gn_chain[-1] * np.conj(trans_fct) / precision + + # sample of X in Fourier space + x_sample = wiener_filter * data + excursion + if params['callback']: + params['callback'](x_sample) + + # sample of Eq. 31 p(gn | x^k, gx^k, y) + gn_chain.append(npr.gamma(data.size / 2, + 2 / uft.image_quad_norm(data - x_sample * + trans_fct))) + + # sample of Eq. 31 p(gx | x^k, gn^k-1, y) + gx_chain.append(npr.gamma((data.size - 1) / 2, + 2 / uft.image_quad_norm(x_sample * reg))) + + # current empirical average + if iteration > params['burnin']: + x_postmean = prev_x_postmean + x_sample + + if iteration > (params['burnin'] + 1): + current = x_postmean / (iteration - params['burnin']) + previous = prev_x_postmean / (iteration - params['burnin'] - 1) + + delta = np.sum(np.abs(current - previous)) / \ + np.sum(np.abs(x_postmean)) / (iteration - params['burnin']) + + prev_x_postmean = x_postmean + + # stop of the algorithm + if (iteration > params['min_iter']) and (delta < params['threshold']): + break + + # Empirical average \approx POSTMEAN Eq. 44 + x_postmean = x_postmean / (iteration - params['burnin']) + x_postmean = uft.uirfft2(x_postmean) + + return (x_postmean, {'noise': gn_chain, 'prior': gx_chain}) + + +def richardson_lucy(data, psf, iterations=50): + """Richardson lucy deconvolution + + This a python conversion of the wikipedia page on this algorithm. See + + http://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution + + Parameters + ---------- + data : ndarray + The data + + psf : ndarray + The impulsionnal response in real space. + + iterations : int + dictionary of gibbs parameters. See below. + + Returns + ------- + im_deconv : ndarray + The deconvolued image + + References + ---------- + .. [2] Richardson, William Hadley, "Bayesian-Based Iterative + Method of Image Restoration". JOSA 62 (1): + 55–59. doi:10.1364/JOSA.62.000055, 1972 + + """ + data = data.astype(np.float) + psf = psf.astype(np.float) + im_deconv = 0.5 * np.ones(data.shape) + psf_mirror = psf[::-1, ::-1] + for iteration in range(iterations): + relative_blur = data / conv2(im_deconv, psf, 'same') + im_deconv *= conv2(relative_blur, psf_mirror, 'same') + + return im_deconv diff --git a/skimage/setup.py b/skimage/setup.py index 14f66897..75b5fa0f 100644 --- a/skimage/setup.py +++ b/skimage/setup.py @@ -12,6 +12,7 @@ def configuration(parent_package='', top_path=None): config.add_subpackage('draw') config.add_subpackage('exposure') config.add_subpackage('feature') + config.add_subpackage('deconvolution') config.add_subpackage('filter') config.add_subpackage('graph') config.add_subpackage('io')