wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-24 13:20:43 +08:00
Code Issues Packages Projects Releases Wiki Activity

PR fixs.

Browse Source
This commit is contained in:
François Orieux
2013-12-10 22:45:15 +01:00
parent 8206695afd
commit e67bdd3471
2 changed files with 65 additions and 41 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/restoration/deconvolution.py
+31 -17
View File
@@ -41,7 +41,7 @@ __keywords__ = "restoration, image, deconvolution"
def wiener(image, psf, balance, reg=None, is_real=True):
"""Wiener-Hunt deconvolution
Return the deconvolution with a wiener-Hunt approach (ie with
Return the deconvolution with a Wiener-Hunt approach (ie with
Fourier diagonalisation).
Parameters
@@ -151,8 +151,9 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None):
"""Unsupervised Wiener-Hunt deconvolution
Return the deconvolution with a Wiener-Hunt approach, where the
hyperparameters are estimated. The algorithm is a stochastic
iterative process (Gibbs sampler) described in [1].
hyperparameters are automatically estimated. The algorithm is a
stochastic iterative process (Gibbs sampler) described in [1]. See
also `wiener` function.
Parameters
----------
@@ -183,7 +184,7 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None):
threshold : float
The stopping criterion: the norm of the difference between to
successive approximated solution (empirical mean of object
samples). 1e-4 by default.
samples, see Notes section). 1e-4 by default.
burnin : int
The number of sample to ignore to start computation of the
mean. 100 by default.
@@ -192,11 +193,12 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None):
max_iter : int
The maximum number of iterations if `threshold` is not
satisfied. 150 by default.
callback : None
callback : callable (None by default)
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. It
has no influence on the algorithm execution.
exists, the current image sample for whatever purpose. The user
can store the sample, or compute other moments than the
mean. It has no influence on the algorithm execution and is
only for inspection.
Examples
--------
@@ -208,6 +210,18 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None):
>>> lena += 0.1 * lena.std() * np.random.standard_normal(lena.shape)
>>> deconvolved_lena = restoration.unsupervised_wiener(lena, psf)
Notes
-----
The estimated image is design as the posterior mean of a
probability law (from a Bayesian analysis). The mean is defined as
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 consider low probable
image useless in the sum. Finally the empirical mean of these
samples give us an estimation of the mean, and an exact
computation with an infinite sample set.
References
----------
.. [1] François Orieux, Jean-François Giovannelli, and Thomas
@@ -223,7 +237,7 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None):
'min_iter': 30, 'burnin': 15, 'callback': None}
params.update(user_params or {})
if not reg:
if reg is None:
reg, _ = uft.laplacian(image.ndim, image.shape)
if not np.iscomplex(reg):
reg = uft.ir2tf(reg, image.shape)
@@ -251,8 +265,7 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None):
# The Fourier transfrom may change the image.size attribut, so we
# store it.
image_size = image.size
image = uft.urfft2(image.astype(np.float))
data_spectrum = uft.urfft2(image.astype(np.float))
# Gibbs sampling
for iteration in range(params['max_iter']):
@@ -261,24 +274,25 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None):
# weighing (correlation in direct space)
precision = gn_chain[-1] * atf2 + gx_chain[-1] * areg2 # Eq. 29
excursion = np.sqrt(0.5) / np.sqrt(precision) * (
np.random.standard_normal(image.shape) +
1j * np.random.standard_normal(image.shape))
np.random.standard_normal(data_spectrum.shape) +
1j * np.random.standard_normal(data_spectrum.shape))
# 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 * image + excursion
x_sample = wiener_filter * data_spectrum + excursion
if params['callback']:
params['callback'](x_sample)
# sample of Eq. 31 p(gn | x^k, gx^k, y)
gn_chain.append(npr.gamma(image_size / 2,
2 / uft.image_quad_norm(image - x_sample *
gn_chain.append(npr.gamma(image.size / 2,
2 / uft.image_quad_norm(data_spectrum -
x_sample *
trans_fct)))
# sample of Eq. 31 p(gx | x^k, gn^k-1, y)
gx_chain.append(npr.gamma((image_size - 1) / 2,
gx_chain.append(npr.gamma((image.size - 1) / 2,
2 / uft.image_quad_norm(x_sample * reg)))
# current empirical average
skimage/restoration/uft.py
+34 -24
View File
@@ -24,7 +24,7 @@
"""Function of unitary fourier transform and utilities
This module implement unitary fourier transform, that is ortho-normal
transform. They are especially and usefull for convolution [1]: they
transform. They are especially and useful for convolution [1]: they
respect the Parseval equality, the value of the null frequency is
equal to
@@ -123,13 +123,12 @@ def ufftn(inarray, dim=None):
Returns
-------
outarray : array-like (same shape than inarray)
outarray : ndarray (same shape than inarray)
The unitary N-D Fourier transform of `inarray`.
"""
if dim is None:
dim = inarray.ndim
outarray = np.fft.fftn(inarray, axes=range(-dim, 0))
return outarray / np.sqrt(np.prod(inarray.shape[-dim:]))
@@ -146,13 +145,12 @@ def uifftn(inarray, dim=None):
Returns
-------
outarray : array-like (same shape than inarray)
outarray : ndarray (same shape than inarray)
The unitary inverse N-D Fourier transform of `inarray`.
"""
if not dim:
if dim is None:
dim = inarray.ndim
outarray = np.fft.ifftn(inarray, axes=range(-dim, 0))
return outarray * np.sqrt(np.prod(inarray.shape[-dim:]))
@@ -172,13 +170,18 @@ def urfftn(inarray, dim=None):
Returns
-------
outarray : array-like (the last dim as N / 2 + 1 lenght)
outarray : ndarray (the last dim as N / 2 + 1 lenght)
The unitary N-D real Fourier transform of `inarray`.
Notes
-----
The `r` function assume an input array of real
values. Consequently, the output have an Hermitian property and
redondant values are not computed and returned.
"""
if not dim:
if dim is None:
dim = inarray.ndim
outarray = np.fft.rfftn(inarray, axes=range(-dim, 0))
return outarray / np.sqrt(np.prod(inarray.shape[-dim:]))
@@ -198,13 +201,18 @@ def uirfftn(inarray, dim=None):
Returns
-------
outarray : array-like (the last dim as (N - 1) *2 lenght)
outarray : ndarray (the last dim as (N - 1) *2 lenght)
The unitary N-D inverse real Fourier transform of `inarray`.
Notes
-----
The `r` function assume an input array of real
values. Consequently, the output have an Hermitian property and
redondant values are not computed and returned.
"""
if not dim:
if dim is None:
dim = inarray.ndim
outarray = np.fft.irfftn(inarray, axes=range(-dim, 0))
return outarray * np.sqrt(np.prod(inarray.shape[-dim:-1]) *
(inarray.shape[-1] - 1) * 2)
@@ -221,7 +229,8 @@ def ufft2(inarray):
Returns
-------
outarray : array-like (same shape than inarray)
outarray : ndarray (same shape than inarray)
The unitary 2-D Fourier transform of `inarray`.
See Also
--------
@@ -242,7 +251,8 @@ def uifft2(inarray):
Returns
-------
outarray : array-like (same shape than inarray)
outarray : ndarray (same shape than inarray)
The unitary 2-D inverse Fourier transform of `inarray`.
See Also
--------
@@ -265,7 +275,8 @@ def urfft2(inarray):
Returns
-------
outarray : array-like (the last dim as (N - 1) *2 lenght)
outarray : ndarray (the last dim as (N - 1) *2 lenght)
The unitary 2-D real Fourier transform of `inarray`.
See Also
--------
@@ -288,7 +299,8 @@ def uirfft2(inarray):
Returns
-------
outarray : array-like (the last dim as (N - 1) *2 lenght)
outarray : ndarray (the last dim as (N - 1) *2 lenght)
The unitary 2-D inverse real Fourier transform of `inarray`.
See Also
--------
@@ -304,12 +316,13 @@ def image_quad_norm(inarray):
Parameters
----------
inarray : array-like
inarray : ndarray
The images are supposed to be in the last two axes
Returns
-------
norm : float
The quadratic norm of `inarray`.
"""
# If there is an hermitian symmetry
@@ -361,7 +374,6 @@ def ir2tf(imp_resp, shape, dim=None, real=True):
"""
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
@@ -372,7 +384,6 @@ def ir2tf(imp_resp, shape, dim=None, real=True):
if i >= imp_resp.ndim - dim
else 0
for i, s in enumerate(imp_resp.shape)])
if real:
return np.fft.rfftn(irpadded, axes=range(-dim, 0))
else:
@@ -409,5 +420,4 @@ def laplacian(ndim, shape):
-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