Merge pull request #1804 from soupault/inpainting

ENH: Inpainting with biharmonic equation
2026-07-03 09:45:23 +08:00 · 2016-01-18 10:56:24 +11:00
parent b7fa0acd86 3ae2f1236b
commit e53d7daaa3
3 changed files with 263 additions and 0 deletions
@@ -0,0 +1,58 @@
+"""
+===========
+Inpainting
+===========
+Inpainting [1]_ is the process of reconstructing lost or deteriorated
+parts of images and videos. 
+
+The reconstruction is supposed to be performed in fully automatic way by
+exploiting the information presented in non-damaged regions.
+
+In this example, we show how the masked pixels get inpainted by
+inpainting algorithm based on 'biharmonic equation'-assumption [2]_ [3]_.
+
+.. [1]  Wikipedia. Inpainting
+        https://en.wikipedia.org/wiki/Inpainting
+.. [2]  Wikipedia. Biharmonic equation
+        https://en.wikipedia.org/wiki/Biharmonic_equation
+.. [3]  N.S.Hoang, S.B.Damelin, "On surface completion and image 
+        inpainting by biharmonic functions: numerical aspects",
+        http://www.ima.umn.edu/~damelin/biharmonic
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from skimage import data, color
+from skimage.restoration import inpaint
+
+image_orig = data.astronaut()
+
+# Create mask with three defect regions: left, middle, right respectively
+mask = np.zeros(image_orig.shape[:-1])
+mask[20:60, 0:20] = 1
+mask[200:300, 150:170] = 1
+mask[50:100, 400:430] = 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
+
+image_result = inpaint.inpaint_biharmonic(image_defect, mask, multichannel=True)
+
+fig, axes = plt.subplots(ncols=3, nrows=1)
+
+axes[0].set_title('Defected image')
+axes[0].imshow(image_orig)
+axes[0].set_xticks([]), axes[0].set_yticks([])
+
+axes[1].set_title('Defect mask')
+axes[1].imshow(mask, cmap=plt.cm.gray)
+axes[1].set_xticks([]), axes[1].set_yticks([])
+
+axes[2].set_title('Inpainted image')
+axes[2].imshow(image_result)
+axes[2].set_xticks([]), axes[2].set_yticks([])
+
+plt.show()
@@ -0,0 +1,139 @@
+from __future__ import division
+
+import numpy as np
+import skimage
+from scipy import sparse
+from scipy.sparse.linalg import spsolve
+from scipy.ndimage.filters import laplace
+
+
+def _get_neighborhood(nd_idx, radius, nd_shape):
+    bounds_lo = (nd_idx - radius).clip(min=0)
+    bounds_hi = (nd_idx + radius + 1).clip(max=nd_shape)
+    return bounds_lo, bounds_hi
+
+
+def _inpaint_biharmonic_single_channel(img, mask, out, limits):
+    # Initialize sparse matrices
+    matrix_unknown = sparse.lil_matrix((np.sum(mask), out.size))
+    matrix_known = sparse.lil_matrix((np.sum(mask), out.size))
+
+    # Find indexes of masked points in flatten array
+    mask_i = np.ravel_multi_index(np.where(mask), mask.shape)
+
+    # Find masked points and prepare them to be easily enumerate over
+    mask_pts = np.array(np.where(mask)).T
+
+    # Iterate over masked points
+    for mask_pt_n, mask_pt_idx in enumerate(mask_pts):
+        # Get bounded neighborhood of selected radius
+        b_lo, b_hi = _get_neighborhood(mask_pt_idx, 2, out.shape)
+
+        # Create biharmonic coefficients ndarray
+        neigh_coef = np.zeros(b_hi - b_lo)
+        neigh_coef[tuple(mask_pt_idx - b_lo)] = 1
+        neigh_coef = laplace(laplace(neigh_coef))
+
+        # Iterate over masked point's neighborhood
+        it_inner = np.nditer(neigh_coef, flags=['multi_index'])
+        for coef in it_inner:
+            if coef == 0:
+                continue
+            tmp_pt_idx = np.add(b_lo, it_inner.multi_index)
+            tmp_pt_i = np.ravel_multi_index(tmp_pt_idx, mask.shape)
+
+            if mask[tuple(tmp_pt_idx)]:
+                matrix_unknown[mask_pt_n, tmp_pt_i] = coef
+            else:
+                matrix_known[mask_pt_n, tmp_pt_i] = coef
+
+    # Prepare diagonal matrix
+    flat_diag_image = sparse.dia_matrix((out.flatten(), np.array([0])),
+                                        shape=(out.size, out.size))
+
+    # Calculate right hand side as a sum of known matrix's columns
+    matrix_known = matrix_known.tocsr()
+    rhs = -(matrix_known * flat_diag_image).sum(axis=1)
+
+    # Solve linear system for masked points
+    matrix_unknown = matrix_unknown[:, mask_i]
+    matrix_unknown = sparse.csr_matrix(matrix_unknown)
+    result = spsolve(matrix_unknown, rhs)
+
+    # Handle enormous values
+    result = np.clip(result, *limits)
+
+    result = result.ravel()
+
+    # Substitute masked points with inpainted versions
+    for mask_pt_n, mask_pt_idx in enumerate(mask_pts):
+        out[tuple(mask_pt_idx)] = result[mask_pt_n]
+
+    return out
+
+
+def inpaint_biharmonic(img, mask, multichannel=False):
+    """Inpaint masked points in image with biharmonic equations.
+
+    Parameters
+    ----------
+    img : (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 'img' channels. Unknown pixels have to be represented with 1,
+        known pixels - with 0.
+    multichannel : boolean, optional
+        If True, the last `img` dimension is considered as a color channel,
+        otherwise as spatial.
+
+    Returns
+    -------
+    out : (M[, N[, ..., P]][, C]) ndarray
+        Input image with masked pixels inpainted.
+
+    Example
+    -------
+    >>> 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)
+
+    References
+    ----------
+    Algorithm is based on:
+    .. [1]  N.S.Hoang, S.B.Damelin, "On surface completion and image inpainting
+            by biharmonic functions: numerical aspects",
+            http://www.ima.umn.edu/~damelin/biharmonic
+    """
+
+    if img.ndim < 1:
+        raise ValueError('Input array has to be at least 1D')
+    
+    img_baseshape = img.shape[:-1] if multichannel else img.shape
+    if img_baseshape != mask.shape:
+        raise ValueError('Input arrays have to be the same shape')
+    
+    if np.ma.isMaskedArray(img):
+        raise TypeError('Masked arrays are not supported')
+
+    img = skimage.img_as_float(img)
+    mask = mask.astype(np.bool)
+
+    if not multichannel:
+        img = img[..., np.newaxis]
+
+    out = np.copy(img)
+
+    for i in range(img.shape[-1]):
+        known_points = img[..., i][~mask]
+        limits = (np.min(known_points), np.max(known_points))
+        _inpaint_biharmonic_single_channel(img[..., i], mask,
+                                           out[..., i], limits)
+
+    if not multichannel:
+        out = out[..., 0]
+
+    return out
@@ -0,0 +1,66 @@
+from __future__ import print_function, division
+
+import numpy as np
+from numpy.testing import (run_module_suite, assert_allclose,
+                           assert_raises)
+from skimage.restoration import inpaint
+
+
+def test_inpaint_biharmonic_2d():
+    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
+    out = inpaint.inpaint_biharmonic(img, mask)
+    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)
+
+
+def test_inpaint_biharmonic_3d():
+    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)
+    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([])
+    assert_raises(ValueError, inpaint.inpaint_biharmonic, img, mask)
+
+    img, mask = np.zeros((2, 2)), np.zeros((4, 1))
+    assert_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))
+    assert_raises(TypeError, inpaint.inpaint_biharmonic, img, mask)
+
+
+if __name__ == '__main__':
+    run_module_suite()