scikit-image · alexjarosch · Oct 20, 2020 · Oct 20, 2020 · Oct 20, 2020 · Oct 20, 2020
diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
@@ -299,6 +299,9 @@ each and every contributor, regardless of their role.
 - Taylor D. Scott
   Simplified _upsampled_dft and extended register_translation to nD images.
 
+- Alexander H. Jarosch
+  Deconvolution algorithm based on natural priors.
+
 - David J. Mellert
   Polar and log-polar warping, nD windows
 

diff --git a/benchmarks/benchmark_restoration.py b/benchmarks/benchmark_restoration.py
@@ -97,6 +97,18 @@ def peakmem_richardson_lucy_f32(self):
         restoration.richardson_lucy(self.volume_f32, self.psf_f32,
                                     iterations=1)
 
+    def time_gaussian_natural_prior_f64(self):
+        restoration.gaussian_natural_prior(self.volume_f64, self.psf_f64)
+
+    def time_gaussian_natural_prior_f32(self):
+        restoration.gaussian_natural_prior(self.volume_f32, self.psf_f32)
+
+    def peakmem_gaussian_natural_prior_f64(self):
+        restoration.gaussian_natural_prior(self.volume_f64, self.psf_f64)
+
+    def peakmem_gaussian_natural_prior_f32(self):
+        restoration.gaussian_natural_prior(self.volume_f32, self.psf_f32)
+
 
 class RollingBall(object):
     """Benchmark Rolling Ball algorithm."""

diff --git a/doc/examples/filters/plot_deconvolution_natural_prior.py b/doc/examples/filters/plot_deconvolution_natural_prior.py
@@ -0,0 +1,59 @@
+"""
+==================================
+Image Deconvolution Natural Priors
+==================================
+In this example, we deconvolve an image using deconvolution with natural
+priors, e.g Gauss, solved in the frequency domain with a new
+deconvolution algorithm ([1]_, [2]_, [3]_).
+
+The algorithm is based on a PSF (Point Spread Function),
+where PSF is described as the impulse response of the
+optical system. The blurred image is sharpened with a smoothness weight,
+which defines the weight of the Gauss prior during deconvolution .
+
+.. [1] http://groups.csail.mit.edu/graphics/CodedAperture/
+.. [2] Levin, A., Fergus, R., Durand, F., & Freeman, W. T. (2007).
+       Deconvolution using natural image priors.
+       Massachusetts Institute of Technology,
+       Computer Science and Artificial Intelligence Laboratory, 3.
+.. [3] Levin, A., Fergus, R., Durand, F., & Freeman, W. T. (2007).
+       Image and depth from a conventional camera with a coded aperture.
+       ACM transactions on graphics (TOG), 26(3), 70-es.
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+<
8000
/span>
+from scipy.signal import convolve2d as conv2
+
+from skimage import color, data, restoration
+
+astro = color.rgb2gray(data.astronaut())
+
+psf = np.ones((5, 5)) / 25
+astro = conv2(astro, psf, 'same')
+# Add Noise to Image
+astro_noisy = astro.copy()
+astro_noisy += (np.random.poisson(lam=25, size=astro.shape) - 10) / 255.
+
+# Restore Image using Richardson-Lucy algorithm
+deconvolved_DNP = restoration.gaussian_natural_prior(astro_noisy, psf)
+
+fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(8, 5))
+plt.gray()
+
+for a in (ax[0], ax[1], ax[2]):
+    a.axis('off')
+
+ax[0].imshow(astro)
+ax[0].set_title('Original Data')
+
+ax[1].imshow(astro_noisy)
+ax[1].set_title('Noisy data')
+
+ax[2].imshow(deconvolved_DNP, vmin=astro_noisy.min(), vmax=astro_noisy.max())
+ax[2].set_title('Restoration using\nnatural Gauss prior')
+
+
+fig.subplots_adjust(wspace=0.02, hspace=0.2,
+                    top=0.9, bottom=0.05, left=0, right=1)
+plt.show()
diff --git a/skimage/data/_registry.py b/skimage/data/_registry.py
@@ -121,6 +121,8 @@
     "restoration/tests/astronaut_rl.npy":
         "3f8373e2c6182a89366e51cef6624e3625deac75fdda1079cbdad2a33322152c",
     "restoration/tests/camera_rl.npy": "fd4f59af84dd471fbbe79ee70c1b7e68a69864c461f0db5ac587e7975363f78f",
+    "restoration/tests/astronaut_dnp.npy":
+        "5a9d319ab7749c369f5a670d684b9d56d619b8c8e6a326b14879baa4bba072f7",
     "restoration/tests/camera_unsup.npy": "401e512f7becec46c3379f34dbd7956a02a4f2e68f802a8d5bd605d355588880",
     "restoration/tests/camera_unsup2.npy": "cc4b31ea2a6ebe07bc7529174677b3e94eae4b7679fd822d3e721a9c8d4854d1",
     "restoration/tests/camera_wiener.npy": "4505ea8b0d63d03250c6d756560d615751b76dd6ffc4a95972fa260c0c84633e",
@@ -152,6 +154,8 @@
     "data/pivchallenge-B-B001_1.tif": "https://gitlab.com/scikit-image/data/-/raw/master/pivchallenge/B/B001_1.tif",
     "data/pivchallenge-B-B001_2.tif": "https://gitlab.com/scikit-image/data/-/raw/master/pivchallenge/B/B001_2.tif",
     "restoration/tests/astronaut_rl.npy": "https://gitlab.com/scikit-image/data/-/raw/master/astronaut_rl.npy",
+    "restoration/tests/astronaut_dnp.npy":
+        "https://gitlab.com/scikit-image/data/-/raw/master/astronaut_dnp.npy",
 }
 
 legacy_registry = {

diff --git a/skimage/restoration/__init__.py b/skimage/restoration/__init__.py
@@ -2,7 +2,8 @@
 
 """
 
