diff --git a/skimage/restoration/deconvolution.py b/skimage/restoration/deconvolution.py index 367af563..f3e37c01 100644 --- a/skimage/restoration/deconvolution.py +++ b/skimage/restoration/deconvolution.py @@ -92,12 +92,14 @@ def wiener(image, psf, balance, reg=None, is_real=True): where :math:`n` is noise, :math:`H` the PSF and :math:`x` the unknown original image, the Wiener filter is - .. math:: \hat x = F^\dag (|\Lambda_H|^2 + \lambda |\Lambda_D|^2) \Lambda_H^\dag F y + .. math:: + \hat x = F^\dag (|\Lambda_H|^2 + \lambda |\Lambda_D|^2) + \Lambda_H^\dag F y where :math:`F` and :math:`F^\dag` are the Fourier and inverse Fourier transfroms respectively, :math:`\Lambda_H` the transfer - function (or the Fourier transfrom of the PSF, see [2]) and - :math:`\Lambda_D` the filter to penalize the restored image + function (or the Fourier transfrom of the PSF, see [Hunt] below) + and :math:`\Lambda_D` the filter to penalize the restored image frequencies (Laplacian by default, that is penalization of high frequency). The parameter :math:`\lambda` tunes the balance between the data (that tends to increase high frequency, even @@ -110,9 +112,9 @@ def wiener(image, psf, balance, reg=None, is_real=True): as high-frequency penalization to compensate noise amplification or so called "explosive" solution. These methods are well interpreted by Bayesian analysis. - + Finally, the use of Fourier space implies a circulant property of - :math:`H`, see [2]. + :math:`H`, see [Hunt]. References ---------- @@ -128,15 +130,14 @@ def wiener(image, psf, balance, reg=None, is_real=True): .. [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 reg is None: - reg, _ = uft.laplacian(image.ndim, image.shape) + reg, _ = uft.laplacian(image.ndim, image.shape, is_real=is_real) if not np.iscomplexobj(reg): - reg = uft.ir2tf(reg, image.shape) + reg = uft.ir2tf(reg, image.shape, is_real=is_real) if psf.shape != reg.shape: - trans_func = uft.ir2tf(psf, image.shape) + trans_func = uft.ir2tf(psf, image.shape, is_real=is_real) else: trans_func = psf @@ -148,7 +149,7 @@ def wiener(image, psf, balance, reg=None, is_real=True): return uft.uifft2(wiener_filter * uft.ufft2(image)) -def unsupervised_wiener(image, psf, reg=None, user_params=None): +def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True): """Unsupervised Wiener-Hunt deconvolution Return the deconvolution with a Wiener-Hunt approach, where the @@ -240,12 +241,12 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None): params.update(user_params or {}) if reg is None: - reg, _ = uft.laplacian(image.ndim, image.shape) + reg, _ = uft.laplacian(image.ndim, image.shape, is_real=is_real) if not np.iscomplexobj(reg): - reg = uft.ir2tf(reg, image.shape) + reg = uft.ir2tf(reg, image.shape, is_real=is_real) if psf.shape != reg.shape: - trans_fct = uft.ir2tf(psf, image.shape) + trans_fct = uft.ir2tf(psf, image.shape, is_real=is_real) else: trans_fct = psf @@ -267,13 +268,16 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None): # The Fourier transfrom may change the image.size attribut, so we # store it. - data_spectrum = uft.urfft2(image.astype(np.float)) + if is_real: + data_spectrum = uft.urfft2(image.astype(np.float)) + else: + data_spectrum = uft.ufft2(image.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) + # weighting (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(data_spectrum.shape) + @@ -316,7 +320,10 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None): # Empirical average \approx POSTMEAN Eq. 44 x_postmean = x_postmean / (iteration - params['burnin']) - x_postmean = uft.uirfft2(x_postmean) + if is_real: + x_postmean = uft.uirfft2(x_postmean) + else: + x_postmean = uft.uifft2(x_postmean) return (x_postmean, {'noise': gn_chain, 'prior': gx_chain}) diff --git a/skimage/restoration/uft.py b/skimage/restoration/uft.py index 39ff2390..8dbd619b 100644 --- a/skimage/restoration/uft.py +++ b/skimage/restoration/uft.py @@ -343,7 +343,7 @@ def image_quad_norm(inarray): ------- norm : float The quadratic norm of `inarray`. - + Examples -------- >>> input = np.ones((5, 5)) @@ -360,7 +360,7 @@ def image_quad_norm(inarray): return np.sum(np.sum(np.abs(inarray)**2, axis=-1), axis=-1) -def ir2tf(imp_resp, shape, dim=None, real=True): +def ir2tf(imp_resp, shape, dim=None, is_real=True): """Compute the transfer function of IR This function make the necessary correct zero-padding, zero @@ -378,7 +378,7 @@ def ir2tf(imp_resp, shape, dim=None, real=True): dim : int, optional The `dim` last axis along wich to compute the transform. All axes by default. - real : boolean (optionnal, default True) + is_real : boolean (optionnal, default True) If True, imp_resp is supposed real and the hermissian property is used with rfftn Fourier transform. @@ -395,6 +395,10 @@ def ir2tf(imp_resp, shape, dim=None, real=True): -------- >>> np.all(np.array([[4, 0], [0, 0]]) == ir2tf(np.ones((2, 2)), (2, 2))) True + >>> ir2tf(np.ones((2, 2)), (512, 512)).shape == (512, 257) + True + >>> ir2tf(np.ones((2, 2)), (512, 512), is_real=False).shape == (512, 512) + True Notes ----- @@ -416,13 +420,13 @@ def ir2tf(imp_resp, shape, dim=None, real=True): irpadded = np.roll(irpadded, shift=-int(np.floor(axis_size / 2)), axis=axis) - if real: + if is_real: return np.fft.rfftn(irpadded, axes=range(-dim, 0)) else: return np.fft.fftn(irpadded, axes=range(-dim, 0)) -def laplacian(ndim, shape): +def laplacian(ndim, shape, is_real=True): """Return the transfert function of the laplacian Laplacian is the second order difference, on line and column. @@ -433,6 +437,10 @@ def laplacian(ndim, shape): The dimension of the laplacian shape : tuple, shape The support on which to compute the transfert 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 transfert + function. Returns ------- @@ -440,7 +448,7 @@ def laplacian(ndim, shape): The transfert function impr : array_like, real The laplacian - + Examples -------- >>> tf, ir = laplacian(2, (32, 32)) @@ -459,4 +467,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 + return ir2tf(impr, shape, is_real=is_real), impr