From 214613c2a8b248331db619ea75298a0e4fdc79a2 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Fri, 26 Sep 2014 18:31:56 +1000 Subject: [PATCH 1/5] Run restoration test suite when run as main --- skimage/restoration/tests/test_restoration.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/skimage/restoration/tests/test_restoration.py b/skimage/restoration/tests/test_restoration.py index 915640af..0e70af56 100644 --- a/skimage/restoration/tests/test_restoration.py +++ b/skimage/restoration/tests/test_restoration.py @@ -62,3 +62,8 @@ def test_richardson_lucy(): path = pjoin(dirname(abspath(__file__)), 'camera_rl.npy') np.testing.assert_allclose(deconvolved, np.load(path), rtol=1e-3) + + +if __name__ == '__main__': + from numpy import testing + testing.run_module_suite() From d50afed18ebdf2b312745f9dc8eb2937cf2fd084 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Fri, 26 Sep 2014 19:09:13 +1000 Subject: [PATCH 2/5] Add test for preserving image shape in deconv I can confirm that this test does not pass in the current master. --- skimage/restoration/tests/test_restoration.py | 24 +++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/skimage/restoration/tests/test_restoration.py b/skimage/restoration/tests/test_restoration.py index 0e70af56..f248f613 100644 --- a/skimage/restoration/tests/test_restoration.py +++ b/skimage/restoration/tests/test_restoration.py @@ -2,6 +2,7 @@ from os.path import abspath, dirname, join as pjoin import numpy as np from scipy.signal import convolve2d +from scipy import ndimage as nd import skimage from skimage.data import camera @@ -53,6 +54,29 @@ def test_unsupervised_wiener(): rtol=1e-3) +def test_image_shape(): + """Test that shape of output image in deconvolution is same as input. + + This addresses issue #1172. + """ + point = np.zeros((5, 5), np.float) + point[2, 2] = 1. + psf = nd.gaussian_filter(point, sigma=1.) + # image shape: (45, 45), as reported in #1172 + image = skimage.img_as_float(camera()[110:155, 225:270]) # just the face + image_conv = nd.convolve(image, psf) + deconv_sup = restoration.wiener(image_conv, psf, 1) + deconv_un = restoration.unsupervised_wiener(image_conv, psf)[0] + # test the shape + np.testing.assert_equal(image.shape, deconv_sup.shape) + np.testing.assert_equal(image.shape, deconv_un.shape) + # test the reconstruction error + sup_relative_error = np.abs(deconv_sup - image) / image + un_relative_error = np.abs(deconv_un - image) / image + np.testing.assert_array_less(np.median(sup_relative_error), 0.1) + np.testing.assert_array_less(np.median(un_relative_error), 0.1) + + def test_richardson_lucy(): psf = np.ones((5, 5)) / 25 data = convolve2d(test_img, psf, 'same') From 2b7bf24c1fb8abfdd32505cf1045d2038a89a1d4 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Fri, 26 Sep 2014 19:11:32 +1000 Subject: [PATCH 3/5] Add shape fix for deconvolution functions --- skimage/restoration/deconvolution.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/skimage/restoration/deconvolution.py b/skimage/restoration/deconvolution.py index 51c2c705..ade2516d 100644 --- a/skimage/restoration/deconvolution.py +++ b/skimage/restoration/deconvolution.py @@ -130,7 +130,8 @@ def wiener(image, psf, balance, reg=None, is_real=True, clip=True): wiener_filter = np.conj(trans_func) / (np.abs(trans_func)**2 + balance * np.abs(reg)**2) if is_real: - deconv = uft.uirfft2(wiener_filter * uft.urfft2(image)) + deconv = uft.uirfft2(wiener_filter * uft.urfft2(image), + shape=image.shape) else: deconv = uft.uifft2(wiener_filter * uft.ufft2(image)) @@ -320,7 +321,7 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True, # Empirical average \approx POSTMEAN Eq. 44 x_postmean = x_postmean / (iteration - params['burnin']) if is_real: - x_postmean = uft.uirfft2(x_postmean) + x_postmean = uft.uirfft2(x_postmean, shape=image.shape) else: x_postmean = uft.uifft2(x_postmean) From 141d459b71bbbe3f0a33024dd35b3e51c4f73d04 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Fri, 26 Sep 2014 19:21:59 +1000 Subject: [PATCH 4/5] Fix small errors in the documentation --- skimage/restoration/deconvolution.py | 47 +++++++++++++--------------- 1 file changed, 22 insertions(+), 25 deletions(-) diff --git a/skimage/restoration/deconvolution.py b/skimage/restoration/deconvolution.py index ade2516d..60be08ca 100644 --- a/skimage/restoration/deconvolution.py +++ b/skimage/restoration/deconvolution.py @@ -36,11 +36,11 @@ def wiener(image, psf, balance, reg=None, is_real=True, clip=True): The regularisation parameter value that tunes the balance between the data adequacy that improve frequency restoration and the prior adequacy that reduce frequency restoration (to - avoid noise artifact). + avoid noise artifacts). reg : ndarray, optional The regularisation operator. The Laplacian by default. It can be an impulse response or a transfer function, as for the - psf. Shape constraint is the same than for the `psf` parameter. + psf. Shape constraint is the same as for the `psf` parameter. is_real : boolean, optional True by default. Specify if ``psf`` and ``reg`` are provided with hermitian hypothesis, that is only half of the frequency @@ -49,14 +49,13 @@ def wiener(image, psf, balance, reg=None, is_real=True, clip=True): provided as transfer function. For the hermitian property see ``uft`` module or ``np.fft.rfftn``. clip : boolean, optional - True by default. If true, pixel value of the result above 1 or - under -1 are thresholded for skimage pipeline - compatibility. + True by default. If True, pixel values of the result above 1 or + under -1 are thresholded for skimage pipeline compatibility. Returns ------- im_deconv : (M, N) ndarray - The deconvolved image + The deconvolved image. Examples -------- @@ -96,8 +95,8 @@ def wiener(image, psf, balance, reg=None, is_real=True, clip=True): prior model. By default, the prior model (Laplacian) introduce image smoothness or pixel correlation. It can also be interpreted as high-frequency penalization to compensate the instability of - the solution wrt. data (sometimes called noise amplification or - "explosive" solution). + the solution with respect to the data (sometimes called noise + amplification or "explosive" solution). Finally, the use of Fourier space implies a circulant property of :math:`H`, see [Hunt]. @@ -144,7 +143,7 @@ def wiener(image, psf, balance, reg=None, is_real=True, clip=True): def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True, clip=True): - """Unsupervised Wiener-Hunt deconvolution + """Unsupervised Wiener-Hunt deconvolution. Return the deconvolution with a Wiener-Hunt approach, where the hyperparameters are automatically estimated. The algorithm is a @@ -154,21 +153,20 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True, Parameters ---------- image : (M, N) ndarray - The input degraded image + The input degraded image. psf : ndarray The impulse response (input image's space) or the transfer function (Fourier space). Both are accepted. The transfer - function is recognize as being complex + function is automatically recognized as being complex (``np.iscomplexobj(psf)``). reg : ndarray, optional The regularisation operator. The Laplacian by default. It can be an impulse response or a transfer function, as for the psf. user_params : dict - dictionary of gibbs parameters. See below. + Dictionary of parameters for the Gibbs sampler. See below. clip : boolean, optional - True by default. If true, pixel value of the result above 1 or - under -1 are thresholded for skimage pipeline - compatibility. + True by default. If true, pixel values of the result above 1 or + under -1 are thresholded for skimage pipeline compatibility. Returns ------- @@ -218,10 +216,10 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True, a sum over all the possible images weighted by their respective probability. Given the size of the problem, the exact sum is not tractable. This algorithm use of MCMC to draw image under the - posterior law. The practical idea is to only draw high probable - image since they have the biggest contribution to the mean. At the - opposite, the lowest probable image are draw less often since - their contribution are low. Finally the empirical mean of these + posterior law. The practical idea is to only draw highly probable + images since they have the biggest contribution to the mean. At the + opposite, the less probable images are drawn less often since + their contribution is low. Finally the empirical mean of these samples give us an estimation of the mean, and an exact computation with an infinite sample set. @@ -338,21 +336,20 @@ def richardson_lucy(image, psf, iterations=50, clip=True): Parameters ---------- image : ndarray - Input degraded image + Input degraded image. psf : ndarray - The point spread function + The point spread function. iterations : int - Number of iterations. This parameter play to role of + Number of iterations. This parameter plays the role of regularisation. clip : boolean, optional True by default. If true, pixel value of the result above 1 or - under -1 are thresholded for skimage pipeline - compatibility. + under -1 are thresholded for skimage pipeline compatibility. Returns ------- im_deconv : ndarray - The deconvolved image + The deconvolved image. Examples -------- From 91f1d807458543d119c2380fd0c36c45320bbbdf Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Fri, 26 Sep 2014 19:54:15 +1000 Subject: [PATCH 5/5] Fix typos and grammatical errors in uft.py docs --- skimage/restoration/uft.py | 153 +++++++++++++++++++------------------ 1 file changed, 78 insertions(+), 75 deletions(-) diff --git a/skimage/restoration/uft.py b/skimage/restoration/uft.py index fedcc6f6..9e452bef 100644 --- a/skimage/restoration/uft.py +++ b/skimage/restoration/uft.py @@ -3,17 +3,16 @@ """Function of unitary fourier transform and utilities -This module implement unitary fourier transform, that is ortho-normal -transform. They are especially and useful for convolution [1]: they -respect the Parseval equality, the value of the null frequency is -equal to +This module implements the unitary fourier transform, also known as +the ortho-normal transform. It is especially useful for convolution +[1], as it respects the Parseval equality. The value of the null +frequency is equal to .. math:: \frac{1}{\sqrt{n}} \sum_i x_i -or the Fourier tranform have the same energy than the original image +so the Fourier tranform has the same energy as the original image (see ``image_quad_norm`` function). The transform is applied from the -last axes for performance reason (c order array). You may use directly -the numpy.fft module for more sophisticated purpose. +last axis for performance (assuming a C-order array input). References ---------- @@ -31,14 +30,14 @@ __keywords__ = "fft, Fourier Transform, orthonormal, unitary" def ufftn(inarray, dim=None): - """N-dim unitary Fourier transform + """N-dimensional unitary Fourier transform. Parameters ---------- inarray : ndarray The array to transform. dim : int, optional - The ``dim`` last axis along wich to compute the transform. All + The last axis along which to compute the transform. All axes by default. Returns @@ -62,14 +61,14 @@ def ufftn(inarray, dim=None): def uifftn(inarray, dim=None): - """N-dim unitary inverse Fourier transform + """N-dimensional 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 + The last axis along which to compute the transform. All axes by default. Returns @@ -93,29 +92,29 @@ def uifftn(inarray, dim=None): def urfftn(inarray, dim=None): - """N-dim real unitary Fourier transform + """N-dimensional real unitary Fourier transform. - This transform consider the Hermitian property of the transform on - real input + This transform considers the Hermitian property of the transform on + real-valued input. Parameters ---------- - inarray : ndarray + inarray : ndarray, shape (M, N, ..., P) The array to transform. dim : int, optional - The ``dim`` last axis along wich to compute the transform. All + The last axis along which to compute the transform. All axes by default. Returns ------- - outarray : ndarray (the last dim as N / 2 + 1 lenght) + outarray : ndarray, shape (M, N, ..., P / 2 + 1) The unitary N-D real Fourier transform of ``inarray``. Notes ----- The ``urfft`` functions assume an input array of real - values. Consequently, the output have an Hermitian property and - redondant values are not computed and returned. + values. Consequently, the output has a Hermitian property and + redundant values are not computed or returned. Examples -------- @@ -133,22 +132,22 @@ def urfftn(inarray, dim=None): def uirfftn(inarray, dim=None, shape=None): - """N-dim real unitary Fourier transform + """N-dimensional inverse real unitary Fourier transform. - This transform consider the Hermitian property of the transform - from complex to real real input. + This transform considers the Hermitian property of the transform + from complex to real input. Parameters ---------- inarray : ndarray The array to transform. dim : int, optional - The ``dim`` last axis along wich to compute the transform. All + The last axis along which to compute the transform. All axes by default. - shape : tuple of int + shape : tuple of int, optional The shape of the output. The shape of ``rfft`` is ambiguous in - case of odd shape. In this case, the parameter must be - used. see ``np.fft.irfftn``. + case of odd-valued input shape. In this case, this parameter + should be provided. See ``np.fft.irfftn``. Returns ------- @@ -157,9 +156,9 @@ def uirfftn(inarray, dim=None, shape=None): Notes ----- - The ``uirfft`` function assume that output array is of real - values. Consequently, the input is assumed of having an Hermitian - property and redondant values are implicit. + The ``uirfft`` function assumes that the output array is + real-valued. Consequently, the input is assumed to have a Hermitian + property and redundant values are implicit. Examples -------- @@ -177,7 +176,7 @@ def uirfftn(inarray, dim=None, shape=None): def ufft2(inarray): - """2-dim unitary Fourier transform + """2-dimensional unitary Fourier transform. Compute the Fourier transform on the last 2 axes. @@ -188,7 +187,7 @@ def ufft2(inarray): Returns ------- - outarray : ndarray (same shape than inarray) + outarray : ndarray (same shape as inarray) The unitary 2-D Fourier transform of ``inarray``. See Also @@ -199,7 +198,8 @@ def ufft2(inarray): -------- >>> input = np.ones((10, 128, 128)) >>> output = ufft2(input) - >>> np.allclose(np.sum(input[1, ...]) / np.sqrt(input[1, ...].size), output[1, 0, 0]) + >>> np.allclose(np.sum(input[1, ...]) / np.sqrt(input[1, ...].size), + ... output[1, 0, 0]) True >>> output.shape (10, 128, 128) @@ -208,7 +208,7 @@ def ufft2(inarray): def uifft2(inarray): - """2-dim inverse unitary Fourier transform + """2-dimensional inverse unitary Fourier transform. Compute the inverse Fourier transform on the last 2 axes. @@ -219,7 +219,7 @@ def uifft2(inarray): Returns ------- - outarray : ndarray (same shape than inarray) + outarray : ndarray (same shape as inarray) The unitary 2-D inverse Fourier transform of ``inarray``. See Also @@ -230,7 +230,8 @@ def uifft2(inarray): -------- >>> input = np.ones((10, 128, 128)) >>> output = uifft2(input) - >>> np.allclose(np.sum(input[1, ...]) / np.sqrt(input[1, ...].size), output[0, 0, 0]) + >>> np.allclose(np.sum(input[1, ...]) / np.sqrt(input[1, ...].size), + ... output[0, 0, 0]) True >>> output.shape (10, 128, 128) @@ -239,20 +240,20 @@ def uifft2(inarray): def urfft2(inarray): - """2-dim real unitary Fourier transform + """2-dimensional 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. + transform considers the Hermitian property of the transform from + complex to real-valued input. Parameters ---------- - inarray : ndarray + inarray : ndarray, shape (M, N, ..., P) The array to transform. Returns ------- - outarray : ndarray (the last dim as (N - 1) *2 lenght) + outarray : ndarray, shape (M, N, ..., 2 * (P - 1)) The unitary 2-D real Fourier transform of ``inarray``. See Also @@ -263,7 +264,8 @@ def urfft2(inarray): -------- >>> input = np.ones((10, 128, 128)) >>> output = urfft2(input) - >>> np.allclose(np.sum(input[1,...]) / np.sqrt(input[1,...].size), output[1, 0, 0]) + >>> np.allclose(np.sum(input[1,...]) / np.sqrt(input[1,...].size), + ... output[1, 0, 0]) True >>> output.shape (10, 128, 65) @@ -272,20 +274,20 @@ def urfft2(inarray): def uirfft2(inarray, shape=None): - """2-dim real unitary Fourier transform + """2-dimensional inverse 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. + This transform considers the Hermitian property of the transform + from complex to real-valued input. Parameters ---------- - inarray : ndarray + inarray : ndarray, shape (M, N, ..., P) The array to transform. Returns ------- - outarray : ndarray (the last dim as (N - 1) *2 lenght) + outarray : ndarray, shape (M, N, ..., 2 * (P - 1)) The unitary 2-D inverse real Fourier transform of ``inarray``. See Also @@ -305,14 +307,16 @@ def uirfft2(inarray, shape=None): def image_quad_norm(inarray): - """Return quadratic norm of images in Fourier space + """Return the quadratic norm of images in Fourier space. - This function detect if the image suppose the hermitian property. + This function detects whether the input image satisfies the + Hermitian property. Parameters ---------- inarray : ndarray - The images are supposed to be in the last two axes + Input image. The image data should reside in the final two + axes. Returns ------- @@ -327,40 +331,40 @@ def image_quad_norm(inarray): >>> image_quad_norm(ufft2(input)) == image_quad_norm(urfft2(input)) True """ - # If there is an hermitian symmetry + # If there is a 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) + 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 ir2tf(imp_resp, shape, dim=None, is_real=True): - """Compute the transfer function of IR + """Compute the transfer function of an impulse response (IR). - This function make the necessary correct zero-padding, zero - convention, correct fft2 etc... to compute the transfer function + This function makes 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. + The impulse responses. shape : tuple of int - A tuple of integer corresponding to the target shape of the - tranfert function. + A tuple of integer corresponding to the target shape of the + tranfer function. dim : int, optional - The ``dim`` last axis along wich to compute the transform. All + The last axis along which to compute the transform. All axes by default. is_real : boolean (optionnal, default True) - If True, imp_resp is supposed real and the hermissian property + If True, imp_resp is supposed real and the Hermitian property is used with rfftn Fourier transform. Returns ------- y : complex ndarray - The tranfert function of shape ``shape``. + The tranfer function of shape ``shape``. See Also -------- @@ -377,11 +381,10 @@ def ir2tf(imp_resp, shape, dim=None, is_real=True): 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. + The input array can be composed of multiple-dimensional IR with + an arbitrary number of IR. The individual IR must be accesed + through the first axes. The last ``dim`` axes contain the space + definition. """ if not dim: dim = imp_resp.ndim @@ -402,27 +405,27 @@ def ir2tf(imp_resp, shape, dim=None, is_real=True): def laplacian(ndim, shape, is_real=True): - """Return the transfer function of the Laplacian + """Return the transfer function of the Laplacian. - Laplacian is the second order difference, on line and column. + Laplacian is the second order difference, on row and column. Parameters ---------- ndim : int - The dimension of the Laplacian + The dimension of the Laplacian. shape : tuple, shape The support on which to compute the transfer function - is_real : boolean (optionnal, default True) - If True, imp_resp is supposed real and the hermissian property - is used with rfftn Fourier transform to return the transfer - function. + is_real : boolean (optional, default True) + If True, imp_resp is assumed to be real-valued and + the Hermitian property is used with rfftn Fourier transform + to return the transfer function. Returns ------- tf : array_like, complex - The transfer function + The transfer function. impr : array_like, real - The Laplacian + The Laplacian. Examples --------