-from .deconvolution import wiener, unsupervised_wiener, richardson_lucy
+from .deconvolution import (wiener, unsupervised_wiener, richardson_lucy,
+                            gaussian_natural_prior)
 from .unwrap import unwrap_phase
 from ._denoise import (denoise_tv_chambolle, denoise_tv_bregman,
                        denoise_bilateral, denoise_wavelet, estimate_sigma)
@@ -16,6 +17,7 @@
 __all__ = ['wiener',
            'unsupervised_wiener',
            'richardson_lucy',
+           'gaussian_natural_prior',
            'unwrap_phase',
            'denoise_tv_bregman',
            'denoise_tv_chambolle',

diff --git a/skimage/restoration/deconvolution.py b/skimage/restoration/deconvolution.py
@@ -6,6 +6,8 @@
 from scipy.signal import convolve
 
 from . import uft
+from .._shared.fft import fftmodule as fft
+from ..util import crop
 
 __keywords__ = "restoration, image, deconvolution"
 
@@ -383,3 +385,105 @@ def richardson_lucy(image, psf, iterations=50, clip=True, filter_epsilon=None):
         im_deconv[im_deconv < -1] = -1
 
     return im_deconv
+
+
+def gaussian_natural_prior(image, psf, smoothness_weight=0.01, clip=True,
+                           pad_width=0):
+    """ Deconvolution using Gaussian natural image priors.
+
+    Parameters
+    ----------
+    image : ndarray
+        Input degraded image.
+    psf : ndarray
+        The point spread function.
+    smoothness_weight : float, optional
+        The smoothness weight defines how much weight the Gaussian prior is
+        given in the process. No prior is used when smoothness_weight=0.
+    clip : boolean, optional
+        True by default. If true, pixel value of the result above 1 or
+        under 0 are thresholded for skimage pipeline compatibility.
+    pad_width : int or tuple of int, optional
+        If pad_width > 0, the image will be extended by `pad_width` along each
+        boundary by use of `numpy.pad` with `mode='reflect'`. This is done to
+        reduce artifacts that can arise near the edges of the deconvolved image
+        due to the periodic nature of the discrete Fourier transform.
+
+    Returns
+    -------
+    im_deconv : ndarray
+        The deconvolved image.
+
+    Examples
+    --------
+    >>> from skimage import img_as_float, data, restoration
+    >>> from scipy.ndimage import convolve
+    >>> camera = img_as_float(data.camera())
+    >>> psf = np.ones((5, 5)) / 25
+    >>> camera = convolve(camera, psf)
+    >>> camera += 0.1 * camera.std() * np.random.standard_normal(camera.shape)
+    >>> deconvolved = restoration.gaussian_natural_prior(camera, psf)
+
+    Notes
+    -----
+    This algorithm solves the deconvolution in the frequency domain. Artefacts
+    like 'ringing' on sharp transitions are strongly reduced in comparison to
+    other deconvolution strategies.
+    The smoothness weight needs to be estimated by trial and error. See [1-3]
+    for guidance.
+
+    References
+    ----------
+    .. [1] http://groups.csail.mit.edu/graphics/CodedAperture/
+    .. [2] Levin, A., Fergus, R., Durand, F., & Freeman, W. T. (2007).
+           Deconvolution using natural image priors.
+           Massachusetts Institute of Technology,
+           Computer Science and Artificial Intelligence Laboratory, 3.
+    .. [3] Levin, A., Fergus, R., Durand, F., & Freeman, W. T. (2007).
+           Image and depth from a conventional camera with a coded aperture.
+           ACM transactions on graphics (TOG), 26(3), 70-es.
+    """
+    float_type = np.promote_types(image.dtype, np.float32)
+    image = image.astype(float_type, copy=False)
+    psf = psf.astype(float_type, copy=False)
+    if image.ndim != psf.ndim:
+        raise ValueError(
+            "psf and image must have an equal number of dimensions"
+        )
+
+    # pad to reduce boundary artifacts
+    image = np.pad(image, pad_width, mode='reflect')
+    shape = image.shape
+
+    # sum of squared first order differences along each axis
+    ndim = psf.ndim
+    G = 0
+    d_shape = [1, ] * ndim
+    for n in range(ndim):
+        if n == ndim - 1:
+            D = fft.rfft([-1, 1], n=shape[n])
+        else:
+            D = fft.fft([-1, 1], n=shape[n])
+        d_shape[n] = D.size
+        D = D.reshape(d_shape)  # reshape 1D array for broadcasting
+        G = G + np.conj(D) * D
+        d_shape[n] = 1
+
+    F = fft.rfftn(psf, s=shape)
+    F_conj = np.conj(F)
+    A = F_conj * F + smoothness_weight * G
+    b = F_conj * fft.rfft2(image)
+
+    X = np.divide(b, A)
+    x = np.real(fft.irfft2(X))
+
+    # recenter the deconvolved image
+    shift = tuple([(s - 1) // 2 for s in psf.shape])
+    im_deconv = np.roll(x, shift=shift, axis=range(psf.ndim))
+    # remove any padding that was previously added
+    im_deconv = crop(im_deconv, pad_width)
+    if clip:
+        im_deconv[im_deconv > 1] = 1
+        im_deconv[im_deconv < 0] = 0
+
+    return im_deconv
diff --git a/skimage/restoration/tests/test_restoration.py b/skimage/restoration/tests/test_restoration.py
@@ -7,6 +7,7 @@
 from skimage.data import astronaut, camera
 from skimage import restoration, util
 from skimage.restoration import uft
+from skimage.metrics import structural_similarity
 
 test_img = util.img_as_float(camera())
 
@@ -62,7 +63,7 @@ def test_image_shape():
     point[2, 2] = 1.
     psf = ndi.gaussian_filter(point, sigma=1.)
     # image shape: (45, 45), as reported in #1172
-    image = util.img_as_float(camera()[65:165, 215:315]) # just the face
+    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]
@@ -87,6 +88,26 @@ def test_richardson_lucy():
     np.testing.assert_allclose(deconvolved, np.load(path), rtol=1e-3)
 
 
+def test_gaussian_natural_prior():
+    ideal = rgb2gray(astronaut())
+    psf = np.ones((5, 5)) / 25
+    data = convolve2d(ideal, psf, 'same')
+    # add some noise to the data
+    np.random.seed(0)
+    data_n = data + 0.1 * data.std() * np.random.standard_normal(data.shape)
+    deconvolved_noise = restoration.gaussian_natural_prior(data_n, psf)
+    # also deconvolve with no noise data
+    deconvolved = restoration.gaussian_natural_prior(data, psf)
+
+    # test for improvement of the image on noise free data
+    ssim_blurred = structural_similarity(ideal, data)
+    ssim_deconv = structural_similarity(ideal, deconvolved)
+    assert ssim_deconv > ssim_blurred
+
+    path = fetch('restoration/tests/astronaut_dnp.npy')
+    np.testing.assert_allclose(deconvolved_noise, np.load(path), rtol=1e-5)
+
+
 @pytest.mark.parametrize('dtype_image', [np.float32, np.float64])
 @pytest.mark.parametrize('dtype_psf', [np.float32, np.float64])
 def test_richardson_lucy_filtered(dtype_image, dtype_psf):