diff --git a/doc/examples/filters/plot_inpaint.py b/doc/examples/filters/plot_inpaint.py new file mode 100644 index 00000000..104694ec --- /dev/null +++ b/doc/examples/filters/plot_inpaint.py @@ -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() diff --git a/skimage/restoration/inpaint.py b/skimage/restoration/inpaint.py new file mode 100644 index 00000000..bf696e04 --- /dev/null +++ b/skimage/restoration/inpaint.py @@ -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 diff --git a/skimage/restoration/tests/test_inpaint.py b/skimage/restoration/tests/test_inpaint.py new file mode 100644 index 00000000..04e8c8c2 --- /dev/null +++ b/skimage/restoration/tests/test_inpaint.py @@ -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()