Merge pull request #1174 from jni/deconv-shape-fix

Fix inconsistent deconvolution output shape. Closes gh-1172.
This commit is contained in:
Stefan van der Walt
2014-09-26 17:05:43 +02:00
3 changed files with 132 additions and 102 deletions
+25 -27
View File
@@ -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].
@@ -130,7 +129,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))
@@ -143,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
@@ -153,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
-------
@@ -217,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.
@@ -320,7 +319,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)
@@ -337,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
--------
@@ -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')
@@ -62,3 +86,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()
+78 -75
View File
@@ -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
--------