wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-21 12:50:27 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #800 from orieux/wiener-deconvolution

Browse Source 
Add Wiener deconvolution
This commit is contained in:
Stefan van der Walt
2013-12-13 18:42:41 -08:00
parent 0d5b822110 9eb878e440
commit 7ecd6946e8
12 changed files with 982 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CONTRIBUTORS.txt
+3
View File
@@ -173,3 +173,6 @@
- Gregor Thalhammer
Phase unwrapping integration
- François Orieux
Image deconvolution http://research.orieux.fr
doc/examples/plot_restoration.py
+60
View File
@@ -0,0 +1,60 @@
# -*- 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 based on
linear models that can't restore sharp edge as much as non-linear
methods (like TV restoration) but are much faster.
Wiener filter
-------------
The inverse filter based on the PSF (Point Spread Function),
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. The algorithm is based on a iterative Gibbs sampler that
draw alternatively samples of posterior conditionnal law of the image,
the noise power and the image frequency power.
.. [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, restoration
lena = color.rgb2gray(data.lena())
from scipy.signal import convolve2d as conv2
psf = np.ones((5, 5)) / 25
lena = conv2(lena, psf, 'same')
lena += 0.1 * lena.std() * np.random.standard_normal(lena.shape)
deconvolved, _ = restoration.unsupervised_wiener(lena, psf)
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8, 5))
plt.gray()
ax[0].imshow(lena, vmin=deconvolved.min(), vmax=deconvolved.max())
ax[0].axis('off')
ax[0].set_title('Data')
ax[1].imshow(deconvolved)
ax[1].axis('off')
ax[1].set_title('Self tuned restoration')
fig.subplots_adjust(wspace=0.02, hspace=0.2,
top=0.9, bottom=0.05, left=0, right=1)
plt.show()
skimage/__init__.py
+2
View File
@@ -29,6 +29,8 @@ measure
Measurement of image properties, e.g., similarity and contours.
morphology
Morphological operations, e.g. opening or skeletonization.
restoration
Deconvolution algorithms.
segmentation
Splitting an image into self-similar regions.
transform
skimage/restoration/__init__.py
+25
View File
@@ -0,0 +1,25 @@
# -*- coding: utf-8 -*-
"""Image restoration module.
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 .deconvolution import wiener, unsupervised_wiener, richardson_lucy
__all__ = ['wiener',
"unsupervised_wiener",
"richardson_lucy"]
skimage/restoration/deconvolution.py
+382
View File
@@ -0,0 +1,382 @@
# -*- coding: utf-8 -*-
# deconvolution.py --- Image deconvolution
"""Implementations restoration functions"""
from __future__ import division
import numpy as np
import numpy.random as npr
from scipy.signal import convolve2d
from . import uft
__keywords__ = "restoration, image, deconvolution"
def wiener(image, psf, balance, reg=None, is_real=True, clip=True):
"""Wiener-Hunt deconvolution
Return the deconvolution with a Wiener-Hunt approach (i.e. with
Fourier diagonalisation).
Parameters
----------
image : (M, N) ndarray
Input degraded image
psf : ndarray
Point Spread Function. This is assumed to be the impulse
response (input image space) if the data-type is real, or the
transfer function (Fourier space) if the data-type is
complex. There is no constraints on the shape of the impulse
response. The transfer function must be of shape `(M, N)` if
`is_real is True`, `(M, N // 2 + 1)` otherwise (see
`np.fft.rfftn`).
balance : float
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).
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.
is_real : boolean, optional
True by default. Specify if ``psf`` and ``reg`` are provided
with hermitian hypothesis, that is only half of the frequency
plane is provided (due to the redundancy of Fourier transform
of real signal). It's apply only if ``psf`` and/or ``reg`` are
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.
Returns
-------
im_deconv : (M, N) ndarray
The deconvolved image
Examples
--------
>>> from skimage import color, data, restoration
>>> lena = color.rgb2gray(data.lena())
>>> from scipy.signal import convolve2d
>>> psf = np.ones((5, 5)) / 25
>>> lena = convolve2d(lena, psf, 'same')
>>> lena += 0.1 * lena.std() * np.random.standard_normal(lena.shape)
>>> deconvolved_lena = restoration.wiener(lena, psf, 1100)
Notes
-----
This function applies the Wiener filter to a noisy and degraded
image by an impulse response (or PSF). If the data model is
.. math:: y = Hx + n
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
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 [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
those coming from noise), and the regularization.
These methods are then specific to a prior model. Consequently,
the application or the true image nature must corresponds to the
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).
Finally, the use of Fourier space implies a circulant property of
:math:`H`, see [Hunt].
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
http://research.orieux.fr/files/papers/OGR-JOSA10.pdf
.. [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, is_real=is_real)
if not np.iscomplexobj(reg):
reg = uft.ir2tf(reg, image.shape, is_real=is_real)
if psf.shape != reg.shape:
trans_func = uft.ir2tf(psf, image.shape, is_real=is_real)
else:
trans_func = psf
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))
else:
deconv = uft.uifft2(wiener_filter * uft.ufft2(image))
if clip:
deconv[deconv > 1] = 1
deconv[deconv < -1] = -1
return deconv
def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True,
clip=True):
"""Unsupervised Wiener-Hunt deconvolution
Return the deconvolution with a Wiener-Hunt approach, where the
hyperparameters are automatically estimated. The algorithm is a
stochastic iterative process (Gibbs sampler) described in the
reference below. See also ``wiener`` function.
Parameters
----------
image : (M, N) ndarray
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
(``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.
clip : boolean, optional
True by default. If true, pixel value of the result above 1 or
under -1 are thresholded for skimage pipeline
compatibility.
Returns
-------
x_postmean : (M, N) ndarray
The deconvolved image (the posterior mean).
chains : dict
The keys ``noise`` and ``prior`` contain the chain list of
noise and prior precision respectively.
Other parameters
----------------
The keys of ``user_params`` are:
threshold : float
The stopping criterion: the norm of the difference between to
successive approximated solution (empirical mean of object
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.
min_iter : int
The minimum number of iterations. 30 by default.
max_iter : int
The maximum number of iterations if ``threshold`` is not
satisfied. 150 by default.
callback : callable (None by default)
A user provided callable to which is passed, if the function
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
--------
>>> from skimage import color, data, restoration
>>> lena = color.rgb2gray(data.lena())
>>> from scipy.signal import convolve2d
>>> psf = np.ones((5, 5)) / 25
>>> lena = convolve2d(lena, psf, 'same')
>>> 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 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
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
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
http://research.orieux.fr/files/papers/OGR-JOSA10.pdf
"""
params = {'threshold': 1e-4, 'max_iter': 200,
'min_iter': 30, 'burnin': 15, 'callback': None}
params.update(user_params or {})
if reg is None:
reg, _ = uft.laplacian(image.ndim, image.shape, is_real=is_real)
if not np.iscomplexobj(reg):
reg = uft.ir2tf(reg, image.shape, is_real=is_real)
if psf.shape != reg.shape:
trans_fct = uft.ir2tf(psf, image.shape, is_real=is_real)
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
# The Fourier transfrom may change the image.size attribut, so we
# store it.
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).
# 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) +
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 * 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(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,
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'])
if is_real:
x_postmean = uft.uirfft2(x_postmean)
else:
x_postmean = uft.uifft2(x_postmean)
if clip:
x_postmean[x_postmean > 1] = 1
x_postmean[x_postmean < -1] = -1
return (x_postmean, {'noise': gn_chain, 'prior': gx_chain})
def richardson_lucy(image, psf, iterations=50, clip=True):
"""Richardson-Lucy deconvolution.
Parameters
----------
image : ndarray
Input degraded image
psf : ndarray
The point spread function
iterations : int
Number of iterations. This parameter play to 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.
Returns
-------
im_deconv : ndarray
The deconvolved image
Examples
--------
>>> from skimage import color, data, restoration
>>> camera = color.rgb2gray(data.camera())
>>> from scipy.signal import convolve2d
>>> psf = np.ones((5, 5)) / 25
>>> camera = convolve2d(camera, psf, 'same')
>>> camera += 0.1 * camera.std() * np.random.standard_normal(camera.shape)
>>> deconvolved = restoration.richardson_lucy(camera, psf, 5)
References
----------
.. [2] http://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution
"""
image = image.astype(np.float)
psf = psf.astype(np.float)
im_deconv = 0.5 * np.ones(image.shape)
psf_mirror = psf[::-1, ::-1]
for _ in range(iterations):
relative_blur = image / convolve2d(im_deconv, psf, 'same')
im_deconv *= convolve2d(relative_blur, psf_mirror, 'same')
if clip:
im_deconv[im_deconv > 1] = 1
im_deconv[im_deconv < -1] = -1
return im_deconv
skimage/restoration/tests/camera_rl.npy
BIN
View File
Binary file not shown.
skimage/restoration/tests/camera_unsup.npy
BIN
View File
Binary file not shown.
skimage/restoration/tests/camera_unsup2.npy
BIN
View File
Binary file not shown.
skimage/restoration/tests/camera_wiener.npy
BIN
View File
Binary file not shown.
skimage/restoration/tests/test_restoration.py
+64
View File
@@ -0,0 +1,64 @@
from os.path import abspath, dirname, join as pjoin
import numpy as np
from scipy.signal import convolve2d
import skimage
from skimage.data import camera
from skimage import restoration
from skimage.restoration import uft
test_img = skimage.img_as_float(camera())
def test_wiener():
psf = np.ones((5, 5)) / 25
data = convolve2d(test_img, psf, 'same')
np.random.seed(0)
data += 0.1 * data.std() * np.random.standard_normal(data.shape)
deconvolved = restoration.wiener(data, psf, 0.05)
path = pjoin(dirname(abspath(__file__)), 'camera_wiener.npy')
np.testing.assert_allclose(deconvolved, np.load(path), rtol=1e-3)
_, laplacian = uft.laplacian(2, data.shape)
otf = uft.ir2tf(psf, data.shape, is_real=False)
deconvolved = restoration.wiener(data, otf, 0.05,
reg=laplacian,
is_real=False)
np.testing.assert_allclose(np.real(deconvolved),
np.load(path),
rtol=1e-3)
def test_unsupervised_wiener():
psf = np.ones((5, 5)) / 25
data = convolve2d(test_img, psf, 'same')
np.random.seed(0)
data += 0.1 * data.std() * np.random.standard_normal(data.shape)
deconvolved, _ = restoration.unsupervised_wiener(data, psf)
path = pjoin(dirname(abspath(__file__)), 'camera_unsup.npy')
np.testing.assert_allclose(deconvolved, np.load(path), rtol=1e-3)
_, laplacian = uft.laplacian(2, data.shape)
otf = uft.ir2tf(psf, data.shape, is_real=False)
np.random.seed(0)
deconvolved = restoration.unsupervised_wiener(
data, otf, reg=laplacian, is_real=False,
user_params={"callback": lambda x: None})[0]
path = pjoin(dirname(abspath(__file__)), 'camera_unsup2.npy')
np.testing.assert_allclose(np.real(deconvolved),
np.load(path),
rtol=1e-3)
def test_richardson_lucy():
psf = np.ones((5, 5)) / 25
data = convolve2d(test_img, psf, 'same')
np.random.seed(0)
data += 0.1 * data.std() * np.random.standard_normal(data.shape)
deconvolved = restoration.richardson_lucy(data, psf, 5)
path = pjoin(dirname(abspath(__file__)), 'camera_rl.npy')
np.testing.assert_allclose(deconvolved, np.load(path), rtol=1e-3)
skimage/restoration/uft.py
+445
View File
@@ -0,0 +1,445 @@
# -*- coding: utf-8 -*-
# uft.py --- Unitary fourier transform
"""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
.. math:: \frac{1}{\sqrt{n}} \sum_i x_i
or the Fourier tranform have the same energy than 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.
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
"""
from __future__ import division, print_function
import numpy as np
__keywords__ = "fft, Fourier Transform, orthonormal, unitary"
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 : ndarray (same shape than inarray)
The unitary N-D Fourier transform of ``inarray``.
Examples
--------
>>> input = np.ones((3, 3, 3))
>>> output = ufftn(input)
>>> np.allclose(np.sum(input) / np.sqrt(input.size), output[0, 0, 0])
True
>>> output.shape
(3, 3, 3)
"""
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:]))
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 : ndarray (same shape than inarray)
The unitary inverse N-D Fourier transform of ``inarray``.
Examples
--------
>>> input = np.ones((3, 3, 3))
>>> output = uifftn(input)
>>> np.allclose(np.sum(input) / np.sqrt(input.size), output[0, 0, 0])
True
>>> output.shape
(3, 3, 3)
"""
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:]))
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 : ndarray (the last dim as N / 2 + 1 lenght)
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.
Examples
--------
>>> input = np.ones((5, 5, 5))
>>> output = urfftn(input)
>>> np.allclose(np.sum(input) / np.sqrt(input.size), output[0, 0, 0])
True
>>> output.shape
(5, 5, 3)
"""
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:]))
def uirfftn(inarray, dim=None, shape=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.
shape : tuple of int
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``.
Returns
-------
outarray : ndarray
The unitary N-D inverse real Fourier transform of ``inarray``.
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.
Examples
--------
>>> input = np.ones((5, 5, 5))
>>> output = uirfftn(urfftn(input), shape=input.shape)
>>> np.allclose(input, output)
True
>>> output.shape
(5, 5, 5)
"""
if dim is None:
dim = inarray.ndim
outarray = np.fft.irfftn(inarray, shape, axes=range(-dim, 0))
return outarray * np.sqrt(np.prod(outarray.shape[-dim:]))
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 : ndarray (same shape than inarray)
The unitary 2-D Fourier transform of ``inarray``.
See Also
--------
uifft2, ufftn, urfftn
Examples
--------
>>> input = np.ones((10, 128, 128))
>>> output = ufft2(input)
>>> np.allclose(np.sum(input[1, ...]) / np.sqrt(input[1, ...].size), output[1, 0, 0])
True
>>> output.shape
(10, 128, 128)
"""
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 : ndarray (same shape than inarray)
The unitary 2-D inverse Fourier transform of ``inarray``.
See Also
--------
uifft2, uifftn, uirfftn
Examples
--------
>>> input = np.ones((10, 128, 128))
>>> output = uifft2(input)
>>> np.allclose(np.sum(input[1, ...]) / np.sqrt(input[1, ...].size), output[0, 0, 0])
True
>>> output.shape
(10, 128, 128)
"""
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 : ndarray (the last dim as (N - 1) *2 lenght)
The unitary 2-D real Fourier transform of ``inarray``.
See Also
--------
ufft2, ufftn, urfftn
Examples
--------
>>> input = np.ones((10, 128, 128))
>>> output = urfft2(input)
>>> np.allclose(np.sum(input[1,...]) / np.sqrt(input[1,...].size), output[1, 0, 0])
True
>>> output.shape
(10, 128, 65)
"""
return urfftn(inarray, 2)
def uirfft2(inarray, shape=None):
"""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 : ndarray (the last dim as (N - 1) *2 lenght)
The unitary 2-D inverse real Fourier transform of ``inarray``.
See Also
--------
urfft2, uifftn, uirfftn
Examples
--------
>>> input = np.ones((10, 128, 128))
>>> output = uirfftn(urfftn(input), shape=input.shape)
>>> np.allclose(input, output)
True
>>> output.shape
(10, 128, 128)
"""
return uirfftn(inarray, 2, shape=shape)
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 : ndarray
The images are supposed to be in the last two axes
Returns
-------
norm : float
The quadratic norm of ``inarray``.
Examples
--------
>>> input = np.ones((5, 5))
>>> image_quad_norm(ufft2(input)) == np.sum(np.abs(input)**2)
True
>>> image_quad_norm(ufft2(input)) == image_quad_norm(urfft2(input))
True
"""
# 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 ir2tf(imp_resp, shape, dim=None, is_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.
is_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
Examples
--------
>>> 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
-----
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
# Roll for zero convention of the fft to avoid the phase
# problem. Work with odd and even size.
for axis, axis_size in enumerate(imp_resp.shape):
if axis >= imp_resp.ndim - dim:
irpadded = np.roll(irpadded,
shift=-int(np.floor(axis_size / 2)),
axis=axis)
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, is_real=True):
"""Return the transfer 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 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.
Returns
-------
tf : array_like, complex
The transfer function
impr : array_like, real
The Laplacian
Examples
--------
>>> tf, ir = laplacian(2, (32, 32))
>>> np.all(ir == np.array([[0, -1, 0], [-1, 4, -1], [0, -1, 0]]))
True
>>> np.all(tf == ir2tf(ir, (32, 32)))
True
"""
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, is_real=is_real), impr
skimage/setup.py
+1
View File
@@ -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('restoration')
config.add_subpackage('filter')
config.add_subpackage('graph')
config.add_subpackage('io')