mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-03 08:19:01 +08:00
Add deconvolution module to skimage.
This module add three function to skimage. The `wiener` function is a simple wiener deconvolution. The `unsupervised_wiener` is a more sophisticated wiener deconvolution with automatic estimation of regularisation parameters. The third function is a literal traduction in python of the rychardson lucy deconvolution of wikipedia.
This commit is contained in:
François Orieux
@@ -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()
|
||||
@@ -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"]
|
||||
Binary file not shown.
Binary file not shown.
@@ -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
|
||||
@@ -0,0 +1,473 @@
|
||||
# -*- coding: utf-8 -*-
|
||||
# uft.py --- Unitary fourier transform
|
||||
|
||||
# Copyright (c) 2011, 2012, 2013 François Orieux <orieux@iap.fr>
|
||||
|
||||
# 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 <orieux@iap.fr>"
|
||||
__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
|
||||
@@ -0,0 +1,325 @@
|
||||
# -*- coding: utf-8 -*-
|
||||
|
||||
# Copyright (c) 2013 François Orieux <orieux@iap.fr>
|
||||
|
||||
# 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 <orieux@iap.fr>"
|
||||
__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
|
||||
@@ -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')
|
||||
|
||||
Reference in New Issue
Block a user