/
deconvolution.py
421 lines (353 loc) · 15.5 KB
/
deconvolution.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
"""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 \left( |\Lambda_H|^2 + \lambda |\Lambda_D|^2 \right)^{-1}
\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, *, rng=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.
rng : {`numpy.random.Generator`, int}, optional
Pseudo-random number generator.
By default, a PCG64 generator is used (see :func:`numpy.random.default_rng`).
If `rng` is an int, it is used to seed the generator.
.. 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(rng)
# 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