Some improvements to new deconvolution code

This commit is contained in:
Christoph Deil
2013-12-10 22:45:09 +01:00
committed by François Orieux
parent 17cc22802c
commit e25ccf61f6
4 changed files with 37 additions and 48 deletions
+2 -1
View File
@@ -1,6 +1,7 @@
# -*- coding: utf-8 -*-
"""Skimage module for image deconvolution
This module implement various algorithm of the literarture for image
This module implement various algorithm of the literature for image
deconvolution.
References
@@ -1,38 +1,33 @@
import warnings
from os.path import abspath, dirname, join as pjoin
import numpy as np
import numpy.testing.assert_array_almost_equal
from scipy.signal import convolve2d
from skimage.data import camera
from skimage import deconvolution
from scipy.signal import convolve2d as conv2
from skimage import data, deconvolution
# Test deconvolution
# ===========================
test_img = data.camera().astype(np.float)
test_img = camera().astype(np.float)
def test_wiener():
psf = np.ones((5, 5))
data = conv2(test_img, psf, 'same')
data = convolve2d(test_img, psf, 'same')
np.random.seed(0)
data += 0.1 * data.std() * np.random.standard_normal(data.shape)
deconvolued = deconvolution.wiener(data, psf, 25)
numpy.testing.assert_array_almost_equal(deconvolued,
np.load("./camera_wiener.npy"))
path = pjoin(dirname(abspath(__file__)), 'camera_wiener.npy')
np.testing.assert_array_almost_equal(deconvolued, np.load(path))
def test_unsupervised_wiener():
psf = np.ones((5, 5))
data = conv2(test_img, psf, 'same')
data = convolve2d(test_img, psf, 'same')
np.random.seed(0)
data += 0.1 * data.std() * np.random.standard_normal(data.shape)
deconvolued, _ = deconvolution.unsupervised_wiener(data, psf)
numpy.testing.assert_array_almost_equal(deconvolued,
np.load("./camera_unsup.npy"))
path = pjoin(dirname(abspath(__file__)), 'camera_unsup.npy')
np.testing.assert_array_almost_equal(deconvolued, np.load(path))
def test_rychardson_lucy():
def test_richardson_lucy():
return True
+6 -6
View File
@@ -1,7 +1,7 @@
# -*- coding: utf-8 -*-
# uft.py --- Unitary fourier transform
# Copyright (c) 2011, 2012, 2013 François Orieux <orieux@iap.fr>
# Copyright (c) 2011, 2012, 2013 Franois Orieux <orieux@iap.fr>
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation files
@@ -63,12 +63,12 @@ except ImportError:
" by using fftw library.")
ANFFTMOD = False
__author__ = "François Orieux"
__author__ = "Franois Orieux"
__copyright__ = "Copyright (C) 2011, 2012, 2013 F. Orieux <orieux@iap.fr>"
__credits__ = ["François Orieux"]
__credits__ = ["Franois Orieux"]
__license__ = "mit"
__version__ = "0.1.0"
__maintainer__ = "François Orieux"
__maintainer__ = "Franois Orieux"
__email__ = "orieux@iap.fr"
__status__ = "development"
__url__ = ""
@@ -427,12 +427,12 @@ def ir2tf(imp_resp, shape, dim=None, real=True):
for i, s in enumerate(imp_resp.shape)])
if real:
if anfft:
if ANFFTMOD:
return anfft.rfftn(irpadded, k=dim)
else:
return np.fft.rfftn(irpadded, axes=range(-dim, 0))
else:
if anfft:
if ANFFTMOD:
return anfft.fftn(irpadded, k=dim)
else:
return np.fft.fftn(irpadded, axes=range(-dim, 0))
+17 -24
View File
@@ -22,12 +22,8 @@
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
# Commentary:
"""Implementations deconvolution functions"""
# code:
from __future__ import division
import numpy as np
@@ -81,15 +77,16 @@ def wiener(data, psf, reg_val, reg=None, real=True):
im_deconv : (M, N) ndarray
The deconvolued data
Exemples
Examples
--------
>>> import numpy as np
>>> from skimage import color, data, deconvolution
>>> lena = color.rgb2gray(data.lena())
>>> from scipy.signal import convolve2d as conv2
>>> from scipy.signal import convolve2d
>>> psf = np.ones((5, 5))
>>> lena = conv2(lena, psf, 'same')
>>> lena = convolve2d(lena, psf, 'same')
>>> lena += 0.1 * lena.std() * np.random.standard_normal(lena.shape)
>>> deconvolued_lena = deconvolution.wiener(lena, psf, 1100)
>>> deconvolved_lena = deconvolution.wiener(lena, psf, 1100)
References
----------
@@ -184,15 +181,16 @@ def unsupervised_wiener(data, psf, reg=None, user_params=None):
exists, the current image sample. This function can be used to
store the sample, or compute other moments than the mean.
Exemples
Examples
--------
>>> import numpy as np
>>> from skimage import color, data, deconvolution
>>> lena = color.rgb2gray(data.lena())
>>> from scipy.signal import convolve2d as conv2
>>> from scipy.signal import convolve2d
>>> psf = np.ones((5, 5))
>>> lena = conv2(lena, psf, 'same')
>>> lena = convolve2d(lena, psf, 'same')
>>> lena += 0.1 * lena.std() * np.random.standard_normal(lena.shape)
>>> deconvolued_lena = deconvolution.unsupervised_wiener(lena, psf)
>>> deconvolved_lena = deconvolution.unsupervised_wiener(lena, psf)
References
----------
@@ -201,7 +199,7 @@ def unsupervised_wiener(data, psf, reg=None, user_params=None):
spread function parameters for Wiener-Hunt deconvolution",
J. Opt. Soc. Am. A 27, 1593-1607 (2010)
http://www.opticsinfobase.org/josaa/abstract.cfm?URI=josaa-27-7-1593
http://www.opticsinfobase.org/josaa/abstract.cfm?URI=josaa-27-7-1593
"""
params = {'threshold': 1e-4, 'max_iter': 200,
'min_iter': 30, 'burnin': 15, 'callback': None}
@@ -285,11 +283,8 @@ def unsupervised_wiener(data, psf, reg=None, user_params=None):
def richardson_lucy(data, psf, iterations=50):
"""Richardson lucy deconvolution
"""Richardson-Lucy deconvolution.
This a python conversion of the wikipedia page on this algorithm. See
http://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution
Parameters
----------
@@ -297,28 +292,26 @@ def richardson_lucy(data, psf, iterations=50):
The data
psf : ndarray
The impulsionnal response in real space.
The point spread function
iterations : int
dictionary of gibbs parameters. See below.
Number of iterations
Returns
-------
im_deconv : ndarray
The deconvolued image
The deconvolved image
References
----------
.. [2] Richardson, William Hadley, "Bayesian-Based Iterative
Method of Image Restoration". JOSA 62 (1):
5559. doi:10.1364/JOSA.62.000055, 1972
.. [2] http://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution
"""
data = data.astype(np.float)
psf = psf.astype(np.float)
im_deconv = 0.5 * np.ones(data.shape)
psf_mirror = psf[::-1, ::-1]
for iteration in range(iterations):
for _ in range(iterations):
relative_blur = data / conv2(im_deconv, psf, 'same')
im_deconv *= conv2(relative_blur, psf_mirror, 'same')