rm CondaPkg environment
This commit is contained in:
@@ -1,32 +0,0 @@
|
||||
"""Image restoration module.
|
||||
|
||||
"""
|
||||
|
||||
from .deconvolution import wiener, unsupervised_wiener, richardson_lucy
|
||||
from .unwrap import unwrap_phase
|
||||
from ._denoise import (denoise_tv_chambolle, denoise_tv_bregman,
|
||||
denoise_bilateral, denoise_wavelet, estimate_sigma)
|
||||
from ._cycle_spin import cycle_spin
|
||||
from .non_local_means import denoise_nl_means
|
||||
from .inpaint import inpaint_biharmonic
|
||||
from .j_invariant import calibrate_denoiser
|
||||
from .rolling_ball import rolling_ball, ball_kernel, ellipsoid_kernel
|
||||
|
||||
|
||||
__all__ = ['wiener',
|
||||
'unsupervised_wiener',
|
||||
'richardson_lucy',
|
||||
'unwrap_phase',
|
||||
'denoise_tv_bregman',
|
||||
'denoise_tv_chambolle',
|
||||
'denoise_bilateral',
|
||||
'denoise_wavelet',
|
||||
'denoise_nl_means',
|
||||
'estimate_sigma',
|
||||
'inpaint_biharmonic',
|
||||
'cycle_spin',
|
||||
'calibrate_denoiser',
|
||||
'rolling_ball',
|
||||
'ellipsoid_kernel',
|
||||
'ball_kernel',
|
||||
]
|
||||
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
@@ -1,152 +0,0 @@
|
||||
from itertools import product
|
||||
import numpy as np
|
||||
from .._shared import utils
|
||||
from .._shared.utils import warn
|
||||
|
||||
try:
|
||||
import dask
|
||||
dask_available = True
|
||||
except ImportError:
|
||||
dask_available = False
|
||||
|
||||
|
||||
def _generate_shifts(ndim, multichannel, max_shifts, shift_steps=1):
|
||||
"""Returns all combinations of shifts in n dimensions over the specified
|
||||
max_shifts and step sizes.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> s = list(_generate_shifts(2, False, max_shifts=(1, 2), shift_steps=1))
|
||||
>>> print(s)
|
||||
[(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2)]
|
||||
"""
|
||||
mc = int(multichannel)
|
||||
if np.isscalar(max_shifts):
|
||||
max_shifts = (max_shifts, ) * (ndim - mc) + (0, ) * mc
|
||||
elif multichannel and len(max_shifts) == ndim - 1:
|
||||
max_shifts = tuple(max_shifts) + (0, )
|
||||
elif len(max_shifts) != ndim:
|
||||
raise ValueError("max_shifts should have length ndim")
|
||||
|
||||
if np.isscalar(shift_steps):
|
||||
shift_steps = (shift_steps, ) * (ndim - mc) + (1, ) * mc
|
||||
elif multichannel and len(shift_steps) == ndim - 1:
|
||||
shift_steps = tuple(shift_steps) + (1, )
|
||||
elif len(shift_steps) != ndim:
|
||||
raise ValueError("max_shifts should have length ndim")
|
||||
|
||||
if any(s < 1 for s in shift_steps):
|
||||
raise ValueError("shift_steps must all be >= 1")
|
||||
|
||||
if multichannel and max_shifts[-1] != 0:
|
||||
raise ValueError(
|
||||
"Multichannel cycle spinning should not have shifts along the "
|
||||
"last axis.")
|
||||
|
||||
return product(*[range(0, s + 1, t) for
|
||||
s, t in zip(max_shifts, shift_steps)])
|
||||
|
||||
|
||||
@utils.channel_as_last_axis()
|
||||
def cycle_spin(x, func, max_shifts, shift_steps=1, num_workers=None,
|
||||
func_kw={}, *, channel_axis=None):
|
||||
"""Cycle spinning (repeatedly apply func to shifted versions of x).
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array-like
|
||||
Data for input to ``func``.
|
||||
func : function
|
||||
A function to apply to circularly shifted versions of ``x``. Should
|
||||
take ``x`` as its first argument. Any additional arguments can be
|
||||
supplied via ``func_kw``.
|
||||
max_shifts : int or tuple
|
||||
If an integer, shifts in ``range(0, max_shifts+1)`` will be used along
|
||||
each axis of ``x``. If a tuple, ``range(0, max_shifts[i]+1)`` will be
|
||||
along axis i.
|
||||
shift_steps : int or tuple, optional
|
||||
The step size for the shifts applied along axis, i, are::
|
||||
``range((0, max_shifts[i]+1, shift_steps[i]))``. If an integer is
|
||||
provided, the same step size is used for all axes.
|
||||
num_workers : int or None, optional
|
||||
The number of parallel threads to use during cycle spinning. If set to
|
||||
``None``, the full set of available cores are used.
|
||||
func_kw : dict, optional
|
||||
Additional keyword arguments to supply to ``func``.
|
||||
channel_axis : int or None, optional
|
||||
If None, the image is assumed to be a grayscale (single channel) image.
|
||||
Otherwise, this parameter indicates which axis of the array corresponds
|
||||
to channels.
|
||||
|
||||
.. versionadded:: 0.19
|
||||
``channel_axis`` was added in 0.19.
|
||||
|
||||
Returns
|
||||
-------
|
||||
avg_y : np.ndarray
|
||||
The output of ``func(x, **func_kw)`` averaged over all combinations of
|
||||
the specified axis shifts.
|
||||
|
||||
Notes
|
||||
-----
|
||||
Cycle spinning was proposed as a way to approach shift-invariance via
|
||||
performing several circular shifts of a shift-variant transform [1]_.
|
||||
|
||||
For a n-level discrete wavelet transforms, one may wish to perform all
|
||||
shifts up to ``max_shifts = 2**n - 1``. In practice, much of the benefit
|
||||
can often be realized with only a small number of shifts per axis.
|
||||
|
||||
For transforms such as the blockwise discrete cosine transform, one may
|
||||
wish to evaluate shifts up to the block size used by the transform.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] R.R. Coifman and D.L. Donoho. "Translation-Invariant De-Noising".
|
||||
Wavelets and Statistics, Lecture Notes in Statistics, vol.103.
|
||||
Springer, New York, 1995, pp.125-150.
|
||||
:DOI:`10.1007/978-1-4612-2544-7_9`
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> import skimage.data
|
||||
>>> from skimage import img_as_float
|
||||
>>> from skimage.restoration import denoise_wavelet, cycle_spin
|
||||
>>> img = img_as_float(skimage.data.camera())
|
||||
>>> sigma = 0.1
|
||||
>>> img = img + sigma * np.random.standard_normal(img.shape)
|
||||
>>> denoised = cycle_spin(img, func=denoise_wavelet,
|
||||
... max_shifts=3)
|
||||
|
||||
"""
|
||||
x = np.asanyarray(x)
|
||||
multichannel = channel_axis is not None
|
||||
all_shifts = _generate_shifts(x.ndim, multichannel, max_shifts,
|
||||
shift_steps)
|
||||
all_shifts = list(all_shifts)
|
||||
roll_axes = tuple(range(x.ndim))
|
||||
|
||||
def _run_one_shift(shift):
|
||||
# shift, apply function, inverse shift
|
||||
xs = np.roll(x, shift, axis=roll_axes)
|
||||
tmp = func(xs, **func_kw)
|
||||
return np.roll(tmp, tuple(-s for s in shift), axis=roll_axes)
|
||||
|
||||
if not dask_available and (num_workers is None or num_workers > 1):
|
||||
num_workers = 1
|
||||
warn('The optional dask dependency is not installed. '
|
||||
'The number of workers is set to 1. To silence '
|
||||
'this warning, install dask or explicitly set `num_workers=1` '
|
||||
'when calling the `cycle_spin` function')
|
||||
# compute a running average across the cycle shifts
|
||||
if num_workers == 1:
|
||||
# serial processing
|
||||
mean = _run_one_shift(all_shifts[0])
|
||||
for shift in all_shifts[1:]:
|
||||
mean += _run_one_shift(shift)
|
||||
mean /= len(all_shifts)
|
||||
else:
|
||||
# multithreaded via dask
|
||||
futures = [dask.delayed(_run_one_shift)(s) for s in all_shifts]
|
||||
mean = sum(futures) / len(futures)
|
||||
mean = mean.compute(num_workers=num_workers)
|
||||
return mean
|
||||
File diff suppressed because it is too large
Load Diff
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
@@ -1,416 +0,0 @@
|
||||
"""Implementations restoration functions"""
|
||||
|
||||
import numpy as np
|
||||
from scipy.signal import convolve
|
||||
|
||||
from .._shared.utils import _supported_float_type
|
||||
from . import uft
|
||||
|
||||
|
||||
def wiener(image, psf, balance, reg=None, is_real=True, clip=True):
|
||||
r"""Wiener-Hunt deconvolution
|
||||
|
||||
Return the deconvolution with a Wiener-Hunt approach (i.e. with
|
||||
Fourier diagonalisation).
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Input degraded image (can be n-dimensional).
|
||||
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
|
||||
`(N1, N2, ..., ND)` if `is_real is True`,
|
||||
`(N1, N2, ..., ND // 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 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 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
|
||||
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 values 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
|
||||
>>> img = color.rgb2gray(data.astronaut())
|
||||
>>> from scipy.signal import convolve2d
|
||||
>>> psf = np.ones((5, 5)) / 25
|
||||
>>> img = convolve2d(img, psf, 'same')
|
||||
>>> rng = np.random.default_rng()
|
||||
>>> img += 0.1 * img.std() * rng.standard_normal(img.shape)
|
||||
>>> deconvolved_img = restoration.wiener(img, psf, 0.1)
|
||||
|
||||
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^\dagger (|\Lambda_H|^2 + \lambda |\Lambda_D|^2)
|
||||
\Lambda_H^\dagger F y
|
||||
|
||||
where :math:`F` and :math:`F^\dagger` are the Fourier and inverse
|
||||
Fourier transforms respectively, :math:`\Lambda_H` the transfer
|
||||
function (or the Fourier transform 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 correspond 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 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 [2]_.
|
||||
|
||||
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)
|
||||
|
||||
https://www.osapublishing.org/josaa/abstract.cfm?URI=josaa-27-7-1593
|
||||
|
||||
https://hal.archives-ouvertes.fr/hal-00674508
|
||||
|
||||
.. [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)
|
||||
float_type = _supported_float_type(image.dtype)
|
||||
image = image.astype(float_type, copy=False)
|
||||
psf = psf.real.astype(float_type, copy=False)
|
||||
reg = reg.real.astype(float_type, copy=False)
|
||||
|
||||
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.uirfftn(wiener_filter * uft.urfftn(image),
|
||||
shape=image.shape)
|
||||
else:
|
||||
deconv = uft.uifftn(wiener_filter * uft.ufftn(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, *, random_state=None):
|
||||
"""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 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, optional
|
||||
Dictionary of parameters for the Gibbs sampler. See below.
|
||||
clip : boolean, optional
|
||||
True by default. If true, pixel values of the result above 1 or
|
||||
under -1 are thresholded for skimage pipeline compatibility.
|
||||
random_state : {None, int, `numpy.random.Generator`}, optional
|
||||
If `random_state` is None the `numpy.random.Generator` singleton is
|
||||
used.
|
||||
If `random_state` is an int, a new ``Generator`` instance is used,
|
||||
seeded with `random_state`.
|
||||
If `random_state` is already a ``Generator`` instance then that
|
||||
instance is used.
|
||||
|
||||
.. versionadded:: 0.19
|
||||
|
||||
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. 15 by default.
|
||||
min_num_iter : int
|
||||
The minimum number of iterations. 30 by default.
|
||||
max_num_iter : int
|
||||
The maximum number of iterations if ``threshold`` is not
|
||||
satisfied. 200 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
|
||||
>>> img = color.rgb2gray(data.astronaut())
|
||||
>>> from scipy.signal import convolve2d
|
||||
>>> psf = np.ones((5, 5)) / 25
|
||||
>>> img = convolve2d(img, psf, 'same')
|
||||
>>> rng = np.random.default_rng()
|
||||
>>> img += 0.1 * img.std() * rng.standard_normal(img.shape)
|
||||
>>> deconvolved_img = restoration.unsupervised_wiener(img, 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 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.
|
||||
|
||||
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)
|
||||
|
||||
https://www.osapublishing.org/josaa/abstract.cfm?URI=josaa-27-7-1593
|
||||
|
||||
https://hal.archives-ouvertes.fr/hal-00674508
|
||||
"""
|
||||
params = {'threshold': 1e-4, 'max_num_iter': 200,
|
||||
'min_num_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)
|
||||
float_type = _supported_float_type(image.dtype)
|
||||
image = image.astype(float_type, copy=False)
|
||||
psf = psf.real.astype(float_type, copy=False)
|
||||
reg = reg.real.astype(float_type, copy=False)
|
||||
|
||||
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, dtype=float_type)
|
||||
# The previous computed mean in the iterative loop
|
||||
prev_x_postmean = np.zeros(trans_fct.shape, dtype=float_type)
|
||||
|
||||
# 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 transform may change the image.size attribute, so we
|
||||
# store it.
|
||||
if is_real:
|
||||
data_spectrum = uft.urfft2(image)
|
||||
else:
|
||||
data_spectrum = uft.ufft2(image)
|
||||
|
||||
rng = np.random.default_rng(random_state)
|
||||
|
||||
# Gibbs sampling
|
||||
for iteration in range(params['max_num_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
|
||||
# Note: Use astype instead of dtype argument to standard_normal to get
|
||||
# similar random values across precisions, as needed for
|
||||
# reference data used by test_unsupervised_wiener.
|
||||
_rand1 = rng.standard_normal(data_spectrum.shape)
|
||||
_rand1 = _rand1.astype(float_type, copy=False)
|
||||
_rand2 = rng.standard_normal(data_spectrum.shape)
|
||||
_rand2 = _rand2.astype(float_type, copy=False)
|
||||
excursion = np.sqrt(0.5 / precision) * (_rand1 + 1j * _rand2)
|
||||
|
||||
# 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(rng.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(rng.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_num_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, shape=image.shape)
|
||||
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, num_iter=50, clip=True, filter_epsilon=None):
|
||||
"""Richardson-Lucy deconvolution.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Input degraded image (can be n-dimensional).
|
||||
psf : ndarray
|
||||
The point spread function.
|
||||
num_iter : int, optional
|
||||
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.
|
||||
filter_epsilon: float, optional
|
||||
Value below which intermediate results become 0 to avoid division
|
||||
by small numbers.
|
||||
|
||||
Returns
|
||||
-------
|
||||
im_deconv : ndarray
|
||||
The deconvolved image.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from skimage import img_as_float, data, restoration
|
||||
>>> camera = img_as_float(data.camera())
|
||||
>>> from scipy.signal import convolve2d
|
||||
>>> psf = np.ones((5, 5)) / 25
|
||||
>>> camera = convolve2d(camera, psf, 'same')
|
||||
>>> rng = np.random.default_rng()
|
||||
>>> camera += 0.1 * camera.std() * rng.standard_normal(camera.shape)
|
||||
>>> deconvolved = restoration.richardson_lucy(camera, psf, 5)
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] https://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution
|
||||
"""
|
||||
float_type = _supported_float_type(image.dtype)
|
||||
image = image.astype(float_type, copy=False)
|
||||
psf = psf.astype(float_type, copy=False)
|
||||
im_deconv = np.full(image.shape, 0.5, dtype=float_type)
|
||||
psf_mirror = np.flip(psf)
|
||||
|
||||
# Small regularization parameter used to avoid 0 divisions
|
||||
eps = 1e-12
|
||||
|
||||
for _ in range(num_iter):
|
||||
conv = convolve(im_deconv, psf, mode='same') + eps
|
||||
if filter_epsilon:
|
||||
relative_blur = np.where(conv < filter_epsilon, 0, image / conv)
|
||||
else:
|
||||
relative_blur = image / conv
|
||||
im_deconv *= convolve(relative_blur, psf_mirror, mode='same')
|
||||
|
||||
if clip:
|
||||
im_deconv[im_deconv > 1] = 1
|
||||
im_deconv[im_deconv < -1] = -1
|
||||
|
||||
return im_deconv
|
||||
@@ -1,332 +0,0 @@
|
||||
import numpy as np
|
||||
from scipy import sparse
|
||||
from scipy.sparse.linalg import spsolve
|
||||
import scipy.ndimage as ndi
|
||||
from scipy.ndimage import laplace
|
||||
|
||||
import skimage
|
||||
from .._shared import utils
|
||||
from ..measure import label
|
||||
from ._inpaint import _build_matrix_inner
|
||||
|
||||
|
||||
def _get_neighborhood(nd_idx, radius, nd_shape):
|
||||
bounds_lo = np.maximum(nd_idx - radius, 0)
|
||||
bounds_hi = np.minimum(nd_idx + radius + 1, nd_shape)
|
||||
return bounds_lo, bounds_hi
|
||||
|
||||
|
||||
def _get_neigh_coef(shape, center, dtype=float):
|
||||
# Create biharmonic coefficients ndarray
|
||||
neigh_coef = np.zeros(shape, dtype=dtype)
|
||||
neigh_coef[center] = 1
|
||||
neigh_coef = laplace(laplace(neigh_coef))
|
||||
|
||||
# extract non-zero locations and values
|
||||
coef_idx = np.where(neigh_coef)
|
||||
coef_vals = neigh_coef[coef_idx]
|
||||
|
||||
coef_idx = np.stack(coef_idx, axis=0)
|
||||
return neigh_coef, coef_idx, coef_vals
|
||||
|
||||
|
||||
def _inpaint_biharmonic_single_region(image, mask, out, neigh_coef_full,
|
||||
coef_vals, raveled_offsets):
|
||||
"""Solve a (sparse) linear system corresponding to biharmonic inpainting.
|
||||
|
||||
This function creates a linear system of the form:
|
||||
|
||||
``A @ u = b``
|
||||
|
||||
where ``A`` is a sparse matrix, ``b`` is a vector enforcing smoothness and
|
||||
boundary constraints and ``u`` is the vector of inpainted values to be
|
||||
(uniquely) determined by solving the linear system.
|
||||
|
||||
``A`` is a sparse matrix of shape (n_mask, n_mask) where `n_mask``
|
||||
corresponds to the number of non-zero values in ``mask`` (i.e. the number
|
||||
of pixels to be inpainted). Each row in A will have a number of non-zero
|
||||
values equal to the number of non-zero values in the biharmonic kernel,
|
||||
``neigh_coef_full``. In practice, biharmonic kernels with reduced extent
|
||||
are used at the image borders. This matrix, ``A`` is the same for all
|
||||
image channels (since the same inpainting mask is currently used for all
|
||||
channels).
|
||||
|
||||
``u`` is a dense matrix of shape ``(n_mask, n_channels)`` and represents
|
||||
the vector of unknown values for each channel.
|
||||
|
||||
``b`` is a dense matrix of shape ``(n_mask, n_channels)`` and represents
|
||||
the desired output of convolving the solution with the biharmonic kernel.
|
||||
At mask locations where there is no overlap with known values, ``b`` will
|
||||
have a value of 0. This enforces the biharmonic smoothness constraint in
|
||||
the interior of inpainting regions. For regions near the boundary that
|
||||
overlap with known values, the entries in ``b`` enforce boundary conditions
|
||||
designed to avoid discontinuity with the known values.
|
||||
"""
|
||||
|
||||
n_channels = out.shape[-1]
|
||||
radius = neigh_coef_full.shape[0] // 2
|
||||
|
||||
edge_mask = np.ones(mask.shape, dtype=bool)
|
||||
edge_mask[(slice(radius, -radius),) * mask.ndim] = 0
|
||||
boundary_mask = edge_mask * mask
|
||||
center_mask = ~edge_mask * mask
|
||||
|
||||
boundary_pts = np.where(boundary_mask)
|
||||
boundary_i = np.flatnonzero(boundary_mask)
|
||||
center_i = np.flatnonzero(center_mask)
|
||||
mask_i = np.concatenate((boundary_i, center_i))
|
||||
|
||||
center_pts = np.where(center_mask)
|
||||
mask_pts = tuple(
|
||||
[np.concatenate((b, c)) for b, c in zip(boundary_pts, center_pts)]
|
||||
)
|
||||
|
||||
# Use convolution to predetermine the number of non-zero entries in the
|
||||
# sparse system matrix.
|
||||
structure = neigh_coef_full != 0
|
||||
tmp = ndi.convolve(mask, structure, output=np.uint8, mode='constant')
|
||||
nnz_matrix = tmp[mask].sum()
|
||||
|
||||
# Need to estimate the number of zeros for the right hand side vector.
|
||||
# The computation below will slightly overestimate the true number of zeros
|
||||
# due to edge effects (the kernel itself gets shrunk in size near the
|
||||
# edges, but that isn't accounted for here). We can trim any excess entries
|
||||
# later.
|
||||
n_mask = np.count_nonzero(mask)
|
||||
n_struct = np.count_nonzero(structure)
|
||||
nnz_rhs_vector_max = n_mask - np.count_nonzero(tmp == n_struct)
|
||||
|
||||
# pre-allocate arrays storing sparse matrix indices and values
|
||||
row_idx_known = np.empty(nnz_rhs_vector_max, dtype=np.intp)
|
||||
data_known = np.zeros((nnz_rhs_vector_max, n_channels), dtype=out.dtype)
|
||||
row_idx_unknown = np.empty(nnz_matrix, dtype=np.intp)
|
||||
col_idx_unknown = np.empty(nnz_matrix, dtype=np.intp)
|
||||
data_unknown = np.empty(nnz_matrix, dtype=out.dtype)
|
||||
|
||||
# cache the various small, non-square Laplacians used near the boundary
|
||||
coef_cache = {}
|
||||
|
||||
# Iterate over masked points near the boundary
|
||||
mask_flat = mask.reshape(-1)
|
||||
out_flat = np.ascontiguousarray(out.reshape((-1, n_channels)))
|
||||
idx_known = 0
|
||||
idx_unknown = 0
|
||||
mask_pt_n = -1
|
||||
boundary_pts = np.stack(boundary_pts, axis=1)
|
||||
for mask_pt_n, nd_idx in enumerate(boundary_pts):
|
||||
# Get bounded neighborhood of selected radius
|
||||
b_lo, b_hi = _get_neighborhood(nd_idx, radius, mask.shape)
|
||||
|
||||
# Create (truncated) biharmonic coefficients ndarray
|
||||
coef_shape = tuple(b_hi - b_lo)
|
||||
coef_center = tuple(nd_idx - b_lo)
|
||||
coef_idx, coefs = coef_cache.get((coef_shape, coef_center),
|
||||
(None, None))
|
||||
if coef_idx is None:
|
||||
_ , coef_idx, coefs = _get_neigh_coef(coef_shape,
|
||||
coef_center,
|
||||
dtype=out.dtype)
|
||||
coef_cache[(coef_shape, coef_center)] = (coef_idx, coefs)
|
||||
|
||||
# compute corresponding 1d indices into the mask
|
||||
coef_idx = coef_idx + b_lo[:, np.newaxis]
|
||||
index1d = np.ravel_multi_index(coef_idx, mask.shape)
|
||||
|
||||
# Iterate over masked point's neighborhood
|
||||
nvals = 0
|
||||
for coef, i in zip(coefs, index1d):
|
||||
if mask_flat[i]:
|
||||
row_idx_unknown[idx_unknown] = mask_pt_n
|
||||
col_idx_unknown[idx_unknown] = i
|
||||
data_unknown[idx_unknown] = coef
|
||||
idx_unknown += 1
|
||||
else:
|
||||
data_known[idx_known, :] -= coef * out_flat[i, :]
|
||||
nvals += 1
|
||||
if nvals:
|
||||
row_idx_known[idx_known] = mask_pt_n
|
||||
idx_known += 1
|
||||
|
||||
# Call an efficient Cython-based implementation for all interior points
|
||||
row_start = mask_pt_n + 1
|
||||
known_start_idx = idx_known
|
||||
unknown_start_idx = idx_unknown
|
||||
nnz_rhs = _build_matrix_inner(
|
||||
# starting indices
|
||||
row_start, known_start_idx, unknown_start_idx,
|
||||
# input arrays
|
||||
center_i, raveled_offsets, coef_vals, mask_flat,
|
||||
out_flat,
|
||||
# output arrays
|
||||
row_idx_known, data_known, row_idx_unknown, col_idx_unknown,
|
||||
data_unknown
|
||||
)
|
||||
|
||||
# trim RHS vector values and indices to the exact length
|
||||
row_idx_known = row_idx_known[:nnz_rhs]
|
||||
data_known = data_known[:nnz_rhs, :]
|
||||
|
||||
# Form sparse matrix of unknown values
|
||||
sp_shape = (n_mask, out.size)
|
||||
matrix_unknown = sparse.coo_matrix(
|
||||
(data_unknown, (row_idx_unknown, col_idx_unknown)), shape=sp_shape
|
||||
).tocsr()
|
||||
|
||||
# Solve linear system for masked points
|
||||
matrix_unknown = matrix_unknown[:, mask_i]
|
||||
|
||||
# dense vectors representing the right hand side for each channel
|
||||
rhs = np.zeros((n_mask, n_channels), dtype=out.dtype)
|
||||
rhs[row_idx_known, :] = data_known
|
||||
|
||||
# set use_umfpack to False so float32 data is supported
|
||||
result = spsolve(matrix_unknown, rhs, use_umfpack=False,
|
||||
permc_spec='MMD_ATA')
|
||||
if result.ndim == 1:
|
||||
result = result[:, np.newaxis]
|
||||
|
||||
out[mask_pts] = result
|
||||
return out
|
||||
|
||||
|
||||
@utils.channel_as_last_axis()
|
||||
def inpaint_biharmonic(image, mask, *,
|
||||
split_into_regions=False, channel_axis=None):
|
||||
"""Inpaint masked points in image with biharmonic equations.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : (M[, N[, ..., P]][, C]) ndarray
|
||||
Input image.
|
||||
mask : (M[, N[, ..., P]]) ndarray
|
||||
Array of pixels to be inpainted. Have to be the same shape as one
|
||||
of the 'image' channels. Unknown pixels have to be represented with 1,
|
||||
known pixels - with 0.
|
||||
split_into_regions : boolean, optional
|
||||
If True, inpainting is performed on a region-by-region basis. This is
|
||||
likely to be slower, but will have reduced memory requirements.
|
||||
channel_axis : int or None, optional
|
||||
If None, the image is assumed to be a grayscale (single channel) image.
|
||||
Otherwise, this parameter indicates which axis of the array corresponds
|
||||
to channels.
|
||||
|
||||
.. versionadded:: 0.19
|
||||
``channel_axis`` was added in 0.19.
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : (M[, N[, ..., P]][, C]) ndarray
|
||||
Input image with masked pixels inpainted.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] S.B.Damelin and N.S.Hoang. "On Surface Completion and Image
|
||||
Inpainting by Biharmonic Functions: Numerical Aspects",
|
||||
International Journal of Mathematics and Mathematical Sciences,
|
||||
Vol. 2018, Article ID 3950312
|
||||
:DOI:`10.1155/2018/3950312`
|
||||
.. [2] C. K. Chui and H. N. Mhaskar, MRA Contextual-Recovery Extension of
|
||||
Smooth Functions on Manifolds, Appl. and Comp. Harmonic Anal.,
|
||||
28 (2010), 104-113,
|
||||
:DOI:`10.1016/j.acha.2009.04.004`
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> img = np.tile(np.square(np.linspace(0, 1, 5)), (5, 1))
|
||||
>>> mask = np.zeros_like(img)
|
||||
>>> mask[2, 2:] = 1
|
||||
>>> mask[1, 3:] = 1
|
||||
>>> mask[0, 4:] = 1
|
||||
>>> out = inpaint_biharmonic(img, mask)
|
||||
"""
|
||||
|
||||
if image.ndim < 1:
|
||||
raise ValueError('Input array has to be at least 1D')
|
||||
|
||||
multichannel = channel_axis is not None
|
||||
img_baseshape = image.shape[:-1] if multichannel else image.shape
|
||||
if img_baseshape != mask.shape:
|
||||
raise ValueError('Input arrays have to be the same shape')
|
||||
|
||||
if np.ma.isMaskedArray(image):
|
||||
raise TypeError('Masked arrays are not supported')
|
||||
|
||||
image = skimage.img_as_float(image)
|
||||
|
||||
# float16->float32 and float128->float64
|
||||
float_dtype = utils._supported_float_type(image.dtype)
|
||||
image = image.astype(float_dtype, copy=False)
|
||||
|
||||
mask = mask.astype(bool, copy=False)
|
||||
if not multichannel:
|
||||
image = image[..., np.newaxis]
|
||||
out = np.copy(image, order='C')
|
||||
|
||||
# Create biharmonic coefficients ndarray
|
||||
radius = 2
|
||||
coef_shape = (2 * radius + 1,) * mask.ndim
|
||||
coef_center = (radius,) * mask.ndim
|
||||
neigh_coef_full, coef_idx, coef_vals = _get_neigh_coef(coef_shape,
|
||||
coef_center,
|
||||
dtype=out.dtype)
|
||||
|
||||
# stride for the last spatial dimension
|
||||
channel_stride_bytes = out.strides[-2]
|
||||
|
||||
# offsets to all neighboring non-zero elements in the footprint
|
||||
offsets = coef_idx - radius
|
||||
|
||||
# determine per-channel intensity limits
|
||||
known_points = image[~mask]
|
||||
limits = (known_points.min(axis=0), known_points.max(axis=0))
|
||||
|
||||
if split_into_regions:
|
||||
# Split inpainting mask into independent regions
|
||||
kernel = ndi.generate_binary_structure(mask.ndim, 1)
|
||||
mask_dilated = ndi.binary_dilation(mask, structure=kernel)
|
||||
mask_labeled = label(mask_dilated)
|
||||
mask_labeled *= mask
|
||||
|
||||
bbox_slices = ndi.find_objects(mask_labeled)
|
||||
|
||||
for idx_region, bb_slice in enumerate(bbox_slices, 1):
|
||||
# expand object bounding boxes by the biharmonic kernel radius
|
||||
roi_sl = tuple(slice(max(sl.start - radius, 0),
|
||||
min(sl.stop + radius, size))
|
||||
for sl, size in zip(bb_slice, mask_labeled.shape))
|
||||
# extract only the region surrounding the label of interest
|
||||
mask_region = mask_labeled[roi_sl] == idx_region
|
||||
# add slice for axes
|
||||
roi_sl += (slice(None), )
|
||||
# copy for contiguity and to account for possible ROI overlap
|
||||
otmp = out[roi_sl].copy()
|
||||
|
||||
# compute raveled offsets for the ROI
|
||||
ostrides = np.array([s // channel_stride_bytes
|
||||
for s in otmp[..., 0].strides])
|
||||
raveled_offsets = np.sum(offsets * ostrides[..., np.newaxis],
|
||||
axis=0)
|
||||
|
||||
_inpaint_biharmonic_single_region(
|
||||
image[roi_sl], mask_region, otmp,
|
||||
neigh_coef_full, coef_vals, raveled_offsets
|
||||
)
|
||||
# assign output to the
|
||||
out[roi_sl] = otmp
|
||||
else:
|
||||
# compute raveled offsets for output image
|
||||
ostrides = np.array([s // channel_stride_bytes
|
||||
for s in out[..., 0].strides])
|
||||
raveled_offsets = np.sum(offsets * ostrides[..., np.newaxis], axis=0)
|
||||
|
||||
_inpaint_biharmonic_single_region(
|
||||
image, mask, out, neigh_coef_full, coef_vals, raveled_offsets
|
||||
)
|
||||
|
||||
# Handle enormous values on a per-channel basis
|
||||
np.clip(out, a_min=limits[0], a_max=limits[1], out=out)
|
||||
|
||||
if not multichannel:
|
||||
out = out[..., 0]
|
||||
|
||||
return out
|
||||
@@ -1,319 +0,0 @@
|
||||
import itertools
|
||||
import functools
|
||||
|
||||
import numpy as np
|
||||
from scipy import ndimage as ndi
|
||||
|
||||
from .._shared.utils import _supported_float_type
|
||||
from ..metrics import mean_squared_error
|
||||
from ..util import img_as_float
|
||||
|
||||
|
||||
|
||||
def _interpolate_image(image, *, multichannel=False):
|
||||
"""Replacing each pixel in ``image`` with the average of its neighbors.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Input data to be interpolated.
|
||||
multichannel : bool, optional
|
||||
Whether the last axis of the image is to be interpreted as multiple
|
||||
channels or another spatial dimension.
|
||||
|
||||
Returns
|
||||
-------
|
||||
interp : ndarray
|
||||
Interpolated version of `image`.
|
||||
"""
|
||||
spatialdims = image.ndim if not multichannel else image.ndim - 1
|
||||
conv_filter = ndi.generate_binary_structure(spatialdims, 1).astype(image.dtype)
|
||||
conv_filter.ravel()[conv_filter.size // 2] = 0
|
||||
conv_filter /= conv_filter.sum()
|
||||
|
||||
if multichannel:
|
||||
interp = np.zeros_like(image)
|
||||
for i in range(image.shape[-1]):
|
||||
interp[..., i] = ndi.convolve(image[..., i], conv_filter,
|
||||
mode='mirror')
|
||||
else:
|
||||
interp = ndi.convolve(image, conv_filter, mode='mirror')
|
||||
return interp
|
||||
|
||||
|
||||
def _generate_grid_slice(shape, *, offset, stride=3):
|
||||
"""Generate slices of uniformly-spaced points in an array.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
shape : tuple of int
|
||||
Shape of the mask.
|
||||
offset : int
|
||||
The offset of the grid of ones. Iterating over ``offset`` will cover
|
||||
the entire array. It should be between 0 and ``stride ** ndim``, not
|
||||
inclusive, where ``ndim = len(shape)``.
|
||||
stride : int, optional
|
||||
The spacing between ones, used in each dimension.
|
||||
|
||||
Returns
|
||||
-------
|
||||
mask : ndarray
|
||||
The mask.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> shape = (4, 4)
|
||||
>>> array = np.zeros(shape, dtype=int)
|
||||
>>> grid_slice = _generate_grid_slice(shape, offset=0, stride=2)
|
||||
>>> array[grid_slice] = 1
|
||||
>>> print(array)
|
||||
[[1 0 1 0]
|
||||
[0 0 0 0]
|
||||
[1 0 1 0]
|
||||
[0 0 0 0]]
|
||||
|
||||
Changing the offset moves the location of the 1s:
|
||||
|
||||
>>> array = np.zeros(shape, dtype=int)
|
||||
>>> grid_slice = _generate_grid_slice(shape, offset=3, stride=2)
|
||||
>>> array[grid_slice] = 1
|
||||
>>> print(array)
|
||||
[[0 0 0 0]
|
||||
[0 1 0 1]
|
||||
[0 0 0 0]
|
||||
[0 1 0 1]]
|
||||
"""
|
||||
phases = np.unravel_index(offset, (stride,) * len(shape))
|
||||
mask = tuple(slice(p, None, stride) for p in phases)
|
||||
|
||||
return mask
|
||||
|
||||
|
||||
def _invariant_denoise(image, denoise_function, *, stride=4,
|
||||
masks=None, denoiser_kwargs=None):
|
||||
"""Apply a J-invariant version of `denoise_function`.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Input data to be denoised (converted using `img_as_float`).
|
||||
denoise_function : function
|
||||
Original denoising function.
|
||||
stride : int, optional
|
||||
Stride used in masking procedure that converts `denoise_function`
|
||||
to J-invariance.
|
||||
masks : list of ndarray, optional
|
||||
Set of masks to use for computing J-invariant output. If `None`,
|
||||
a full set of masks covering the image will be used.
|
||||
denoiser_kwargs:
|
||||
Keyword arguments passed to `denoise_function`.
|
||||
|
||||
Returns
|
||||
-------
|
||||
output : ndarray
|
||||
Denoised image, of same shape as `image`.
|
||||
"""
|
||||
image = img_as_float(image)
|
||||
|
||||
# promote float16->float32 if needed
|
||||
float_dtype = _supported_float_type(image.dtype)
|
||||
image = image.astype(float_dtype, copy=False)
|
||||
|
||||
if denoiser_kwargs is None:
|
||||
denoiser_kwargs = {}
|
||||
|
||||
multichannel = denoiser_kwargs.get('channel_axis', None) is not None
|
||||
interp = _interpolate_image(image, multichannel=multichannel)
|
||||
output = np.zeros_like(image)
|
||||
|
||||
if masks is None:
|
||||
spatialdims = image.ndim if not multichannel else image.ndim - 1
|
||||
n_masks = stride ** spatialdims
|
||||
masks = (_generate_grid_slice(image.shape[:spatialdims],
|
||||
offset=idx, stride=stride)
|
||||
for idx in range(n_masks))
|
||||
|
||||
for mask in masks:
|
||||
input_image = image.copy()
|
||||
input_image[mask] = interp[mask]
|
||||
output[mask] = denoise_function(input_image, **denoiser_kwargs)[mask]
|
||||
return output
|
||||
|
||||
|
||||
def _product_from_dict(dictionary):
|
||||
"""Utility function to convert parameter ranges to parameter combinations.
|
||||
|
||||
Converts a dict of lists into a list of dicts whose values consist of the
|
||||
cartesian product of the values in the original dict.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
dictionary : dict of lists
|
||||
Dictionary of lists to be multiplied.
|
||||
|
||||
Yields
|
||||
------
|
||||
selections : dicts of values
|
||||
Dicts containing individual combinations of the values in the input
|
||||
dict.
|
||||
"""
|
||||
keys = dictionary.keys()
|
||||
for element in itertools.product(*dictionary.values()):
|
||||
yield dict(zip(keys, element))
|
||||
|
||||
|
||||
def calibrate_denoiser(image, denoise_function, denoise_parameters, *,
|
||||
stride=4, approximate_loss=True,
|
||||
extra_output=False):
|
||||
"""Calibrate a denoising function and return optimal J-invariant version.
|
||||
|
||||
The returned function is partially evaluated with optimal parameter values
|
||||
set for denoising the input image.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Input data to be denoised (converted using `img_as_float`).
|
||||
denoise_function : function
|
||||
Denoising function to be calibrated.
|
||||
denoise_parameters : dict of list
|
||||
Ranges of parameters for `denoise_function` to be calibrated over.
|
||||
stride : int, optional
|
||||
Stride used in masking procedure that converts `denoise_function`
|
||||
to J-invariance.
|
||||
approximate_loss : bool, optional
|
||||
Whether to approximate the self-supervised loss used to evaluate the
|
||||
denoiser by only computing it on one masked version of the image.
|
||||
If False, the runtime will be a factor of `stride**image.ndim` longer.
|
||||
extra_output : bool, optional
|
||||
If True, return parameters and losses in addition to the calibrated
|
||||
denoising function
|
||||
|
||||
Returns
|
||||
-------
|
||||
best_denoise_function : function
|
||||
The optimal J-invariant version of `denoise_function`.
|
||||
|
||||
If `extra_output` is True, the following tuple is also returned:
|
||||
|
||||
(parameters_tested, losses) : tuple (list of dict, list of int)
|
||||
List of parameters tested for `denoise_function`, as a dictionary of
|
||||
kwargs
|
||||
Self-supervised loss for each set of parameters in `parameters_tested`.
|
||||
|
||||
|
||||
Notes
|
||||
-----
|
||||
|
||||
The calibration procedure uses a self-supervised mean-square-error loss
|
||||
to evaluate the performance of J-invariant versions of `denoise_function`.
|
||||
The minimizer of the self-supervised loss is also the minimizer of the
|
||||
ground-truth loss (i.e., the true MSE error) [1]. The returned function
|
||||
can be used on the original noisy image, or other images with similar
|
||||
characteristics.
|
||||
|
||||
Increasing the stride increases the performance of `best_denoise_function`
|
||||
at the expense of increasing its runtime. It has no effect on the runtime
|
||||
of the calibration.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] J. Batson & L. Royer. Noise2Self: Blind Denoising by Self-Supervision,
|
||||
International Conference on Machine Learning, p. 524-533 (2019).
|
||||
|
||||
Examples
|
||||
--------
|
||||
|
||||
>>> from skimage import color, data
|
||||
>>> from skimage.restoration import denoise_wavelet
|
||||
>>> import numpy as np
|
||||
>>> img = color.rgb2gray(data.astronaut()[:50, :50])
|
||||
>>> rng = np.random.default_rng()
|
||||
>>> noisy = img + 0.5 * img.std() * rng.standard_normal(img.shape)
|
||||
>>> parameters = {'sigma': np.arange(0.1, 0.4, 0.02)}
|
||||
>>> denoising_function = calibrate_denoiser(noisy, denoise_wavelet,
|
||||
... denoise_parameters=parameters)
|
||||
>>> denoised_img = denoising_function(img)
|
||||
|
||||
"""
|
||||
parameters_tested, losses = _calibrate_denoiser_search(
|
||||
image, denoise_function,
|
||||
denoise_parameters=denoise_parameters,
|
||||
stride=stride,
|
||||
approximate_loss=approximate_loss
|
||||
)
|
||||
|
||||
idx = np.argmin(losses)
|
||||
best_parameters = parameters_tested[idx]
|
||||
|
||||
best_denoise_function = functools.partial(
|
||||
_invariant_denoise,
|
||||
denoise_function=denoise_function,
|
||||
stride=stride,
|
||||
denoiser_kwargs=best_parameters,
|
||||
)
|
||||
|
||||
if extra_output:
|
||||
return best_denoise_function, (parameters_tested, losses)
|
||||
else:
|
||||
return best_denoise_function
|
||||
|
||||
|
||||
def _calibrate_denoiser_search(image, denoise_function, denoise_parameters, *,
|
||||
stride=4, approximate_loss=True):
|
||||
"""Return a parameter search history with losses for a denoise function.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Input data to be denoised (converted using `img_as_float`).
|
||||
denoise_function : function
|
||||
Denoising function to be calibrated.
|
||||
denoise_parameters : dict of list
|
||||
Ranges of parameters for `denoise_function` to be calibrated over.
|
||||
stride : int, optional
|
||||
Stride used in masking procedure that converts `denoise_function`
|
||||
to J-invariance.
|
||||
approximate_loss : bool, optional
|
||||
Whether to approximate the self-supervised loss used to evaluate the
|
||||
denoiser by only computing it on one masked version of the image.
|
||||
If False, the runtime will be a factor of `stride**image.ndim` longer.
|
||||
|
||||
Returns
|
||||
-------
|
||||
parameters_tested : list of dict
|
||||
List of parameters tested for `denoise_function`, as a dictionary of
|
||||
kwargs.
|
||||
losses : list of int
|
||||
Self-supervised loss for each set of parameters in `parameters_tested`.
|
||||
"""
|
||||
image = img_as_float(image)
|
||||
parameters_tested = list(_product_from_dict(denoise_parameters))
|
||||
losses = []
|
||||
|
||||
for denoiser_kwargs in parameters_tested:
|
||||
multichannel = denoiser_kwargs.get('channel_axis', None) is not None
|
||||
if not approximate_loss:
|
||||
denoised = _invariant_denoise(
|
||||
image, denoise_function,
|
||||
stride=stride,
|
||||
denoiser_kwargs=denoiser_kwargs
|
||||
)
|
||||
loss = mean_squared_error(image, denoised)
|
||||
else:
|
||||
spatialdims = image.ndim if not multichannel else image.ndim - 1
|
||||
n_masks = stride ** spatialdims
|
||||
mask = _generate_grid_slice(image.shape[:spatialdims],
|
||||
offset=n_masks // 2, stride=stride)
|
||||
|
||||
masked_denoised = _invariant_denoise(
|
||||
image, denoise_function,
|
||||
masks=[mask],
|
||||
denoiser_kwargs=denoiser_kwargs
|
||||
)
|
||||
|
||||
loss = mean_squared_error(image[mask], masked_denoised[mask])
|
||||
|
||||
losses.append(loss)
|
||||
|
||||
return parameters_tested, losses
|
||||
@@ -1,178 +0,0 @@
|
||||
import numpy as np
|
||||
|
||||
from .._shared import utils
|
||||
from .._shared.utils import convert_to_float
|
||||
from ._nl_means_denoising import (_nl_means_denoising_2d,
|
||||
_nl_means_denoising_3d,
|
||||
_fast_nl_means_denoising_2d,
|
||||
_fast_nl_means_denoising_3d,
|
||||
_fast_nl_means_denoising_4d)
|
||||
|
||||
|
||||
@utils.channel_as_last_axis()
|
||||
def denoise_nl_means(image, patch_size=7, patch_distance=11, h=0.1,
|
||||
fast_mode=True, sigma=0., *,
|
||||
preserve_range=False, channel_axis=None):
|
||||
"""Perform non-local means denoising on 2D-4D grayscale or RGB images.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : 2D or 3D ndarray
|
||||
Input image to be denoised, which can be 2D or 3D, and grayscale
|
||||
or RGB (for 2D images only, see ``channel_axis`` parameter). There can
|
||||
be any number of channels (does not strictly have to be RGB).
|
||||
patch_size : int, optional
|
||||
Size of patches used for denoising.
|
||||
patch_distance : int, optional
|
||||
Maximal distance in pixels where to search patches used for denoising.
|
||||
h : float, optional
|
||||
Cut-off distance (in gray levels). The higher h, the more permissive
|
||||
one is in accepting patches. A higher h results in a smoother image,
|
||||
at the expense of blurring features. For a Gaussian noise of standard
|
||||
deviation sigma, a rule of thumb is to choose the value of h to be
|
||||
sigma of slightly less.
|
||||
fast_mode : bool, optional
|
||||
If True (default value), a fast version of the non-local means
|
||||
algorithm is used. If False, the original version of non-local means is
|
||||
used. See the Notes section for more details about the algorithms.
|
||||
sigma : float, optional
|
||||
The standard deviation of the (Gaussian) noise. If provided, a more
|
||||
robust computation of patch weights is computed that takes the expected
|
||||
noise variance into account (see Notes below).
|
||||
preserve_range : bool, optional
|
||||
Whether to keep the original range of values. Otherwise, the input
|
||||
image is converted according to the conventions of `img_as_float`.
|
||||
Also see https://scikit-image.org/docs/dev/user_guide/data_types.html
|
||||
channel_axis : int or None, optional
|
||||
If None, the image is assumed to be a grayscale (single channel) image.
|
||||
Otherwise, this parameter indicates which axis of the array corresponds
|
||||
to channels.
|
||||
|
||||
.. versionadded:: 0.19
|
||||
``channel_axis`` was added in 0.19.
|
||||
|
||||
Returns
|
||||
-------
|
||||
result : ndarray
|
||||
Denoised image, of same shape as `image`.
|
||||
|
||||
Notes
|
||||
-----
|
||||
|
||||
The non-local means algorithm is well suited for denoising images with
|
||||
specific textures. The principle of the algorithm is to average the value
|
||||
of a given pixel with values of other pixels in a limited neighborhood,
|
||||
provided that the *patches* centered on the other pixels are similar enough
|
||||
to the patch centered on the pixel of interest.
|
||||
|
||||
In the original version of the algorithm [1]_, corresponding to
|
||||
``fast=False``, the computational complexity is::
|
||||
|
||||
image.size * patch_size ** image.ndim * patch_distance ** image.ndim
|
||||
|
||||
Hence, changing the size of patches or their maximal distance has a
|
||||
strong effect on computing times, especially for 3-D images.
|
||||
|
||||
However, the default behavior corresponds to ``fast_mode=True``, for which
|
||||
another version of non-local means [2]_ is used, corresponding to a
|
||||
complexity of::
|
||||
|
||||
image.size * patch_distance ** image.ndim
|
||||
|
||||
The computing time depends only weakly on the patch size, thanks to
|
||||
the computation of the integral of patches distances for a given
|
||||
shift, that reduces the number of operations [1]_. Therefore, this
|
||||
algorithm executes faster than the classic algorithm
|
||||
(``fast_mode=False``), at the expense of using twice as much memory.
|
||||
This implementation has been proven to be more efficient compared to
|
||||
other alternatives, see e.g. [3]_.
|
||||
|
||||
Compared to the classic algorithm, all pixels of a patch contribute
|
||||
to the distance to another patch with the same weight, no matter
|
||||
their distance to the center of the patch. This coarser computation
|
||||
of the distance can result in a slightly poorer denoising
|
||||
performance. Moreover, for small images (images with a linear size
|
||||
that is only a few times the patch size), the classic algorithm can
|
||||
be faster due to boundary effects.
|
||||
|
||||
The image is padded using the `reflect` mode of `skimage.util.pad`
|
||||
before denoising.
|
||||
|
||||
If the noise standard deviation, `sigma`, is provided a more robust
|
||||
computation of patch weights is used. Subtracting the known noise variance
|
||||
from the computed patch distances improves the estimates of patch
|
||||
similarity, giving a moderate improvement to denoising performance [4]_.
|
||||
It was also mentioned as an option for the fast variant of the algorithm in
|
||||
[3]_.
|
||||
|
||||
When `sigma` is provided, a smaller `h` should typically be used to
|
||||
avoid oversmoothing. The optimal value for `h` depends on the image
|
||||
content and noise level, but a reasonable starting point is
|
||||
``h = 0.8 * sigma`` when `fast_mode` is `True`, or ``h = 0.6 * sigma`` when
|
||||
`fast_mode` is `False`.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] A. Buades, B. Coll, & J-M. Morel. A non-local algorithm for image
|
||||
denoising. In CVPR 2005, Vol. 2, pp. 60-65, IEEE.
|
||||
:DOI:`10.1109/CVPR.2005.38`
|
||||
|
||||
.. [2] J. Darbon, A. Cunha, T.F. Chan, S. Osher, and G.J. Jensen, Fast
|
||||
nonlocal filtering applied to electron cryomicroscopy, in 5th IEEE
|
||||
International Symposium on Biomedical Imaging: From Nano to Macro,
|
||||
2008, pp. 1331-1334.
|
||||
:DOI:`10.1109/ISBI.2008.4541250`
|
||||
|
||||
.. [3] Jacques Froment. Parameter-Free Fast Pixelwise Non-Local Means
|
||||
Denoising. Image Processing On Line, 2014, vol. 4, pp. 300-326.
|
||||
:DOI:`10.5201/ipol.2014.120`
|
||||
|
||||
.. [4] A. Buades, B. Coll, & J-M. Morel. Non-Local Means Denoising.
|
||||
Image Processing On Line, 2011, vol. 1, pp. 208-212.
|
||||
:DOI:`10.5201/ipol.2011.bcm_nlm`
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> a = np.zeros((40, 40))
|
||||
>>> a[10:-10, 10:-10] = 1.
|
||||
>>> rng = np.random.default_rng()
|
||||
>>> a += 0.3 * rng.standard_normal(a.shape)
|
||||
>>> denoised_a = denoise_nl_means(a, 7, 5, 0.1)
|
||||
"""
|
||||
if channel_axis is None:
|
||||
multichannel = False
|
||||
image = image[..., np.newaxis]
|
||||
else:
|
||||
multichannel = True
|
||||
|
||||
ndim_no_channel = image.ndim - 1
|
||||
if (ndim_no_channel < 2) or (ndim_no_channel > 4):
|
||||
raise NotImplementedError(
|
||||
"Non-local means denoising is only implemented for 2D, "
|
||||
"3D or 4D grayscale or multichannel images.")
|
||||
|
||||
image = convert_to_float(image, preserve_range)
|
||||
if not image.flags.c_contiguous:
|
||||
image = np.ascontiguousarray(image)
|
||||
|
||||
kwargs = dict(s=patch_size, d=patch_distance, h=h, var=sigma * sigma)
|
||||
if ndim_no_channel == 2:
|
||||
nlm_func = (_fast_nl_means_denoising_2d if fast_mode else
|
||||
_nl_means_denoising_2d)
|
||||
elif ndim_no_channel == 3:
|
||||
if multichannel and not fast_mode:
|
||||
raise NotImplementedError(
|
||||
"Multichannel 3D requires fast_mode to be True.")
|
||||
if fast_mode:
|
||||
nlm_func = _fast_nl_means_denoising_3d
|
||||
else:
|
||||
# have to drop the size 1 channel axis for slow mode
|
||||
image = image[..., 0]
|
||||
nlm_func = _nl_means_denoising_3d
|
||||
elif ndim_no_channel == 4:
|
||||
if fast_mode:
|
||||
nlm_func = _fast_nl_means_denoising_4d
|
||||
else:
|
||||
raise NotImplementedError("4D requires fast_mode to be True.")
|
||||
dn = np.asarray(nlm_func(image, **kwargs))
|
||||
return dn
|
||||
@@ -1,195 +0,0 @@
|
||||
import numpy as np
|
||||
|
||||
from .._shared.utils import _supported_float_type
|
||||
from ._rolling_ball_cy import apply_kernel, apply_kernel_nan
|
||||
|
||||
|
||||
def rolling_ball(image, *, radius=100, kernel=None,
|
||||
nansafe=False, num_threads=None):
|
||||
"""Estimate background intensity by rolling/translating a kernel.
|
||||
|
||||
This rolling ball algorithm estimates background intensity for a
|
||||
ndimage in case of uneven exposure. It is a generalization of the
|
||||
frequently used rolling ball algorithm [1]_.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
The image to be filtered.
|
||||
radius : int, optional
|
||||
Radius of a ball shaped kernel to be rolled/translated in the image.
|
||||
Used if ``kernel = None``.
|
||||
kernel : ndarray, optional
|
||||
The kernel to be rolled/translated in the image. It must have the
|
||||
same number of dimensions as ``image``. Kernel is filled with the
|
||||
intensity of the kernel at that position.
|
||||
nansafe: bool, optional
|
||||
If ``False`` (default) assumes that none of the values in ``image``
|
||||
are ``np.nan``, and uses a faster implementation.
|
||||
num_threads: int, optional
|
||||
The maximum number of threads to use. If ``None`` use the OpenMP
|
||||
default value; typically equal to the maximum number of virtual cores.
|
||||
Note: This is an upper limit to the number of threads. The exact number
|
||||
is determined by the system's OpenMP library.
|
||||
|
||||
Returns
|
||||
-------
|
||||
background : ndarray
|
||||
The estimated background of the image.
|
||||
|
||||
Notes
|
||||
-----
|
||||
For the pixel that has its background intensity estimated (without loss
|
||||
of generality at ``center``) the rolling ball method centers ``kernel``
|
||||
under it and raises the kernel until the surface touches the image umbra
|
||||
at some ``pos=(y,x)``. The background intensity is then estimated
|
||||
using the image intensity at that position (``image[pos]``) plus the
|
||||
difference of ``kernel[center] - kernel[pos]``.
|
||||
|
||||
This algorithm assumes that dark pixels correspond to the background. If
|
||||
you have a bright background, invert the image before passing it to the
|
||||
function, e.g., using `utils.invert`. See the gallery example for details.
|
||||
|
||||
This algorithm is sensitive to noise (in particular salt-and-pepper
|
||||
noise). If this is a problem in your image, you can apply mild
|
||||
gaussian smoothing before passing the image to this function.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] Sternberg, Stanley R. "Biomedical image processing." Computer 1
|
||||
(1983): 22-34. :DOI:`10.1109/MC.1983.1654163`
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> import numpy as np
|
||||
>>> from skimage import data
|
||||
>>> from skimage.restoration import rolling_ball
|
||||
>>> image = data.coins()
|
||||
>>> background = rolling_ball(data.coins())
|
||||
>>> filtered_image = image - background
|
||||
|
||||
|
||||
>>> import numpy as np
|
||||
>>> from skimage import data
|
||||
>>> from skimage.restoration import rolling_ball, ellipsoid_kernel
|
||||
>>> image = data.coins()
|
||||
>>> kernel = ellipsoid_kernel((101, 101), 75)
|
||||
>>> background = rolling_ball(data.coins(), kernel=kernel)
|
||||
>>> filtered_image = image - background
|
||||
"""
|
||||
|
||||
image = np.asarray(image)
|
||||
float_type = _supported_float_type(image.dtype)
|
||||
img = image.astype(float_type, copy=False)
|
||||
|
||||
if num_threads is None:
|
||||
num_threads = 0
|
||||
|
||||
if kernel is None:
|
||||
kernel = ball_kernel(radius, image.ndim)
|
||||
|
||||
kernel = kernel.astype(float_type)
|
||||
kernel_shape = np.asarray(kernel.shape)
|
||||
kernel_center = (kernel_shape // 2)
|
||||
center_intensity = kernel[tuple(kernel_center)]
|
||||
|
||||
intensity_difference = center_intensity - kernel
|
||||
intensity_difference[kernel == np.Inf] = np.Inf
|
||||
intensity_difference = intensity_difference.astype(img.dtype)
|
||||
intensity_difference = intensity_difference.reshape(-1)
|
||||
|
||||
img = np.pad(img, kernel_center[:, np.newaxis],
|
||||
constant_values=np.Inf, mode="constant")
|
||||
|
||||
func = apply_kernel_nan if nansafe else apply_kernel
|
||||
background = func(
|
||||
img.reshape(-1),
|
||||
intensity_difference,
|
||||
np.zeros_like(image, dtype=img.dtype).reshape(-1),
|
||||
np.array(image.shape, dtype=np.intp),
|
||||
np.array(img.shape, dtype=np.intp),
|
||||
kernel_shape.astype(np.intp),
|
||||
num_threads
|
||||
)
|
||||
|
||||
background = background.astype(image.dtype, copy=False)
|
||||
|
||||
return background
|
||||
|
||||
|
||||
def ball_kernel(radius, ndim):
|
||||
"""Create a ball kernel for restoration.rolling_ball.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
radius : int
|
||||
Radius of the ball.
|
||||
ndim : int
|
||||
Number of dimensions of the ball. ``ndim`` should match the
|
||||
dimensionality of the image the kernel will be applied to.
|
||||
|
||||
Returns
|
||||
-------
|
||||
kernel : ndarray
|
||||
The kernel containing the surface intensity of the top half
|
||||
of the ellipsoid.
|
||||
|
||||
See Also
|
||||
--------
|
||||
rolling_ball
|
||||
"""
|
||||
|
||||
kernel_coords = np.stack(
|
||||
np.meshgrid(
|
||||
*[np.arange(-x, x + 1) for x in [np.ceil(radius)] * ndim],
|
||||
indexing='ij'
|
||||
),
|
||||
axis=-1
|
||||
)
|
||||
|
||||
sum_of_squares = np.sum(kernel_coords ** 2, axis=-1)
|
||||
distance_from_center = np.sqrt(sum_of_squares)
|
||||
kernel = np.sqrt(np.clip(radius ** 2 - sum_of_squares, 0, None))
|
||||
kernel[distance_from_center > radius] = np.Inf
|
||||
|
||||
return kernel
|
||||
|
||||
|
||||
def ellipsoid_kernel(shape, intensity):
|
||||
"""Create an ellipoid kernel for restoration.rolling_ball.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
shape : arraylike
|
||||
Length of the principal axis of the ellipsoid (excluding
|
||||
the intensity axis). The kernel needs to have the same
|
||||
dimensionality as the image it will be applied to.
|
||||
intensity : int
|
||||
Length of the intensity axis of the ellipsoid.
|
||||
|
||||
Returns
|
||||
-------
|
||||
kernel : ndarray
|
||||
The kernel containing the surface intensity of the top half
|
||||
of the ellipsoid.
|
||||
|
||||
See Also
|
||||
--------
|
||||
rolling_ball
|
||||
"""
|
||||
|
||||
shape = np.asarray(shape)
|
||||
semi_axis = np.clip(shape // 2, 1, None)
|
||||
|
||||
kernel_coords = np.stack(
|
||||
np.meshgrid(
|
||||
*[np.arange(-x, x + 1) for x in semi_axis],
|
||||
indexing='ij'
|
||||
),
|
||||
axis=-1)
|
||||
|
||||
intensity_scaling = 1 - np.sum((kernel_coords / semi_axis) ** 2, axis=-1)
|
||||
kernel = intensity * np.sqrt(np.clip(intensity_scaling, 0, None))
|
||||
kernel[intensity_scaling < 0] = np.Inf
|
||||
|
||||
return kernel
|
||||
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
File diff suppressed because it is too large
Load Diff
@@ -1,173 +0,0 @@
|
||||
import numpy as np
|
||||
|
||||
from skimage import data, img_as_float
|
||||
from skimage._shared import testing
|
||||
from skimage._shared.testing import assert_allclose
|
||||
from skimage._shared.utils import _supported_float_type
|
||||
from skimage.color import rgb2gray
|
||||
from skimage.metrics import mean_squared_error, normalized_root_mse
|
||||
from skimage.morphology import binary_dilation, disk
|
||||
from skimage.restoration import inpaint
|
||||
|
||||
|
||||
@testing.parametrize('dtype', [np.float16, np.float32, np.float64])
|
||||
@testing.parametrize('split_into_regions', [False, True])
|
||||
def test_inpaint_biharmonic_2d(dtype, split_into_regions):
|
||||
img = np.tile(np.square(np.linspace(0, 1, 5, dtype=dtype)), (5, 1))
|
||||
mask = np.zeros_like(img)
|
||||
mask[2, 2:] = 1
|
||||
mask[1, 3:] = 1
|
||||
mask[0, 4:] = 1
|
||||
img[np.where(mask)] = 0
|
||||
out = inpaint.inpaint_biharmonic(img, mask,
|
||||
split_into_regions=split_into_regions)
|
||||
assert out.dtype == _supported_float_type(img)
|
||||
|
||||
ref = np.array(
|
||||
[[0., 0.0625, 0.25000000, 0.5625000, 0.73925058],
|
||||
[0., 0.0625, 0.25000000, 0.5478048, 0.76557821],
|
||||
[0., 0.0625, 0.25842878, 0.5623079, 0.85927796],
|
||||
[0., 0.0625, 0.25000000, 0.5625000, 1.00000000],
|
||||
[0., 0.0625, 0.25000000, 0.5625000, 1.00000000]]
|
||||
)
|
||||
rtol = 1e-7 if dtype == np.float64 else 1e-6
|
||||
assert_allclose(ref, out, rtol=rtol)
|
||||
|
||||
|
||||
@testing.parametrize('channel_axis', [0, 1, -1])
|
||||
def test_inpaint_biharmonic_2d_color(channel_axis):
|
||||
img = img_as_float(data.astronaut()[:64, :64])
|
||||
|
||||
mask = np.zeros(img.shape[:2], dtype=bool)
|
||||
mask[8:16, :16] = 1
|
||||
img_defect = img * ~mask[..., np.newaxis]
|
||||
mse_defect = mean_squared_error(img, img_defect)
|
||||
|
||||
img_defect = np.moveaxis(img_defect, -1, channel_axis)
|
||||
img_restored = inpaint.inpaint_biharmonic(img_defect, mask,
|
||||
channel_axis=channel_axis)
|
||||
img_restored = np.moveaxis(img_restored, channel_axis, -1)
|
||||
mse_restored = mean_squared_error(img, img_restored)
|
||||
|
||||
assert mse_restored < 0.01 * mse_defect
|
||||
|
||||
|
||||
@testing.parametrize('dtype', [np.float32, np.float64])
|
||||
def test_inpaint_biharmonic_2d_float_dtypes(dtype):
|
||||
img = np.tile(np.square(np.linspace(0, 1, 5)), (5, 1))
|
||||
mask = np.zeros_like(img)
|
||||
mask[2, 2:] = 1
|
||||
mask[1, 3:] = 1
|
||||
mask[0, 4:] = 1
|
||||
img[np.where(mask)] = 0
|
||||
img = img.astype(dtype, copy=False)
|
||||
out = inpaint.inpaint_biharmonic(img, mask)
|
||||
assert out.dtype == img.dtype
|
||||
ref = np.array(
|
||||
[[0., 0.0625, 0.25000000, 0.5625000, 0.73925058],
|
||||
[0., 0.0625, 0.25000000, 0.5478048, 0.76557821],
|
||||
[0., 0.0625, 0.25842878, 0.5623079, 0.85927796],
|
||||
[0., 0.0625, 0.25000000, 0.5625000, 1.00000000],
|
||||
[0., 0.0625, 0.25000000, 0.5625000, 1.00000000]]
|
||||
)
|
||||
assert_allclose(ref, out, rtol=1e-5)
|
||||
|
||||
|
||||
@testing.parametrize('split_into_regions', [False, True])
|
||||
def test_inpaint_biharmonic_3d(split_into_regions):
|
||||
img = np.tile(np.square(np.linspace(0, 1, 5)), (5, 1))
|
||||
img = np.dstack((img, img.T))
|
||||
mask = np.zeros_like(img)
|
||||
mask[2, 2:, :] = 1
|
||||
mask[1, 3:, :] = 1
|
||||
mask[0, 4:, :] = 1
|
||||
img[np.where(mask)] = 0
|
||||
out = inpaint.inpaint_biharmonic(img, mask,
|
||||
split_into_regions=split_into_regions)
|
||||
ref = np.dstack((
|
||||
np.array(
|
||||
[[0.0000, 0.0625, 0.25000000, 0.56250000, 0.53752796],
|
||||
[0.0000, 0.0625, 0.25000000, 0.44443780, 0.53762210],
|
||||
[0.0000, 0.0625, 0.23693666, 0.46621112, 0.68615592],
|
||||
[0.0000, 0.0625, 0.25000000, 0.56250000, 1.00000000],
|
||||
[0.0000, 0.0625, 0.25000000, 0.56250000, 1.00000000]]),
|
||||
np.array(
|
||||
[[0.0000, 0.0000, 0.00000000, 0.00000000, 0.19621902],
|
||||
[0.0625, 0.0625, 0.06250000, 0.17470756, 0.30140091],
|
||||
[0.2500, 0.2500, 0.27241289, 0.35155440, 0.43068654],
|
||||
[0.5625, 0.5625, 0.56250000, 0.56250000, 0.56250000],
|
||||
[1.0000, 1.0000, 1.00000000, 1.00000000, 1.00000000]])
|
||||
))
|
||||
assert_allclose(ref, out)
|
||||
|
||||
|
||||
def test_invalid_input():
|
||||
img, mask = np.zeros([]), np.zeros([])
|
||||
with testing.raises(ValueError):
|
||||
inpaint.inpaint_biharmonic(img, mask)
|
||||
|
||||
img, mask = np.zeros((2, 2)), np.zeros((4, 1))
|
||||
with testing.raises(ValueError):
|
||||
inpaint.inpaint_biharmonic(img, mask)
|
||||
|
||||
img = np.ma.array(np.zeros((2, 2)), mask=[[0, 0], [0, 0]])
|
||||
mask = np.zeros((2, 2))
|
||||
with testing.raises(TypeError):
|
||||
inpaint.inpaint_biharmonic(img, mask)
|
||||
|
||||
|
||||
@testing.parametrize('dtype', [np.uint8, np.float32, np.float64])
|
||||
@testing.parametrize('order', ['C', 'F'])
|
||||
@testing.parametrize('channel_axis', [None, -1])
|
||||
@testing.parametrize('split_into_regions', [False, True])
|
||||
def test_inpaint_nrmse(dtype, order, channel_axis, split_into_regions):
|
||||
image_orig = data.astronaut()[:, :200]
|
||||
float_dtype = np.float32 if dtype == np.float32 else np.float64
|
||||
image_orig = image_orig.astype(float_dtype, copy=False)
|
||||
|
||||
# Create mask with six block defect regions
|
||||
mask = np.zeros(image_orig.shape[:-1], dtype=bool)
|
||||
mask[20:50, 3:20] = 1
|
||||
mask[165:180, 90:155] = 1
|
||||
mask[40:60, 170:195] = 1
|
||||
mask[-60:-40, 170:195] = 1
|
||||
mask[-180:-165, 90:155] = 1
|
||||
mask[-50:-20, :20] = 1
|
||||
|
||||
# add a few long, narrow defects
|
||||
mask[200:205, -200:] = 1
|
||||
mask[150:255, 20:22] = 1
|
||||
mask[365:368, 60:130] = 1
|
||||
|
||||
# add randomly positioned small point-like defects
|
||||
rstate = np.random.default_rng(0)
|
||||
for radius in [0, 2, 4]:
|
||||
# larger defects are less common
|
||||
thresh = 3.25 + 0.25 * radius # larger defects less common
|
||||
tmp_mask = rstate.standard_normal(image_orig.shape[:-1]) > thresh
|
||||
if radius > 0:
|
||||
tmp_mask = binary_dilation(tmp_mask, disk(radius, dtype=bool))
|
||||
mask[tmp_mask] = 1
|
||||
|
||||
# Defect image over the same region in each color channel
|
||||
image_defect = image_orig.copy()
|
||||
for layer in range(image_defect.shape[-1]):
|
||||
image_defect[np.where(mask)] = 0
|
||||
|
||||
if channel_axis is None:
|
||||
image_orig = rgb2gray(image_orig)
|
||||
image_defect = rgb2gray(image_defect)
|
||||
|
||||
image_orig = image_orig.astype(dtype, copy=False)
|
||||
image_defect = image_defect.astype(dtype, copy=False)
|
||||
|
||||
image_defect = np.asarray(image_defect, order=order)
|
||||
image_result = inpaint.inpaint_biharmonic(
|
||||
image_defect, mask, channel_axis=channel_axis,
|
||||
split_into_regions=split_into_regions
|
||||
)
|
||||
assert image_result.dtype == float_dtype
|
||||
|
||||
nrmse_defect = normalized_root_mse(image_orig, image_defect)
|
||||
nrmse_result = normalized_root_mse(img_as_float(image_orig), image_result)
|
||||
assert nrmse_result < 0.2 * nrmse_defect
|
||||
@@ -1,87 +0,0 @@
|
||||
import functools
|
||||
import numpy as np
|
||||
import pytest
|
||||
|
||||
from skimage._shared.testing import assert_
|
||||
from skimage._shared.utils import _supported_float_type
|
||||
from skimage.data import binary_blobs
|
||||
from skimage.data import camera, chelsea
|
||||
from skimage.metrics import mean_squared_error as mse
|
||||
from skimage.restoration import (calibrate_denoiser,
|
||||
denoise_wavelet)
|
||||
from skimage.restoration.j_invariant import _invariant_denoise
|
||||
from skimage.util import img_as_float, random_noise
|
||||
|
||||
test_img = img_as_float(camera())
|
||||
test_img_color = img_as_float(chelsea())
|
||||
test_img_3d = img_as_float(binary_blobs(64, n_dim=3)) / 2
|
||||
noisy_img = random_noise(test_img, mode='gaussian', var=0.01)
|
||||
noisy_img_color = random_noise(test_img_color, mode='gaussian', var=0.01)
|
||||
noisy_img_3d = random_noise(test_img_3d, mode='gaussian', var=0.1)
|
||||
|
||||
_denoise_wavelet = functools.partial(denoise_wavelet, rescale_sigma=True)
|
||||
|
||||
|
||||
def test_invariant_denoise():
|
||||
denoised_img = _invariant_denoise(noisy_img, _denoise_wavelet)
|
||||
|
||||
denoised_mse = mse(denoised_img, test_img)
|
||||
original_mse = mse(noisy_img, test_img)
|
||||
assert_(denoised_mse < original_mse)
|
||||
|
||||
|
||||
@pytest.mark.parametrize('dtype', [np.float16, np.float32, np.float64])
|
||||
def test_invariant_denoise_color(dtype):
|
||||
denoised_img_color = _invariant_denoise(
|
||||
noisy_img_color.astype(dtype), _denoise_wavelet,
|
||||
denoiser_kwargs=dict(channel_axis=-1))
|
||||
denoised_mse = mse(denoised_img_color, test_img_color)
|
||||
original_mse = mse(noisy_img_color, test_img_color)
|
||||
assert denoised_mse < original_mse
|
||||
assert denoised_img_color.dtype == _supported_float_type(dtype)
|
||||
|
||||
|
||||
def test_invariant_denoise_3d():
|
||||
denoised_img_3d = _invariant_denoise(noisy_img_3d, _denoise_wavelet)
|
||||
|
||||
denoised_mse = mse(denoised_img_3d, test_img_3d)
|
||||
original_mse = mse(noisy_img_3d, test_img_3d)
|
||||
assert_(denoised_mse < original_mse)
|
||||
|
||||
|
||||
def test_calibrate_denoiser_extra_output():
|
||||
parameter_ranges = {'sigma': np.linspace(0.1, 1, 5) / 2}
|
||||
_, (parameters_tested, losses) = calibrate_denoiser(
|
||||
noisy_img,
|
||||
_denoise_wavelet,
|
||||
denoise_parameters=parameter_ranges,
|
||||
extra_output=True
|
||||
)
|
||||
|
||||
all_denoised = [_invariant_denoise(noisy_img, _denoise_wavelet,
|
||||
denoiser_kwargs=denoiser_kwargs)
|
||||
for denoiser_kwargs in parameters_tested]
|
||||
|
||||
ground_truth_losses = [mse(img, test_img) for img in all_denoised]
|
||||
assert_(np.argmin(losses) == np.argmin(ground_truth_losses))
|
||||
|
||||
|
||||
def test_calibrate_denoiser():
|
||||
parameter_ranges = {'sigma': np.linspace(0.1, 1, 5) / 2}
|
||||
|
||||
denoiser = calibrate_denoiser(noisy_img, _denoise_wavelet,
|
||||
denoise_parameters=parameter_ranges)
|
||||
|
||||
denoised_mse = mse(denoiser(noisy_img), test_img)
|
||||
original_mse = mse(noisy_img, test_img)
|
||||
assert_(denoised_mse < original_mse)
|
||||
|
||||
|
||||
def test_input_image_not_modified():
|
||||
input_image = noisy_img.copy()
|
||||
|
||||
parameter_ranges = {'sigma': np.random.random(5) / 2}
|
||||
calibrate_denoiser(input_image, _denoise_wavelet,
|
||||
denoise_parameters=parameter_ranges)
|
||||
|
||||
assert_(np.all(noisy_img == input_image))
|
||||
@@ -1,177 +0,0 @@
|
||||
import numpy as np
|
||||
import pytest
|
||||
from scipy import ndimage as ndi
|
||||
from scipy.signal import convolve2d, convolve
|
||||
|
||||
from skimage import restoration, util
|
||||
from skimage._shared import filters
|
||||
from skimage._shared.testing import fetch
|
||||
from skimage._shared.utils import _supported_float_type
|
||||
from skimage.color import rgb2gray
|
||||
from skimage.data import astronaut, camera
|
||||
from skimage.restoration import uft
|
||||
|
||||
test_img = util.img_as_float(camera())
|
||||
|
||||
|
||||
def _get_rtol_atol(dtype):
|
||||
rtol = 1e-3
|
||||
atol = 0
|
||||
if dtype == np.float16:
|
||||
rtol = 1e-2
|
||||
atol = 1e-3
|
||||
elif dtype == np.float32:
|
||||
atol = 1e-5
|
||||
return rtol, atol
|
||||
|
||||
|
||||
@pytest.mark.parametrize('dtype', [np.float16, np.float32, np.float64])
|
||||
@pytest.mark.parametrize('ndim', [1, 2, 3])
|
||||
def test_wiener(dtype, ndim):
|
||||
"""
|
||||
currently only performs pixelwise comparison to
|
||||
precomputed result in 2d case.
|
||||
"""
|
||||
|
||||
psf = np.ones([5] * ndim, dtype=dtype) / 5 ** ndim
|
||||
np.random.seed(0)
|
||||
# for ndim == 2 use camera (to compare to presaved result)
|
||||
if ndim != 2:
|
||||
test_img = np.random.randint(0, 100, [50] * ndim)
|
||||
else:
|
||||
test_img = util.img_as_float(camera())
|
||||
|
||||
data = convolve(test_img, psf, 'same')
|
||||
data += 0.1 * data.std() * np.random.standard_normal(data.shape)
|
||||
data = data.astype(dtype, copy=False)
|
||||
deconvolved = restoration.wiener(data, psf, 0.05)
|
||||
assert deconvolved.dtype == _supported_float_type(dtype)
|
||||
|
||||
if ndim == 2:
|
||||
rtol, atol = _get_rtol_atol(dtype)
|
||||
path = fetch('restoration/tests/camera_wiener.npy')
|
||||
np.testing.assert_allclose(deconvolved, np.load(path),
|
||||
rtol=rtol, atol=atol)
|
||||
|
||||
_, laplacian = uft.laplacian(ndim, data.shape)
|
||||
otf = uft.ir2tf(psf, data.shape, is_real=False)
|
||||
assert otf.real.dtype == _supported_float_type(dtype)
|
||||
deconvolved = restoration.wiener(data, otf, 0.05,
|
||||
reg=laplacian,
|
||||
is_real=False)
|
||||
assert deconvolved.real.dtype == _supported_float_type(dtype)
|
||||
if ndim == 2:
|
||||
np.testing.assert_allclose(np.real(deconvolved),
|
||||
np.load(path),
|
||||
rtol=rtol, atol=atol)
|
||||
|
||||
|
||||
@pytest.mark.parametrize('dtype', [np.float16, np.float32, np.float64])
|
||||
def test_unsupervised_wiener(dtype):
|
||||
psf = np.ones((5, 5), dtype=dtype) / 25
|
||||
data = convolve2d(test_img, psf, 'same')
|
||||
seed = 16829302
|
||||
# keep old-style RandomState here for compatibility with previously stored
|
||||
# reference data in camera_unsup.npy and camera_unsup2.npy
|
||||
rng = np.random.RandomState(seed)
|
||||
data += 0.1 * data.std() * rng.standard_normal(data.shape)
|
||||
data = data.astype(dtype, copy=False)
|
||||
deconvolved, _ = restoration.unsupervised_wiener(data, psf,
|
||||
random_state=seed)
|
||||
float_type = _supported_float_type(dtype)
|
||||
assert deconvolved.dtype == float_type
|
||||
|
||||
rtol, atol = _get_rtol_atol(dtype)
|
||||
path = fetch('restoration/tests/camera_unsup.npy')
|
||||
np.testing.assert_allclose(deconvolved, np.load(path), rtol=rtol,
|
||||
atol=atol)
|
||||
|
||||
_, laplacian = uft.laplacian(2, data.shape)
|
||||
otf = uft.ir2tf(psf, data.shape, is_real=False)
|
||||
assert otf.real.dtype == _supported_float_type(dtype)
|
||||
deconvolved2 = restoration.unsupervised_wiener(
|
||||
data, otf, reg=laplacian, is_real=False,
|
||||
user_params={
|
||||
"callback": lambda x: None,
|
||||
"max_num_iter": 200,
|
||||
"min_num_iter": 30,
|
||||
},
|
||||
random_state=seed)[0]
|
||||
assert deconvolved2.real.dtype == float_type
|
||||
path = fetch('restoration/tests/camera_unsup2.npy')
|
||||
np.testing.assert_allclose(np.real(deconvolved2),
|
||||
np.load(path),
|
||||
rtol=rtol, atol=atol)
|
||||
|
||||
|
||||
def test_unsupervised_wiener_deprecated_user_param():
|
||||
psf = np.ones((5, 5), dtype=float) / 25
|
||||
data = convolve2d(test_img, psf, 'same')
|
||||
otf = uft.ir2tf(psf, data.shape, is_real=False)
|
||||
_, laplacian = uft.laplacian(2, data.shape)
|
||||
restoration.unsupervised_wiener(
|
||||
data, otf, reg=laplacian, is_real=False,
|
||||
user_params={"min_num_iter": 30}, random_state=5
|
||||
)
|
||||
|
||||
|
||||
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), float)
|
||||
point[2, 2] = 1.
|
||||
psf = filters.gaussian(point, sigma=1., mode='reflect')
|
||||
# image shape: (45, 45), as reported in #1172
|
||||
image = util.img_as_float(camera()[65:165, 215:315]) # just the face
|
||||
image_conv = ndi.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)
|
||||
|
||||
|
||||
@pytest.mark.parametrize('ndim', [1, 2, 3])
|
||||
def test_richardson_lucy(ndim):
|
||||
psf = np.ones([5] * ndim, dtype=float) / 5 ** ndim
|
||||
if ndim != 2:
|
||||
test_img = np.random.randint(0, 100, [30] * ndim)
|
||||
else:
|
||||
test_img = util.img_as_float(camera())
|
||||
data = convolve(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, num_iter=5)
|
||||
|
||||
if ndim == 2:
|
||||
path = fetch('restoration/tests/camera_rl.npy')
|
||||
np.testing.assert_allclose(deconvolved, np.load(path), rtol=1e-3)
|
||||
|
||||
|
||||
@pytest.mark.parametrize('dtype_image', [np.float16, np.float32, np.float64])
|
||||
@pytest.mark.parametrize('dtype_psf', [np.float32, np.float64])
|
||||
def test_richardson_lucy_filtered(dtype_image, dtype_psf):
|
||||
if dtype_image == np.float64:
|
||||
atol = 1e-8
|
||||
else:
|
||||
atol = 1e-5
|
||||
test_img_astro = rgb2gray(astronaut())
|
||||
|
||||
psf = np.ones((5, 5), dtype=dtype_psf) / 25
|
||||
data = convolve2d(test_img_astro, psf, 'same')
|
||||
data = data.astype(dtype_image, copy=False)
|
||||
|
||||
deconvolved = restoration.richardson_lucy(data, psf, 5,
|
||||
filter_epsilon=1e-6)
|
||||
assert deconvolved.dtype == _supported_float_type(data.dtype)
|
||||
|
||||
path = fetch('restoration/tests/astronaut_rl.npy')
|
||||
np.testing.assert_allclose(deconvolved, np.load(path), rtol=1e-3,
|
||||
atol=atol)
|
||||
@@ -1,103 +0,0 @@
|
||||
"""
|
||||
Tests for Rolling Ball Filter
|
||||
(skimage.restoration.rolling_ball)
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import pytest
|
||||
|
||||
from skimage import data
|
||||
from skimage.restoration import rolling_ball
|
||||
from skimage.restoration.rolling_ball import ellipsoid_kernel
|
||||
|
||||
|
||||
@pytest.mark.parametrize(
|
||||
'dtype',
|
||||
[np.uint8, np.int32, np.float16, np.float32, np.float64]
|
||||
)
|
||||
def test_ellipsoid_const(dtype):
|
||||
img = 155 * np.ones((100, 100), dtype=dtype)
|
||||
kernel = ellipsoid_kernel((25, 53), 50)
|
||||
background = rolling_ball(img, kernel=kernel)
|
||||
assert np.allclose(img - background, np.zeros_like(img))
|
||||
assert background.dtype == img.dtype
|
||||
|
||||
|
||||
def test_nan_const():
|
||||
img = 123 * np.ones((100, 100), dtype=float)
|
||||
img[20, 20] = np.nan
|
||||
img[50, 53] = np.nan
|
||||
|
||||
kernel_shape = (10, 10)
|
||||
x = np.arange(-kernel_shape[1] // 2,
|
||||
kernel_shape[1] // 2 + 1)[np.newaxis, :]
|
||||
y = np.arange(-kernel_shape[0] // 2,
|
||||
kernel_shape[0] // 2 + 1)[:, np.newaxis]
|
||||
expected_img = np.zeros_like(img)
|
||||
expected_img[y + 20, x + 20] = np.nan
|
||||
expected_img[y + 50, x + 53] = np.nan
|
||||
kernel = ellipsoid_kernel(kernel_shape, 100)
|
||||
background = rolling_ball(
|
||||
img,
|
||||
kernel=kernel,
|
||||
nansafe=True
|
||||
)
|
||||
assert np.allclose(img - background, expected_img, equal_nan=True)
|
||||
|
||||
|
||||
@pytest.mark.parametrize("radius", [1, 2.5, 10.346, 50])
|
||||
def test_const_image(radius):
|
||||
# infinite plane light source at top left corner
|
||||
img = 23 * np.ones((100, 100), dtype=np.uint8)
|
||||
background = rolling_ball(img, radius=radius)
|
||||
assert np.allclose(img - background, np.zeros_like(img))
|
||||
|
||||
|
||||
def test_radial_gradient():
|
||||
# spot light source at top left corner
|
||||
spot_radius = 50
|
||||
x, y = np.meshgrid(range(5), range(5))
|
||||
img = np.sqrt(np.clip(spot_radius ** 2 - y ** 2 - x ** 2, 0, None))
|
||||
|
||||
background = rolling_ball(img, radius=5)
|
||||
assert np.allclose(img - background, np.zeros_like(img))
|
||||
|
||||
|
||||
def test_linear_gradient():
|
||||
# linear light source centered at top left corner
|
||||
x, y = np.meshgrid(range(100), range(100))
|
||||
img = (y * 20 + x * 20)
|
||||
|
||||
expected_img = 19 * np.ones_like(img)
|
||||
expected_img[0, 0] = 0
|
||||
|
||||
background = rolling_ball(img, radius=1)
|
||||
assert np.allclose(img - background, expected_img)
|
||||
|
||||
|
||||
@pytest.mark.parametrize("radius", [2, 10, 12.5, 50])
|
||||
def test_preserve_peaks(radius):
|
||||
x, y = np.meshgrid(range(100), range(100))
|
||||
img = 0 * x + 0 * y + 10
|
||||
img[10, 10] = 20
|
||||
img[20, 20] = 35
|
||||
img[45, 26] = 156
|
||||
|
||||
expected_img = img - 10
|
||||
background = rolling_ball(img, radius=radius)
|
||||
assert np.allclose(img - background, expected_img)
|
||||
|
||||
|
||||
@pytest.mark.parametrize("num_threads", [None, 1, 2])
|
||||
def test_threads(num_threads):
|
||||
# not testing if we use multiple threads
|
||||
# just checking if the API throws an exception
|
||||
img = 23 * np.ones((100, 100), dtype=np.uint8)
|
||||
rolling_ball(img, radius=10, num_threads=num_threads)
|
||||
rolling_ball(img, radius=10, nansafe=True, num_threads=num_threads)
|
||||
|
||||
|
||||
def test_ndim():
|
||||
image = data.cells3d()[:5, 1, ...]
|
||||
kernel = ellipsoid_kernel((3, 100, 100), 100)
|
||||
rolling_ball(image, kernel=kernel)
|
||||
@@ -1,217 +0,0 @@
|
||||
import numpy as np
|
||||
from skimage.restoration import unwrap_phase
|
||||
import sys
|
||||
|
||||
from skimage._shared import testing
|
||||
from skimage._shared.testing import (assert_array_almost_equal_nulp,
|
||||
assert_almost_equal, assert_array_equal,
|
||||
assert_, skipif)
|
||||
from skimage._shared._warnings import expected_warnings
|
||||
|
||||
|
||||
def assert_phase_almost_equal(a, b, *args, **kwargs):
|
||||
"""An assert_almost_equal insensitive to phase shifts of n*2*pi."""
|
||||
shift = 2 * np.pi * np.round((b.mean() - a.mean()) / (2 * np.pi))
|
||||
with expected_warnings([r'invalid value encountered|\A\Z',
|
||||
r'divide by zero encountered|\A\Z']):
|
||||
print('assert_phase_allclose, abs', np.max(np.abs(a - (b - shift))))
|
||||
print('assert_phase_allclose, rel',
|
||||
np.max(np.abs((a - (b - shift)) / a)))
|
||||
if np.ma.isMaskedArray(a):
|
||||
assert_(np.ma.isMaskedArray(b))
|
||||
assert_array_equal(a.mask, b.mask)
|
||||
assert_(a.fill_value == b.fill_value)
|
||||
au = np.asarray(a)
|
||||
bu = np.asarray(b)
|
||||
with expected_warnings([r'invalid value encountered|\A\Z',
|
||||
r'divide by zero encountered|\A\Z']):
|
||||
print('assert_phase_allclose, no mask, abs',
|
||||
np.max(np.abs(au - (bu - shift))))
|
||||
print('assert_phase_allclose, no mask, rel',
|
||||
np.max(np.abs((au - (bu - shift)) / au)))
|
||||
assert_array_almost_equal_nulp(a + shift, b, *args, **kwargs)
|
||||
|
||||
|
||||
def check_unwrap(image, mask=None):
|
||||
image_wrapped = np.angle(np.exp(1j * image))
|
||||
if mask is not None:
|
||||
print('Testing a masked image')
|
||||
image = np.ma.array(image, mask=mask, fill_value=0.5)
|
||||
image_wrapped = np.ma.array(image_wrapped, mask=mask, fill_value=0.5)
|
||||
image_unwrapped = unwrap_phase(image_wrapped, seed=0)
|
||||
assert_phase_almost_equal(image_unwrapped, image)
|
||||
|
||||
|
||||
def test_unwrap_1d():
|
||||
image = np.linspace(0, 10 * np.pi, 100)
|
||||
check_unwrap(image)
|
||||
# Masked arrays are not allowed in 1D
|
||||
with testing.raises(ValueError):
|
||||
check_unwrap(image, True)
|
||||
# wrap_around is not allowed in 1D
|
||||
with testing.raises(ValueError):
|
||||
unwrap_phase(image, True, seed=0)
|
||||
|
||||
|
||||
@testing.parametrize("check_with_mask", (False, True))
|
||||
def test_unwrap_2d(check_with_mask):
|
||||
mask = None
|
||||
x, y = np.ogrid[:8, :16]
|
||||
image = 2 * np.pi * (x * 0.2 + y * 0.1)
|
||||
if check_with_mask:
|
||||
mask = np.zeros(image.shape, dtype=bool)
|
||||
mask[4:6, 4:8] = True
|
||||
check_unwrap(image, mask)
|
||||
|
||||
|
||||
@testing.parametrize("check_with_mask", (False, True))
|
||||
def test_unwrap_3d(check_with_mask):
|
||||
mask = None
|
||||
x, y, z = np.ogrid[:8, :12, :16]
|
||||
image = 2 * np.pi * (x * 0.2 + y * 0.1 + z * 0.05)
|
||||
if check_with_mask:
|
||||
mask = np.zeros(image.shape, dtype=bool)
|
||||
mask[4:6, 4:6, 1:3] = True
|
||||
check_unwrap(image, mask)
|
||||
|
||||
|
||||
def check_wrap_around(ndim, axis):
|
||||
# create a ramp, but with the last pixel along axis equalling the first
|
||||
elements = 100
|
||||
ramp = np.linspace(0, 12 * np.pi, elements)
|
||||
ramp[-1] = ramp[0]
|
||||
image = ramp.reshape(tuple([elements if n == axis else 1
|
||||
for n in range(ndim)]))
|
||||
image_wrapped = np.angle(np.exp(1j * image))
|
||||
|
||||
index_first = tuple([0] * ndim)
|
||||
index_last = tuple([-1 if n == axis else 0 for n in range(ndim)])
|
||||
# unwrap the image without wrap around
|
||||
# We do not want warnings about length 1 dimensions
|
||||
with expected_warnings([r'Image has a length 1 dimension|\A\Z']):
|
||||
image_unwrap_no_wrap_around = unwrap_phase(image_wrapped, seed=0)
|
||||
print('endpoints without wrap_around:',
|
||||
image_unwrap_no_wrap_around[index_first],
|
||||
image_unwrap_no_wrap_around[index_last])
|
||||
# without wrap around, the endpoints of the image should differ
|
||||
assert_(abs(image_unwrap_no_wrap_around[index_first] -
|
||||
image_unwrap_no_wrap_around[index_last]) > np.pi)
|
||||
# unwrap the image with wrap around
|
||||
wrap_around = [n == axis for n in range(ndim)]
|
||||
# We do not want warnings about length 1 dimensions
|
||||
with expected_warnings([r'Image has a length 1 dimension.|\A\Z']):
|
||||
image_unwrap_wrap_around = unwrap_phase(image_wrapped, wrap_around,
|
||||
seed=0)
|
||||
print('endpoints with wrap_around:',
|
||||
image_unwrap_wrap_around[index_first],
|
||||
image_unwrap_wrap_around[index_last])
|
||||
# with wrap around, the endpoints of the image should be equal
|
||||
assert_almost_equal(image_unwrap_wrap_around[index_first],
|
||||
image_unwrap_wrap_around[index_last])
|
||||
|
||||
|
||||
dim_axis = [(ndim, axis) for ndim in (2, 3) for axis in range(ndim)]
|
||||
|
||||
|
||||
@skipif(sys.version_info[:2] == (3, 4),
|
||||
reason="Doesn't work with python 3.4. See issue #3079")
|
||||
@testing.parametrize("ndim, axis", dim_axis)
|
||||
def test_wrap_around(ndim, axis):
|
||||
check_wrap_around(ndim, axis)
|
||||
|
||||
|
||||
def test_mask():
|
||||
length = 100
|
||||
ramps = [np.linspace(0, 4 * np.pi, length),
|
||||
np.linspace(0, 8 * np.pi, length),
|
||||
np.linspace(0, 6 * np.pi, length)]
|
||||
image = np.vstack(ramps)
|
||||
mask_1d = np.ones((length,), dtype=bool)
|
||||
mask_1d[0] = mask_1d[-1] = False
|
||||
for i in range(len(ramps)):
|
||||
# mask all ramps but the i'th one
|
||||
mask = np.zeros(image.shape, dtype=bool)
|
||||
mask |= mask_1d.reshape(1, -1)
|
||||
mask[i, :] = False # unmask i'th ramp
|
||||
image_wrapped = np.ma.array(np.angle(np.exp(1j * image)), mask=mask)
|
||||
image_unwrapped = unwrap_phase(image_wrapped)
|
||||
image_unwrapped -= image_unwrapped[0, 0] # remove phase shift
|
||||
# The end of the unwrapped array should have value equal to the
|
||||
# endpoint of the unmasked ramp
|
||||
assert_array_almost_equal_nulp(image_unwrapped[:, -1], image[i, -1])
|
||||
assert_(np.ma.isMaskedArray(image_unwrapped))
|
||||
|
||||
# Same tests, but forcing use of the 3D unwrapper by reshaping
|
||||
with expected_warnings(['length 1 dimension']):
|
||||
shape = (1,) + image_wrapped.shape
|
||||
image_wrapped_3d = image_wrapped.reshape(shape)
|
||||
image_unwrapped_3d = unwrap_phase(image_wrapped_3d)
|
||||
# remove phase shift
|
||||
image_unwrapped_3d -= image_unwrapped_3d[0, 0, 0]
|
||||
assert_array_almost_equal_nulp(image_unwrapped_3d[:, :, -1],
|
||||
image[i, -1])
|
||||
|
||||
|
||||
def test_invalid_input():
|
||||
with testing.raises(ValueError):
|
||||
unwrap_phase(np.zeros([]))
|
||||
with testing.raises(ValueError):
|
||||
unwrap_phase(np.zeros((1, 1, 1, 1)))
|
||||
with testing.raises(ValueError):
|
||||
unwrap_phase(np.zeros((1, 1)), 3 * [False])
|
||||
with testing.raises(ValueError):
|
||||
unwrap_phase(np.zeros((1, 1)), 'False')
|
||||
|
||||
|
||||
def test_unwrap_3d_middle_wrap_around():
|
||||
# Segmentation fault in 3D unwrap phase with middle dimension connected
|
||||
# GitHub issue #1171
|
||||
image = np.zeros((20, 30, 40), dtype=np.float32)
|
||||
unwrap = unwrap_phase(image, wrap_around=[False, True, False])
|
||||
assert_(np.all(unwrap == 0))
|
||||
|
||||
|
||||
def test_unwrap_2d_compressed_mask():
|
||||
# ValueError when image is masked array with a compressed mask (no masked
|
||||
# elements). GitHub issue #1346
|
||||
image = np.ma.zeros((10, 10))
|
||||
unwrap = unwrap_phase(image)
|
||||
assert_(np.all(unwrap == 0))
|
||||
|
||||
|
||||
def test_unwrap_2d_all_masked():
|
||||
# Segmentation fault when image is masked array with a all elements masked
|
||||
# GitHub issue #1347
|
||||
# all elements masked
|
||||
image = np.ma.zeros((10, 10))
|
||||
image[:] = np.ma.masked
|
||||
unwrap = unwrap_phase(image)
|
||||
assert_(np.ma.isMaskedArray(unwrap))
|
||||
assert_(np.all(unwrap.mask))
|
||||
|
||||
# 1 unmasked element, still zero edges
|
||||
image = np.ma.zeros((10, 10))
|
||||
image[:] = np.ma.masked
|
||||
image[0, 0] = 0
|
||||
unwrap = unwrap_phase(image)
|
||||
assert_(np.ma.isMaskedArray(unwrap))
|
||||
assert_(np.sum(unwrap.mask) == 99) # all but one masked
|
||||
assert_(unwrap[0, 0] == 0)
|
||||
|
||||
|
||||
def test_unwrap_3d_all_masked():
|
||||
# all elements masked
|
||||
image = np.ma.zeros((10, 10, 10))
|
||||
image[:] = np.ma.masked
|
||||
unwrap = unwrap_phase(image)
|
||||
assert_(np.ma.isMaskedArray(unwrap))
|
||||
assert_(np.all(unwrap.mask))
|
||||
|
||||
# 1 unmasked element, still zero edges
|
||||
image = np.ma.zeros((10, 10, 10))
|
||||
image[:] = np.ma.masked
|
||||
image[0, 0, 0] = 0
|
||||
unwrap = unwrap_phase(image)
|
||||
assert_(np.ma.isMaskedArray(unwrap))
|
||||
assert_(np.sum(unwrap.mask) == 999) # all but one masked
|
||||
assert_(unwrap[0, 0, 0] == 0)
|
||||
@@ -1,454 +0,0 @@
|
||||
r"""Function of unitary fourier transform (uft) and utilities
|
||||
|
||||
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
|
||||
|
||||
so the Fourier transform has the same energy as the original image
|
||||
(see ``image_quad_norm`` function). The transform is applied from the
|
||||
last axis for performance (assuming a C-order array input).
|
||||
|
||||
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
|
||||
|
||||
"""
|
||||
|
||||
|
||||
import numpy as np
|
||||
import scipy.fft as fft
|
||||
|
||||
from .._shared.utils import _supported_float_type
|
||||
|
||||
|
||||
|
||||
def ufftn(inarray, dim=None):
|
||||
"""N-dimensional unitary Fourier transform.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
inarray : ndarray
|
||||
The array to transform.
|
||||
dim : int, optional
|
||||
The last axis along which 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 = fft.fftn(inarray, axes=range(-dim, 0), norm='ortho')
|
||||
return outarray
|
||||
|
||||
|
||||
def uifftn(inarray, dim=None):
|
||||
"""N-dimensional unitary inverse Fourier transform.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
inarray : ndarray
|
||||
The array to transform.
|
||||
dim : int, optional
|
||||
The last axis along which 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 = fft.ifftn(inarray, axes=range(-dim, 0), norm='ortho')
|
||||
return outarray
|
||||
|
||||
|
||||
def urfftn(inarray, dim=None):
|
||||
"""N-dimensional real unitary Fourier transform.
|
||||
|
||||
This transform considers the Hermitian property of the transform on
|
||||
real-valued input.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
inarray : ndarray, shape (M, N, ..., P)
|
||||
The array to transform.
|
||||
dim : int, optional
|
||||
The last axis along which to compute the transform. All
|
||||
axes by default.
|
||||
|
||||
Returns
|
||||
-------
|
||||
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 has a Hermitian property and
|
||||
redundant values are not computed or 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 = fft.rfftn(inarray, axes=range(-dim, 0), norm='ortho')
|
||||
return outarray
|
||||
|
||||
|
||||
def uirfftn(inarray, dim=None, shape=None):
|
||||
"""N-dimensional inverse real unitary Fourier transform.
|
||||
|
||||
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 last axis along which to compute the transform. All
|
||||
axes by default.
|
||||
shape : tuple of int, optional
|
||||
The shape of the output. The shape of ``rfft`` is ambiguous in
|
||||
case of odd-valued input shape. In this case, this parameter
|
||||
should be provided. See ``np.fft.irfftn``.
|
||||
|
||||
Returns
|
||||
-------
|
||||
outarray : ndarray
|
||||
The unitary N-D inverse real Fourier transform of ``inarray``.
|
||||
|
||||
Notes
|
||||
-----
|
||||
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
|
||||
--------
|
||||
>>> 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 = fft.irfftn(inarray, shape, axes=range(-dim, 0), norm='ortho')
|
||||
return outarray
|
||||
|
||||
|
||||
def ufft2(inarray):
|
||||
"""2-dimensional unitary Fourier transform.
|
||||
|
||||
Compute the Fourier transform on the last 2 axes.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
inarray : ndarray
|
||||
The array to transform.
|
||||
|
||||
Returns
|
||||
-------
|
||||
outarray : ndarray (same shape as 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-dimensional 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 as 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-dimensional real unitary Fourier transform
|
||||
|
||||
Compute the real Fourier transform on the last 2 axes. This
|
||||
transform considers the Hermitian property of the transform from
|
||||
complex to real-valued input.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
inarray : ndarray, shape (M, N, ..., P)
|
||||
The array to transform.
|
||||
|
||||
Returns
|
||||
-------
|
||||
outarray : ndarray, shape (M, N, ..., 2 * (P - 1))
|
||||
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-dimensional inverse real unitary Fourier transform.
|
||||
|
||||
Compute the real inverse Fourier transform on the last 2 axes.
|
||||
This transform considers the Hermitian property of the transform
|
||||
from complex to real-valued input.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
inarray : ndarray, shape (M, N, ..., P)
|
||||
The array to transform.
|
||||
shape : tuple of int, optional
|
||||
The shape of the output. The shape of ``rfft`` is ambiguous in
|
||||
case of odd-valued input shape. In this case, this parameter
|
||||
should be provided. See ``np.fft.irfftn``.
|
||||
|
||||
Returns
|
||||
-------
|
||||
outarray : ndarray, shape (M, N, ..., 2 * (P - 1))
|
||||
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 the quadratic norm of images in Fourier space.
|
||||
|
||||
This function detects whether the input image satisfies the
|
||||
Hermitian property.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
inarray : ndarray
|
||||
Input image. The image data should reside in the final 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 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))
|
||||
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 an impulse response (IR).
|
||||
|
||||
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 impulse responses.
|
||||
shape : tuple of int
|
||||
A tuple of integer corresponding to the target shape of the
|
||||
transfer function.
|
||||
dim : int, optional
|
||||
The last axis along which to compute the transform. All
|
||||
axes by default.
|
||||
is_real : boolean, optional
|
||||
If True (default), imp_resp is supposed real and the Hermitian property
|
||||
is used with rfftn Fourier transform.
|
||||
|
||||
Returns
|
||||
-------
|
||||
y : complex ndarray
|
||||
The transfer 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-dimensional IR with
|
||||
an arbitrary number of IR. The individual IR must be accessed
|
||||
through the first axes. The last ``dim`` axes contain the space
|
||||
definition.
|
||||
"""
|
||||
if not dim:
|
||||
dim = imp_resp.ndim
|
||||
# Zero padding and fill
|
||||
irpadded_dtype = _supported_float_type(imp_resp)
|
||||
irpadded = np.zeros(shape, dtype=irpadded_dtype)
|
||||
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)
|
||||
|
||||
func = fft.rfftn if is_real else fft.fftn
|
||||
out = func(irpadded, axes=(range(-dim, 0)))
|
||||
|
||||
# TODO: remove .astype call once SciPy >= 1.4 is required
|
||||
cplx_dtype = np.promote_types(irpadded_dtype, np.complex64)
|
||||
return out.astype(cplx_dtype, copy=False)
|
||||
|
||||
|
||||
def laplacian(ndim, shape, is_real=True):
|
||||
"""Return the transfer function of the Laplacian.
|
||||
|
||||
Laplacian is the second order difference, on row and column.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
ndim : int
|
||||
The dimension of the Laplacian.
|
||||
shape : tuple
|
||||
The support on which to compute the transfer function.
|
||||
is_real : boolean, optional
|
||||
If True (default), 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.
|
||||
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
|
||||
@@ -1,113 +0,0 @@
|
||||
import numpy as np
|
||||
|
||||
from .._shared.utils import warn
|
||||
|
||||
from ._unwrap_1d import unwrap_1d
|
||||
from ._unwrap_2d import unwrap_2d
|
||||
from ._unwrap_3d import unwrap_3d
|
||||
|
||||
|
||||
def unwrap_phase(image, wrap_around=False, seed=None):
|
||||
'''Recover the original from a wrapped phase image.
|
||||
|
||||
From an image wrapped to lie in the interval [-pi, pi), recover the
|
||||
original, unwrapped image.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : 1D, 2D or 3D ndarray of floats, optionally a masked array
|
||||
The values should be in the range [-pi, pi). If a masked array is
|
||||
provided, the masked entries will not be changed, and their values
|
||||
will not be used to guide the unwrapping of neighboring, unmasked
|
||||
values. Masked 1D arrays are not allowed, and will raise a
|
||||
`ValueError`.
|
||||
wrap_around : bool or sequence of bool, optional
|
||||
When an element of the sequence is `True`, the unwrapping process
|
||||
will regard the edges along the corresponding axis of the image to be
|
||||
connected and use this connectivity to guide the phase unwrapping
|
||||
process. If only a single boolean is given, it will apply to all axes.
|
||||
Wrap around is not supported for 1D arrays.
|
||||
seed : int, optional
|
||||
Unwrapping 2D or 3D images uses random initialization. This sets the
|
||||
seed of the PRNG to achieve deterministic behavior.
|
||||
|
||||
Returns
|
||||
-------
|
||||
image_unwrapped : array_like, double
|
||||
Unwrapped image of the same shape as the input. If the input `image`
|
||||
was a masked array, the mask will be preserved.
|
||||
|
||||
Raises
|
||||
------
|
||||
ValueError
|
||||
If called with a masked 1D array or called with a 1D array and
|
||||
``wrap_around=True``.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> c0, c1 = np.ogrid[-1:1:128j, -1:1:128j]
|
||||
>>> image = 12 * np.pi * np.exp(-(c0**2 + c1**2))
|
||||
>>> image_wrapped = np.angle(np.exp(1j * image))
|
||||
>>> image_unwrapped = unwrap_phase(image_wrapped)
|
||||
>>> np.std(image_unwrapped - image) < 1e-6 # A constant offset is normal
|
||||
True
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] Miguel Arevallilo Herraez, David R. Burton, Michael J. Lalor,
|
||||
and Munther A. Gdeisat, "Fast two-dimensional phase-unwrapping
|
||||
algorithm based on sorting by reliability following a noncontinuous
|
||||
path", Journal Applied Optics, Vol. 41, No. 35 (2002) 7437,
|
||||
.. [2] Abdul-Rahman, H., Gdeisat, M., Burton, D., & Lalor, M., "Fast
|
||||
three-dimensional phase-unwrapping algorithm based on sorting by
|
||||
reliability following a non-continuous path. In W. Osten,
|
||||
C. Gorecki, & E. L. Novak (Eds.), Optical Metrology (2005) 32--40,
|
||||
International Society for Optics and Photonics.
|
||||
'''
|
||||
if image.ndim not in (1, 2, 3):
|
||||
raise ValueError('Image must be 1, 2, or 3 dimensional')
|
||||
if isinstance(wrap_around, bool):
|
||||
wrap_around = [wrap_around] * image.ndim
|
||||
elif (hasattr(wrap_around, '__getitem__')
|
||||
and not isinstance(wrap_around, str)):
|
||||
if len(wrap_around) != image.ndim:
|
||||
raise ValueError('Length of `wrap_around` must equal the '
|
||||
'dimensionality of image')
|
||||
wrap_around = [bool(wa) for wa in wrap_around]
|
||||
else:
|
||||
raise ValueError('`wrap_around` must be a bool or a sequence with '
|
||||
'length equal to the dimensionality of image')
|
||||
if image.ndim == 1:
|
||||
if np.ma.isMaskedArray(image):
|
||||
raise ValueError('1D masked images cannot be unwrapped')
|
||||
if wrap_around[0]:
|
||||
raise ValueError('`wrap_around` is not supported for 1D images')
|
||||
if image.ndim in (2, 3) and 1 in image.shape:
|
||||
warn('Image has a length 1 dimension. Consider using an '
|
||||
'array of lower dimensionality to use a more efficient '
|
||||
'algorithm')
|
||||
|
||||
if np.ma.isMaskedArray(image):
|
||||
mask = np.require(np.ma.getmaskarray(image), np.uint8, ['C'])
|
||||
else:
|
||||
mask = np.zeros_like(image, dtype=np.uint8, order='C')
|
||||
|
||||
image_not_masked = np.asarray(
|
||||
np.ma.getdata(image), dtype=np.float64, order='C')
|
||||
image_unwrapped = np.empty_like(image, dtype=np.float64, order='C',
|
||||
subok=False)
|
||||
|
||||
if image.ndim == 1:
|
||||
unwrap_1d(image_not_masked, image_unwrapped)
|
||||
elif image.ndim == 2:
|
||||
unwrap_2d(image_not_masked, mask, image_unwrapped,
|
||||
wrap_around, seed)
|
||||
elif image.ndim == 3:
|
||||
unwrap_3d(image_not_masked, mask, image_unwrapped,
|
||||
wrap_around, seed)
|
||||
|
||||
if np.ma.isMaskedArray(image):
|
||||
return np.ma.array(image_unwrapped, mask=mask,
|
||||
fill_value=image.fill_value)
|
||||
else:
|
||||
return image_unwrapped
|
||||
Reference in New Issue
Block a user