From 5aabc199aa1ad32147bf7a37122725bce7a821c0 Mon Sep 17 00:00:00 2001 From: scottsievert Date: Mon, 30 Nov 2015 23:25:44 -0600 Subject: [PATCH 1/6] uses fftconvolve instead of convolve2d --- skimage/restoration/deconvolution.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/skimage/restoration/deconvolution.py b/skimage/restoration/deconvolution.py index 632d678846f..7cb2e3d703b 100644 --- a/skimage/restoration/deconvolution.py +++ b/skimage/restoration/deconvolution.py @@ -7,7 +7,7 @@ import numpy as np import numpy.random as npr -from scipy.signal import convolve2d +from scipy.signal import fftconvolve from . import uft @@ -370,8 +370,8 @@ def richardson_lucy(image, psf, iterations=50, clip=True): im_deconv = 0.5 * np.ones(image.shape) psf_mirror = psf[::-1, ::-1] for _ in range(iterations): - relative_blur = image / convolve2d(im_deconv, psf, 'same') - im_deconv *= convolve2d(relative_blur, psf_mirror, 'same') + relative_blur = image / fftconvolve(im_deconv, psf, 'same') + im_deconv *= fftconvolve(relative_blur, psf_mirror, 'same') if clip: im_deconv[im_deconv > 1] = 1 From 1be8336599d3f451a65e97a6c4a8aea7bcb8558c Mon Sep 17 00:00:00 2001 From: scottsievert Date: Thu, 10 Dec 2015 19:47:02 -0600 Subject: [PATCH 2/6] adds decision when to use fftconvolve/convolve2d --- skimage/restoration/deconvolution.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/skimage/restoration/deconvolution.py b/skimage/restoration/deconvolution.py index 7cb2e3d703b..6d922aa242e 100644 --- a/skimage/restoration/deconvolution.py +++ b/skimage/restoration/deconvolution.py @@ -365,13 +365,22 @@ def richardson_lucy(image, psf, iterations=50, clip=True): ---------- .. [1] http://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution """ + direct_time = lambda n, m, k, l: k*l * n*m + def fft_time(m, n, k, l): + return m*np.log(m) + n*np.log(n) + k*np.log(k) + l*np.log(l) + + time_ratio = 71.468 * fft_time(*image.shape, *psf.shape) + time_ratio /= direct_time(*image.shape, *psf.shape) + convolve_method = fftconvolve if time_ratio <= 1 else convolve2d + image = image.astype(np.float) psf = psf.astype(np.float) im_deconv = 0.5 * np.ones(image.shape) psf_mirror = psf[::-1, ::-1] + for _ in range(iterations): - relative_blur = image / fftconvolve(im_deconv, psf, 'same') - im_deconv *= fftconvolve(relative_blur, psf_mirror, 'same') + relative_blur = image / convolve_method(im_deconv, psf, 'same') + im_deconv *= convolve_method(relative_blur, psf_mirror, 'same') if clip: im_deconv[im_deconv > 1] = 1 From 3f3b705ebdbba42a430db96621778df4a882b836 Mon Sep 17 00:00:00 2001 From: scottsievert Date: Sat, 12 Dec 2015 19:26:10 -0600 Subject: [PATCH 3/6] fixes import, tuple addition --- skimage/restoration/deconvolution.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/skimage/restoration/deconvolution.py b/skimage/restoration/deconvolution.py index 6d922aa242e..182591276e2 100644 --- a/skimage/restoration/deconvolution.py +++ b/skimage/restoration/deconvolution.py @@ -7,7 +7,7 @@ import numpy as np import numpy.random as npr -from scipy.signal import fftconvolve +from scipy.signal import fftconvolve, convolve2d from . import uft @@ -369,8 +369,8 @@ def richardson_lucy(image, psf, iterations=50, clip=True): def fft_time(m, n, k, l): return m*np.log(m) + n*np.log(n) + k*np.log(k) + l*np.log(l) - time_ratio = 71.468 * fft_time(*image.shape, *psf.shape) - time_ratio /= direct_time(*image.shape, *psf.shape) + time_ratio = 71.468 * fft_time(*(image.shape + psf.shape)) + time_ratio /= direct_time(*(image.shape + psf.shape)) convolve_method = fftconvolve if time_ratio <= 1 else convolve2d image = image.astype(np.float) From 3fd47091dd85a8834a6fe5b2d51112d6e9815ade Mon Sep 17 00:00:00 2001 From: scottsievert Date: Sat, 12 Dec 2015 21:15:01 -0600 Subject: [PATCH 4/6] adds comment explaining time_ratio decision --- skimage/restoration/deconvolution.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/skimage/restoration/deconvolution.py b/skimage/restoration/deconvolution.py index 182591276e2..426766f7f72 100644 --- a/skimage/restoration/deconvolution.py +++ b/skimage/restoration/deconvolution.py @@ -369,6 +369,8 @@ def richardson_lucy(image, psf, iterations=50, clip=True): def fft_time(m, n, k, l): return m*np.log(m) + n*np.log(n) + k*np.log(k) + l*np.log(l) + # see whether the fourier transform convolution method or the direct + # convolution method is faster (discussed in scikit-image PR #1792) time_ratio = 71.468 * fft_time(*(image.shape + psf.shape)) time_ratio /= direct_time(*(image.shape + psf.shape)) convolve_method = fftconvolve if time_ratio <= 1 else convolve2d From 4d6c9e7b845e9e412bc94a02c7ab602823087605 Mon Sep 17 00:00:00 2001 From: scottsievert Date: Sun, 13 Dec 2015 14:57:52 -0600 Subject: [PATCH 5/6] now N dimensional, changes constant, cleans and comments --- skimage/restoration/deconvolution.py | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/skimage/restoration/deconvolution.py b/skimage/restoration/deconvolution.py index 426766f7f72..c70030c5acf 100644 --- a/skimage/restoration/deconvolution.py +++ b/skimage/restoration/deconvolution.py @@ -7,7 +7,7 @@ import numpy as np import numpy.random as npr -from scipy.signal import fftconvolve, convolve2d +from scipy.signal import fftconvolve, convolve from . import uft @@ -336,7 +336,7 @@ def richardson_lucy(image, psf, iterations=50, clip=True): Parameters ---------- image : ndarray - Input degraded image. + Input degraded image (can be N dimensional). psf : ndarray The point spread function. iterations : int @@ -365,15 +365,23 @@ def richardson_lucy(image, psf, iterations=50, clip=True): ---------- .. [1] http://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution """ - direct_time = lambda n, m, k, l: k*l * n*m - def fft_time(m, n, k, l): - return m*np.log(m) + n*np.log(n) + k*np.log(k) + l*np.log(l) + # compute the times for direct convolution and the fft method. The fft is of + # complexity O(N log(N)) for each dimension and the direct method does + # straight arithmetic (and is O(n*k) to add n elements k times) + def direct_time(img_shape, kernel_shape): + return np.prod(img_shape + kernel_shape) + def fft_time(img_shape, kernel_shape): + return np.sum([n*np.log(n) for n in img_shape+kernel_shape]) # see whether the fourier transform convolution method or the direct # convolution method is faster (discussed in scikit-image PR #1792) - time_ratio = 71.468 * fft_time(*(image.shape + psf.shape)) - time_ratio /= direct_time(*(image.shape + psf.shape)) - convolve_method = fftconvolve if time_ratio <= 1 else convolve2d + time_ratio = 40.032 * fft_time(image.shape, psf.shape)) + time_ratio /= direct_time(image.shape, psf.shape) + + if time_ratio <= 1 or len(image.shape) > 2: + convolve_method = fftconvolve + else: + convolve_method = convolve image = image.astype(np.float) psf = psf.astype(np.float) From 2892e90abc1afd5271a8f0633b4a4f81a2ef3cda Mon Sep 17 00:00:00 2001 From: scottsievert Date: Sun, 13 Dec 2015 16:59:59 -0600 Subject: [PATCH 6/6] removes paren --- skimage/restoration/deconvolution.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/restoration/deconvolution.py b/skimage/restoration/deconvolution.py index c70030c5acf..3cfeb269d14 100644 --- a/skimage/restoration/deconvolution.py +++ b/skimage/restoration/deconvolution.py @@ -375,7 +375,7 @@ def fft_time(img_shape, kernel_shape): # see whether the fourier transform convolution method or the direct # convolution method is faster (discussed in scikit-image PR #1792) - time_ratio = 40.032 * fft_time(image.shape, psf.shape)) + time_ratio = 40.032 * fft_time(image.shape, psf.shape) time_ratio /= direct_time(image.shape, psf.shape) if time_ratio <= 1 or len(image.shape) > 2: