From 1f3721fcbe72698ff47fd436b7fb82c8138b6644 Mon Sep 17 00:00:00 2001
From: Egor Panfilov <multicolor.mood@gmail.com>
Date: Sun, 6 Dec 2015 15:52:45 +0300
Subject: [PATCH 01/13] ENH: Inpainting with biharmonic equation

---
 skimage/restoration/inpaint.py            | 167 ++++++++++++++++++++++
 skimage/restoration/tests/test_inpaint.py |  37 +++++
 2 files changed, 204 insertions(+)
 create mode 100644 skimage/restoration/inpaint.py
 create mode 100644 skimage/restoration/tests/test_inpaint.py

diff --git a/skimage/restoration/inpaint.py b/skimage/restoration/inpaint.py
new file mode 100644
index 00000000..602e1062
--- /dev/null
+++ b/skimage/restoration/inpaint.py
@@ -0,0 +1,167 @@
+from __future__ import print_function, division
+
+import numpy as np
+import skimage
+from scipy import sparse
+from scipy.sparse.linalg import spsolve
+
+
+def inpaint_biharmonic(img, mask):
+    """Inpaint masked points in image with biharmonic equations.
+
+    Parameters
+    ----------
+    img : 2-D np.array
+        Input image.
+    mask : 2-D np.array
+        Array of pixels to be inpainted. Have to be the same size as 'img'.
+        Unknown pixels has to be represented with 1, known pixels - with 0.
+
+    Returns
+    -------
+    out : 2-D np.array
+        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
+
+    Realization is based on:
+    .. [2]  John D'Errico,
+        http://www.mathworks.com/matlabcentral/fileexchange/4551-inpaint-nans,
+            method 3
+    """
+
+    if img.ndim != 2 or mask.ndim != 2:
+        raise ValueError('Only 2-dimensional arrays are supported')
+    if img.shape != mask.shape:
+        raise ValueError('Input arrays have to be the same shape')
+    if np.ma.isMaskedArray(img):
+        raise TypeError('Masked arrays are not supported')
+
+    # TODO: add sufficient conditions (if any)
+
+    img = skimage.img_as_float(img)
+    mask = mask.astype(np.bool)
+
+    out = np.copy(img)
+    out_h, out_w = out.shape
+    out_l = out.size
+
+    def _in_bounds(idx):
+        if len(idx) == 1:
+            return 0 <= idx <= out_l - 1
+        else:
+            return (0 <= idx[0] <= out_h - 1) and (0 <= idx[1] <= out_w - 1)
+
+    # Find indexes of masked points in flatten array
+    mask_mn = np.array(np.where(mask)).T
+    mask_i = np.ravel_multi_index(np.where(mask), mask.shape)
+
+    # Initialize sparse matrix
+    # TODO: only points required for computation might be considered
+    matrix_unknown = sparse.lil_matrix((np.sum(mask), out.size), dtype=np.int32)
+    matrix_known = sparse.lil_matrix((np.sum(mask), out.size), dtype=np.int32)
+
+    # INFO: kernels can be reworked using scipy.signal.convolve2d
+    #       and np.array([0, 1, 0], [1, -4, 1], [0, 1, 0])
+
+    # 1 stage. Find points 2 or more pixels far from bounds
+    kernel = [1, 2, -8, 2, 1, -8, 20, -8, 1, 2, -8, 2, 1]
+    offset = [-2 * out_w, -out_w - 1, -out_w, -out_w + 1,
+              -2, -1, 0, 1, 2, out_w - 1, out_w, out_w + 1, 2 * out_w]
+
+    for idx, (i, (m, n)) in enumerate(zip(mask_i, mask_mn)):
+        if 2 <= m <= out_h - 3 and 2 <= n <= out_w - 3:
+            for k, o in zip(kernel, offset):
+                if i + o in mask_i:
+                    matrix_unknown[idx, i + o] = k
+                else:
+                    matrix_known[idx, i + o] = k
+
+    # 2 stage. Find points 1 pixel far from bounds
+    kernel = [1, 1, -4, 1, 1]
+    offset = [-out_w, -1, 0, 1, out_w]
+
+    for idx, (i, (m, n)) in enumerate(zip(mask_i, mask_mn)):
+        if (m in [1, out_h - 2] and 1 <= n <= out_h - 2) or \
+           (n in [1, out_w - 2] and 1 <= m <= out_w - 2):
+            for k, o in zip(kernel, offset):
+                if i + o in mask_i:
+                    matrix_unknown[idx, i + o] = k
+                else:
+                    matrix_known[idx, i + o] = k
+
+    # 3 stage. Find points on the horizontal bounds
+    kernel = [1, -2, 1]
+    offset = [-1, 0, 1]
+
+    for idx, (i, (m, n)) in enumerate(zip(mask_i, mask_mn)):
+        if m in [0, out_h - 1] and 1 <= n <= out_w - 2:
+            for k, o in zip(kernel, offset):
+                if i + o in mask_i:
+                    matrix_unknown[idx, i + o] = k
+                else:
+                    matrix_known[idx, i + o] = k
+
+    # 4 stage. Find points on the vertical bounds
+    kernel = [1, -2, 1]
+    offset = [-out_w, 0, out_w]
+
+    for idx, (i, (m, n)) in enumerate(zip(mask_i, mask_mn)):
+        if n in [0, out_w - 1] and 1 <= m <= out_h - 2:
+            for k, o in zip(kernel, offset):
+                if i + o in mask_i:
+                    matrix_unknown[idx, i + o] = k
+                else:
+                    matrix_known[idx, i + o] = k
+
+    # 5 stage. Find corner points if any
+    kernel = [1, 1, -2, 1, 1]
+    offset = [-out_w, -1, 0, 1, out_w]
+    offset_mn = [(-1, 0), (0, -1), (0, 0), (0, 1), (1, 0)]
+
+    for idx, (i, (m, n)) in enumerate(zip(mask_i, mask_mn)):
+        if m in [0, out_h - 1] and n in [0, out_w - 1]:
+            for k, o_mn in zip(kernel, offset_mn):
+                if _in_bounds((m + o_mn[0], n + o_mn[1])):
+                    o = offset[offset_mn.index(o_mn)]
+                    if i + o in mask_i:
+                        matrix_unknown[idx, i + o] = k
+                    else:
+                        matrix_known[idx, i + o] = k
+
+    # 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 columns
+    matrix_known = matrix_known.tocsr()
+    rhs = -(matrix_known * flat_diag_image).sum(axis=1)
+
+    # Solve linear system over defect points
+    matrix_unknown = matrix_unknown[:, mask_i]
+    matrix_unknown = sparse.csr_matrix(matrix_unknown)
+    result = spsolve(matrix_unknown, rhs)
+
+    # Handle enormous values
+    result[np.where(result < -1)] = -1
+    result[np.where(result > 1)] = 1
+
+    # Put calculated points into the image
+    for idx, (m, n) in enumerate(mask_mn):
+        out[m, n] = result[idx]
+
+    return out
diff --git a/skimage/restoration/tests/test_inpaint.py b/skimage/restoration/tests/test_inpaint.py
new file mode 100644
index 00000000..106c848c
--- /dev/null
+++ b/skimage/restoration/tests/test_inpaint.py
@@ -0,0 +1,37 @@
+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():
+    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.inpaint_biharmonic(img, mask)
+    ref = [[0., 0.0625, 0.25, 0.5625, 0.671875],
+           [0., 0.0625, 0.25, 0.5390625, 0.78125],
+           [0., 0.0625, 0.2578125, 0.5625, 0.890625],
+           [0., 0.0625, 0.25, 0.5625, 1.],
+           [0., 0.0625, 0.25, 0.5625, 1.]]
+    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()

From dc05d56c8921e844ee38e8a351991f3b4c23ecb7 Mon Sep 17 00:00:00 2001
From: Egor Panfilov <multicolor.mood@gmail.com>
Date: Tue, 8 Dec 2015 00:56:32 +0300
Subject: [PATCH 02/13] DOC: Added gallery example for inpainting

---
 doc/examples/plot_inpaint.py | 58 ++++++++++++++++++++++++++++++++++++
 1 file changed, 58 insertions(+)
 create mode 100644 doc/examples/plot_inpaint.py

diff --git a/doc/examples/plot_inpaint.py b/doc/examples/plot_inpaint.py
new file mode 100644
index 00000000..ed79722b
--- /dev/null
+++ b/doc/examples/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]_.
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from skimage import data, color
+from skimage.restoration import inpaint
+
+image_orig = color.rgb2gray(data.astronaut())
+
+# Create mask with three defect regions: left, middle, right respectively
+mask = np.zeros_like(image_orig)
+mask[20:60, 0:20] = 1
+mask[200:300, 150:170] = 1
+mask[50:100, 400:430] = 1
+
+image_defect = image_orig.copy()
+image_defect[np.where(mask)] = 0
+
+image_result = inpaint.inpaint_biharmonic(image_defect, mask)
+
+fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, nrows=1, figsize=(10, 6))
+
+ax1.set_title('Defected image')
+ax1.imshow(image_orig, cmap=plt.cm.gray, interpolation='nearest')
+ax1.set_xticks([]), ax1.set_yticks([])
+
+ax2.set_title('Defect mask')
+ax2.imshow(mask, cmap=plt.cm.gray, interpolation='nearest')
+ax2.set_xticks([]), ax2.set_yticks([])
+
+ax3.set_title('Inpainted image')
+ax3.imshow(image_result, cmap=plt.cm.gray, interpolation='nearest')
+ax3.set_xticks([]), ax3.set_yticks([])
+
+plt.show()
+
+"""
+.. [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
+"""

From add5b8b4d6b4309429086fe7e04068e45e68fb9b Mon Sep 17 00:00:00 2001
From: Egor <multicolor.mood@gmail.com>
Date: Tue, 8 Dec 2015 13:39:42 +0300
Subject: [PATCH 03/13] FIX/TST: Fixed bug with defected bounding frame, added
 testcase

---
 skimage/restoration/inpaint.py            | 42 +++++++++--------------
 skimage/restoration/tests/test_inpaint.py | 28 ++++++++++++---
 2 files changed, 39 insertions(+), 31 deletions(-)

diff --git a/skimage/restoration/inpaint.py b/skimage/restoration/inpaint.py
index 602e1062..813ad6bb 100644
--- a/skimage/restoration/inpaint.py
+++ b/skimage/restoration/inpaint.py
@@ -51,8 +51,6 @@ def inpaint_biharmonic(img, mask):
     if np.ma.isMaskedArray(img):
         raise TypeError('Masked arrays are not supported')
 
-    # TODO: add sufficient conditions (if any)
-
     img = skimage.img_as_float(img)
     mask = mask.astype(np.bool)
 
@@ -78,7 +76,7 @@ def inpaint_biharmonic(img, mask):
     # INFO: kernels can be reworked using scipy.signal.convolve2d
     #       and np.array([0, 1, 0], [1, -4, 1], [0, 1, 0])
 
-    # 1 stage. Find points 2 or more pixels far from bounds
+    # 1 stage. Find points 2 or more pixels far from edges
     kernel = [1, 2, -8, 2, 1, -8, 20, -8, 1, 2, -8, 2, 1]
     offset = [-2 * out_w, -out_w - 1, -out_w, -out_w + 1,
               -2, -1, 0, 1, 2, out_w - 1, out_w, out_w + 1, 2 * out_w]
@@ -91,7 +89,7 @@ def inpaint_biharmonic(img, mask):
                 else:
                     matrix_known[idx, i + o] = k
 
-    # 2 stage. Find points 1 pixel far from bounds
+    # 2 stage. Find points 1 pixel far from edges
     kernel = [1, 1, -4, 1, 1]
     offset = [-out_w, -1, 0, 1, out_w]
 
@@ -104,31 +102,23 @@ def inpaint_biharmonic(img, mask):
                 else:
                     matrix_known[idx, i + o] = k
 
-    # 3 stage. Find points on the horizontal bounds
-    kernel = [1, -2, 1]
-    offset = [-1, 0, 1]
+    # 3 stage. Find points on the edges
+    kernel = [1, 1, -3, 1, 1]
+    offset = [-out_w, -1, 0, 1, out_w]
+    offset_mn = [(-1, 0), (0, -1), (0, 0), (0, 1), (1, 0)]
 
     for idx, (i, (m, n)) in enumerate(zip(mask_i, mask_mn)):
-        if m in [0, out_h - 1] and 1 <= n <= out_w - 2:
-            for k, o in zip(kernel, offset):
-                if i + o in mask_i:
-                    matrix_unknown[idx, i + o] = k
-                else:
-                    matrix_known[idx, i + o] = k
+        if (m in [0, out_h - 1] and 1 <= n <= out_w - 2) or \
+           (n in [0, out_w - 1] and 1 <= m <= out_h - 2):
+            for k, o_mn in zip(kernel, offset_mn):
+                if _in_bounds((m + o_mn[0], n + o_mn[1])):
+                    o = offset[offset_mn.index(o_mn)]
+                    if i + o in mask_i:
+                        matrix_unknown[idx, i + o] = k
+                    else:
+                        matrix_known[idx, i + o] = k
 
-    # 4 stage. Find points on the vertical bounds
-    kernel = [1, -2, 1]
-    offset = [-out_w, 0, out_w]
-
-    for idx, (i, (m, n)) in enumerate(zip(mask_i, mask_mn)):
-        if n in [0, out_w - 1] and 1 <= m <= out_h - 2:
-            for k, o in zip(kernel, offset):
-                if i + o in mask_i:
-                    matrix_unknown[idx, i + o] = k
-                else:
-                    matrix_known[idx, i + o] = k
-
-    # 5 stage. Find corner points if any
+    # 4 stage. Find corner points
     kernel = [1, 1, -2, 1, 1]
     offset = [-out_w, -1, 0, 1, out_w]
     offset_mn = [(-1, 0), (0, -1), (0, 0), (0, 1), (1, 0)]
diff --git a/skimage/restoration/tests/test_inpaint.py b/skimage/restoration/tests/test_inpaint.py
index 106c848c..39d7f0a7 100644
--- a/skimage/restoration/tests/test_inpaint.py
+++ b/skimage/restoration/tests/test_inpaint.py
@@ -13,11 +13,29 @@ def test_inpaint_biharmonic():
     mask[1, 3:] = 1
     mask[0, 4:] = 1
     out = inpaint.inpaint_biharmonic(img, mask)
-    ref = [[0., 0.0625, 0.25, 0.5625, 0.671875],
-           [0., 0.0625, 0.25, 0.5390625, 0.78125],
-           [0., 0.0625, 0.2578125, 0.5625, 0.890625],
-           [0., 0.0625, 0.25, 0.5625, 1.],
-           [0., 0.0625, 0.25, 0.5625, 1.]]
+    ref = np.array(
+        [[0., 0.0625, 0.25, 0.5625, 0.671875],
+         [0., 0.0625, 0.25, 0.5390625, 0.78125],
+         [0., 0.0625, 0.2578125, 0.5625, 0.890625],
+         [0., 0.0625, 0.25, 0.5625, 1.],
+         [0., 0.0625, 0.25, 0.5625, 1.]]
+    )
+    assert_allclose(ref, out)
+
+
+def test_inpaint_edges():
+    img = np.tile(np.square(np.linspace(0, 1, 5)), (5, 1))
+    mask = np.zeros_like(img)
+    mask[[0, -1], :] = 1
+    mask[:, [0, -1]] = 1
+    out = inpaint.inpaint_biharmonic(img, mask)
+    ref = np.array(
+        [[0.12199519, 0.15599245, 0.28348214, 0.44445398, 0.48737981],
+         [0.08799794, 0.0625, 0.25, 0.5625, 0.53030563],
+         [0.07949863, 0.0625, 0.25, 0.5625, 0.54103709],
+         [0.08799794, 0.0625, 0.25, 0.5625, 0.53030563],
+         [0.12199519, 0.15599245, 0.28348214, 0.44445398, 0.48737981]]
+    )
     assert_allclose(ref, out)
 
 

From 92f3114c9ea8c9d3433bb49a6f5cde9b67223aa0 Mon Sep 17 00:00:00 2001
From: Egor <multicolor.mood@gmail.com>
Date: Tue, 8 Dec 2015 14:05:11 +0300
Subject: [PATCH 04/13] FIX: Fixed test for biharmonic in accordance to new
 kernel for edges

---
 skimage/restoration/tests/test_inpaint.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/skimage/restoration/tests/test_inpaint.py b/skimage/restoration/tests/test_inpaint.py
index 39d7f0a7..1fe3e766 100644
--- a/skimage/restoration/tests/test_inpaint.py
+++ b/skimage/restoration/tests/test_inpaint.py
@@ -14,9 +14,9 @@ def test_inpaint_biharmonic():
     mask[0, 4:] = 1
     out = inpaint.inpaint_biharmonic(img, mask)
     ref = np.array(
-        [[0., 0.0625, 0.25, 0.5625, 0.671875],
-         [0., 0.0625, 0.25, 0.5390625, 0.78125],
-         [0., 0.0625, 0.2578125, 0.5625, 0.890625],
+        [[0., 0.0625, 0.25, 0.5625, 0.56947314],
+         [0., 0.0625, 0.25, 0.47029959, 0.57644628],
+         [0., 0.0625, 0.24664256, 0.49225207, 0.68956612],
          [0., 0.0625, 0.25, 0.5625, 1.],
          [0., 0.0625, 0.25, 0.5625, 1.]]
     )

From 0e6e08d0c09f3e80050e099c2d1ff72ba349c05b Mon Sep 17 00:00:00 2001
From: Egor Panfilov <multicolor.mood@gmail.com>
Date: Sun, 13 Dec 2015 12:42:48 +0300
Subject: [PATCH 05/13] TST: Attempt to fix Travis bug with old packages

---
 skimage/restoration/inpaint.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/skimage/restoration/inpaint.py b/skimage/restoration/inpaint.py
index 813ad6bb..b2c6ee9d 100644
--- a/skimage/restoration/inpaint.py
+++ b/skimage/restoration/inpaint.py
@@ -150,6 +150,8 @@ def inpaint_biharmonic(img, mask):
     result[np.where(result < -1)] = -1
     result[np.where(result > 1)] = 1
 
+    result = result.ravel()
+    
     # Put calculated points into the image
     for idx, (m, n) in enumerate(mask_mn):
         out[m, n] = result[idx]

From 1181c0dce840cb21f58f73220e4ddbd9dfbc37a2 Mon Sep 17 00:00:00 2001
From: Egor Panfilov <multicolor.mood@gmail.com>
Date: Sun, 13 Dec 2015 13:06:27 +0300
Subject: [PATCH 06/13] DOC: Moved references to common docstring

---
 doc/examples/plot_inpaint.py | 18 ++++++++----------
 1 file changed, 8 insertions(+), 10 deletions(-)

diff --git a/doc/examples/plot_inpaint.py b/doc/examples/plot_inpaint.py
index ed79722b..7ccc73cc 100644
--- a/doc/examples/plot_inpaint.py
+++ b/doc/examples/plot_inpaint.py
@@ -10,6 +10,14 @@ 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
@@ -46,13 +54,3 @@ ax3.imshow(image_result, cmap=plt.cm.gray, interpolation='nearest')
 ax3.set_xticks([]), ax3.set_yticks([])
 
 plt.show()
-
-"""
-.. [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
-"""

From 2fcf2eafe7eb7c6d81247d1d31db718912cf3a35 Mon Sep 17 00:00:00 2001
From: Egor Panfilov <multicolor.mood@gmail.com>
Date: Thu, 24 Dec 2015 01:31:31 +0300
Subject: [PATCH 07/13] ENH: Added support for multichannel arrays, fixed float
 scale to [0;1]

---
 skimage/restoration/inpaint.py | 234 ++++++++++++++++++---------------
 1 file changed, 127 insertions(+), 107 deletions(-)

diff --git a/skimage/restoration/inpaint.py b/skimage/restoration/inpaint.py
index b2c6ee9d..e41c68b9 100644
--- a/skimage/restoration/inpaint.py
+++ b/skimage/restoration/inpaint.py
@@ -6,16 +6,18 @@ from scipy import sparse
 from scipy.sparse.linalg import spsolve
 
 
-def inpaint_biharmonic(img, mask):
+def inpaint_biharmonic(img, mask, multichannel=False):
     """Inpaint masked points in image with biharmonic equations.
 
     Parameters
     ----------
-    img : 2-D np.array
+    img : 2D{+color} np.array
         Input image.
-    mask : 2-D np.array
+    mask : 2D np.array
         Array of pixels to be inpainted. Have to be the same size as 'img'.
         Unknown pixels has to be represented with 1, known pixels - with 0.
+    multichannel : boolean, optional
+        Defines if the last array dimension is a color dimension,
 
     Returns
     -------
@@ -43,117 +45,135 @@ def inpaint_biharmonic(img, mask):
         http://www.mathworks.com/matlabcentral/fileexchange/4551-inpaint-nans,
             method 3
     """
+    
+    def _inpaint(img, mask):
+        out = np.copy(img)
+        out_h, out_w = out.shape
+        out_l = out.size
 
-    if img.ndim != 2 or mask.ndim != 2:
-        raise ValueError('Only 2-dimensional arrays are supported')
-    if img.shape != mask.shape:
+        def _in_bounds(idx):
+            if len(idx) == 1:
+                return 0 <= idx <= out_l - 1
+            else:
+                return (0 <= idx[0] <= out_h - 1) and (0 <= idx[1] <= out_w - 1)
+        
+        # Find indexes of masked points in flatten array
+        mask_mn = np.array(np.where(mask)).T
+        mask_i = np.ravel_multi_index(np.where(mask), mask.shape)
+
+        # Initialize sparse matrix
+        # TODO: only points required for computation might be considered
+        matrix_unknown = sparse.lil_matrix((np.sum(mask), out.size), dtype=np.int32)
+        matrix_known = sparse.lil_matrix((np.sum(mask), out.size), dtype=np.int32)
+
+        # INFO: kernels can be reworked using scipy.signal.convolve2d
+        #       and np.array([0, 1, 0], [1, -4, 1], [0, 1, 0])
+
+        # 1 stage. Find points 2 or more pixels far from edges
+        kernel = [1, 2, -8, 2, 1, -8, 20, -8, 1, 2, -8, 2, 1]
+        offset = [-2 * out_w, -out_w - 1, -out_w, -out_w + 1,
+                  -2, -1, 0, 1, 2, out_w - 1, out_w, out_w + 1, 2 * out_w]
+
+        for idx, (i, (m, n)) in enumerate(zip(mask_i, mask_mn)):
+            if 2 <= m <= out_h - 3 and 2 <= n <= out_w - 3:
+                for k, o in zip(kernel, offset):
+                    if i + o in mask_i:
+                        matrix_unknown[idx, i + o] = k
+                    else:
+                        matrix_known[idx, i + o] = k
+
+        # 2 stage. Find points 1 pixel far from edges
+        kernel = [1, 1, -4, 1, 1]
+        offset = [-out_w, -1, 0, 1, out_w]
+
+        for idx, (i, (m, n)) in enumerate(zip(mask_i, mask_mn)):
+            if (m in [1, out_h - 2] and 1 <= n <= out_h - 2) or \
+               (n in [1, out_w - 2] and 1 <= m <= out_w - 2):
+                for k, o in zip(kernel, offset):
+                    if i + o in mask_i:
+                        matrix_unknown[idx, i + o] = k
+                    else:
+                        matrix_known[idx, i + o] = k
+
+        # 3 stage. Find points on the edges
+        kernel = [1, 1, -3, 1, 1]
+        offset = [-out_w, -1, 0, 1, out_w]
+        offset_mn = [(-1, 0), (0, -1), (0, 0), (0, 1), (1, 0)]
+
+        for idx, (i, (m, n)) in enumerate(zip(mask_i, mask_mn)):
+            if (m in [0, out_h - 1] and 1 <= n <= out_w - 2) or \
+               (n in [0, out_w - 1] and 1 <= m <= out_h - 2):
+                for k, o_mn in zip(kernel, offset_mn):
+                    if _in_bounds((m + o_mn[0], n + o_mn[1])):
+                        o = offset[offset_mn.index(o_mn)]
+                        if i + o in mask_i:
+                            matrix_unknown[idx, i + o] = k
+                        else:
+                            matrix_known[idx, i + o] = k
+
+        # 4 stage. Find corner points
+        kernel = [1, 1, -2, 1, 1]
+        offset = [-out_w, -1, 0, 1, out_w]
+        offset_mn = [(-1, 0), (0, -1), (0, 0), (0, 1), (1, 0)]
+
+        for idx, (i, (m, n)) in enumerate(zip(mask_i, mask_mn)):
+            if m in [0, out_h - 1] and n in [0, out_w - 1]:
+                for k, o_mn in zip(kernel, offset_mn):
+                    if _in_bounds((m + o_mn[0], n + o_mn[1])):
+                        o = offset[offset_mn.index(o_mn)]
+                        if i + o in mask_i:
+                            matrix_unknown[idx, i + o] = k
+                        else:
+                            matrix_known[idx, i + o] = k
+
+        # 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 columns
+        matrix_known = matrix_known.tocsr()
+        rhs = -(matrix_known * flat_diag_image).sum(axis=1)
+
+        # Solve linear system over defect points
+        matrix_unknown = matrix_unknown[:, mask_i]
+        matrix_unknown = sparse.csr_matrix(matrix_unknown)
+        result = spsolve(matrix_unknown, rhs)
+
+        # Handle enormous values
+        # TODO: consider images in [-1:1] scale
+        result[np.where(result < 0)] = 0
+        result[np.where(result > 1)] = 1
+
+        result = result.ravel()
+
+        # Put calculated points into the image
+        for idx, (m, n) in enumerate(mask_mn):
+            out[m, n] = result[idx]
+
+        return out
+
+    if img.ndim != (2 + multichannel) or mask.ndim != 2:
+        raise ValueError('Only 2D{+color} arrays are supported')
+
+    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)
-
-    out = np.copy(img)
-    out_h, out_w = out.shape
-    out_l = out.size
-
-    def _in_bounds(idx):
-        if len(idx) == 1:
-            return 0 <= idx <= out_l - 1
-        else:
-            return (0 <= idx[0] <= out_h - 1) and (0 <= idx[1] <= out_w - 1)
-
-    # Find indexes of masked points in flatten array
-    mask_mn = np.array(np.where(mask)).T
-    mask_i = np.ravel_multi_index(np.where(mask), mask.shape)
-
-    # Initialize sparse matrix
-    # TODO: only points required for computation might be considered
-    matrix_unknown = sparse.lil_matrix((np.sum(mask), out.size), dtype=np.int32)
-    matrix_known = sparse.lil_matrix((np.sum(mask), out.size), dtype=np.int32)
-
-    # INFO: kernels can be reworked using scipy.signal.convolve2d
-    #       and np.array([0, 1, 0], [1, -4, 1], [0, 1, 0])
-
-    # 1 stage. Find points 2 or more pixels far from edges
-    kernel = [1, 2, -8, 2, 1, -8, 20, -8, 1, 2, -8, 2, 1]
-    offset = [-2 * out_w, -out_w - 1, -out_w, -out_w + 1,
-              -2, -1, 0, 1, 2, out_w - 1, out_w, out_w + 1, 2 * out_w]
-
-    for idx, (i, (m, n)) in enumerate(zip(mask_i, mask_mn)):
-        if 2 <= m <= out_h - 3 and 2 <= n <= out_w - 3:
-            for k, o in zip(kernel, offset):
-                if i + o in mask_i:
-                    matrix_unknown[idx, i + o] = k
-                else:
-                    matrix_known[idx, i + o] = k
-
-    # 2 stage. Find points 1 pixel far from edges
-    kernel = [1, 1, -4, 1, 1]
-    offset = [-out_w, -1, 0, 1, out_w]
-
-    for idx, (i, (m, n)) in enumerate(zip(mask_i, mask_mn)):
-        if (m in [1, out_h - 2] and 1 <= n <= out_h - 2) or \
-           (n in [1, out_w - 2] and 1 <= m <= out_w - 2):
-            for k, o in zip(kernel, offset):
-                if i + o in mask_i:
-                    matrix_unknown[idx, i + o] = k
-                else:
-                    matrix_known[idx, i + o] = k
-
-    # 3 stage. Find points on the edges
-    kernel = [1, 1, -3, 1, 1]
-    offset = [-out_w, -1, 0, 1, out_w]
-    offset_mn = [(-1, 0), (0, -1), (0, 0), (0, 1), (1, 0)]
-
-    for idx, (i, (m, n)) in enumerate(zip(mask_i, mask_mn)):
-        if (m in [0, out_h - 1] and 1 <= n <= out_w - 2) or \
-           (n in [0, out_w - 1] and 1 <= m <= out_h - 2):
-            for k, o_mn in zip(kernel, offset_mn):
-                if _in_bounds((m + o_mn[0], n + o_mn[1])):
-                    o = offset[offset_mn.index(o_mn)]
-                    if i + o in mask_i:
-                        matrix_unknown[idx, i + o] = k
-                    else:
-                        matrix_known[idx, i + o] = k
-
-    # 4 stage. Find corner points
-    kernel = [1, 1, -2, 1, 1]
-    offset = [-out_w, -1, 0, 1, out_w]
-    offset_mn = [(-1, 0), (0, -1), (0, 0), (0, 1), (1, 0)]
-
-    for idx, (i, (m, n)) in enumerate(zip(mask_i, mask_mn)):
-        if m in [0, out_h - 1] and n in [0, out_w - 1]:
-            for k, o_mn in zip(kernel, offset_mn):
-                if _in_bounds((m + o_mn[0], n + o_mn[1])):
-                    o = offset[offset_mn.index(o_mn)]
-                    if i + o in mask_i:
-                        matrix_unknown[idx, i + o] = k
-                    else:
-                        matrix_known[idx, i + o] = k
-
-    # 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 columns
-    matrix_known = matrix_known.tocsr()
-    rhs = -(matrix_known * flat_diag_image).sum(axis=1)
-
-    # Solve linear system over defect points
-    matrix_unknown = matrix_unknown[:, mask_i]
-    matrix_unknown = sparse.csr_matrix(matrix_unknown)
-    result = spsolve(matrix_unknown, rhs)
-
-    # Handle enormous values
-    result[np.where(result < -1)] = -1
-    result[np.where(result > 1)] = 1
-
-    result = result.ravel()
     
-    # Put calculated points into the image
-    for idx, (m, n) in enumerate(mask_mn):
-        out[m, n] = result[idx]
+    if not multichannel:
+        img = img.reshape(img.shape + (1,))
+
+    out = np.zeros_like(img)
+    
+    for i in range(img.shape[-1]):
+        out[..., i] = _inpaint(img[..., i], mask)
+
+    if not multichannel:
+        out = out.reshape(out.shape[:-1])
 
     return out

From 0cc4402c195ac845a07a9308086fab74cded5b55 Mon Sep 17 00:00:00 2001
From: Egor Panfilov <multicolor.mood@gmail.com>
Date: Thu, 24 Dec 2015 02:27:18 +0300
Subject: [PATCH 08/13] DOC: Fixed docstring, moved to assert_nD

---
 skimage/restoration/inpaint.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/skimage/restoration/inpaint.py b/skimage/restoration/inpaint.py
index e41c68b9..b71ad324 100644
--- a/skimage/restoration/inpaint.py
+++ b/skimage/restoration/inpaint.py
@@ -2,6 +2,7 @@ from __future__ import print_function, division
 
 import numpy as np
 import skimage
+from .._shared.utils import assert_nD
 from scipy import sparse
 from scipy.sparse.linalg import spsolve
 
@@ -14,14 +15,15 @@ def inpaint_biharmonic(img, mask, multichannel=False):
     img : 2D{+color} np.array
         Input image.
     mask : 2D np.array
-        Array of pixels to be inpainted. Have to be the same size as 'img'.
-        Unknown pixels has to be represented with 1, known pixels - with 0.
+        Array of pixels to be inpainted. Have to be the same size as one
+        of the 'img' channels. Unknown pixels has to be represented with 1, 
+        known pixels - with 0.
     multichannel : boolean, optional
-        Defines if the last array dimension is a color dimension,
+        Defines if the last `img` dimension is a color dimension.
 
     Returns
     -------
-    out : 2-D np.array
+    out : 2D{+color} np.array
         Input image with masked pixels inpainted.
 
     Example
@@ -152,8 +154,8 @@ def inpaint_biharmonic(img, mask, multichannel=False):
 
         return out
 
-    if img.ndim != (2 + multichannel) or mask.ndim != 2:
-        raise ValueError('Only 2D{+color} arrays are supported')
+    assert_nD(img, 2 + multichannel)
+    assert_nD(mask, 2)
 
     img_baseshape = img.shape[:-1] if multichannel else img.shape
     if img_baseshape != mask.shape:

From f2f698b1901a761dc3c06ffe0f2476858fc13579 Mon Sep 17 00:00:00 2001
From: Egor Panfilov <multicolor.mood@gmail.com>
Date: Thu, 24 Dec 2015 22:32:56 +0300
Subject: [PATCH 09/13] DOC: Moved gallery example to category

---
 doc/examples/{ => filters}/plot_inpaint.py | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)
 rename doc/examples/{ => filters}/plot_inpaint.py (71%)

diff --git a/doc/examples/plot_inpaint.py b/doc/examples/filters/plot_inpaint.py
similarity index 71%
rename from doc/examples/plot_inpaint.py
rename to doc/examples/filters/plot_inpaint.py
index 7ccc73cc..ba863295 100644
--- a/doc/examples/plot_inpaint.py
+++ b/doc/examples/filters/plot_inpaint.py
@@ -39,18 +39,18 @@ image_defect[np.where(mask)] = 0
 
 image_result = inpaint.inpaint_biharmonic(image_defect, mask)
 
-fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, nrows=1, figsize=(10, 6))
+fig, axes = plt.subplots(ncols=3, nrows=1)
 
-ax1.set_title('Defected image')
-ax1.imshow(image_orig, cmap=plt.cm.gray, interpolation='nearest')
-ax1.set_xticks([]), ax1.set_yticks([])
+axes[0].set_title('Defected image')
+axes[0].imshow(image_orig, cmap=plt.cm.gray, interpolation='nearest')
+axes[0].set_xticks([]), axes[0].set_yticks([])
 
-ax2.set_title('Defect mask')
-ax2.imshow(mask, cmap=plt.cm.gray, interpolation='nearest')
-ax2.set_xticks([]), ax2.set_yticks([])
+axes[1].set_title('Defect mask')
+axes[1].imshow(mask, cmap=plt.cm.gray, interpolation='nearest')
+axes[1].set_xticks([]), axes[1].set_yticks([])
 
-ax3.set_title('Inpainted image')
-ax3.imshow(image_result, cmap=plt.cm.gray, interpolation='nearest')
-ax3.set_xticks([]), ax3.set_yticks([])
+axes[2].set_title('Inpainted image')
+axes[2].imshow(image_result, cmap=plt.cm.gray, interpolation='nearest')
+axes[2].set_xticks([]), axes[2].set_yticks([])
 
 plt.show()

From 89784631e79ff3358c267b3e807ef34ea4e9abb8 Mon Sep 17 00:00:00 2001
From: Egor Panfilov <multicolor.mood@gmail.com>
Date: Wed, 30 Dec 2015 22:45:27 +0300
Subject: [PATCH 10/13] ENH: Massively reworked code to enable nD-support

---
 doc/examples/filters/plot_inpaint.py      |  10 +-
 skimage/restoration/inpaint.py            | 127 +++++++---------------
 skimage/restoration/tests/test_inpaint.py |  26 +----
 3 files changed, 51 insertions(+), 112 deletions(-)

diff --git a/doc/examples/filters/plot_inpaint.py b/doc/examples/filters/plot_inpaint.py
index ba863295..ce661843 100644
--- a/doc/examples/filters/plot_inpaint.py
+++ b/doc/examples/filters/plot_inpaint.py
@@ -26,18 +26,20 @@ import matplotlib.pyplot as plt
 from skimage import data, color
 from skimage.restoration import inpaint
 
-image_orig = color.rgb2gray(data.astronaut())
+image_orig = data.astronaut()
 
 # Create mask with three defect regions: left, middle, right respectively
-mask = np.zeros_like(image_orig)
+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()
-image_defect[np.where(mask)] = 0
+for layer in range(image_defect.shape[-1]):
+    image_defect[np.where(mask)] = 0
 
-image_result = inpaint.inpaint_biharmonic(image_defect, mask)
+image_result = inpaint.inpaint_biharmonic(image_defect, mask, multichannel=True)
 
 fig, axes = plt.subplots(ncols=3, nrows=1)
 
diff --git a/skimage/restoration/inpaint.py b/skimage/restoration/inpaint.py
index b71ad324..f0c417b3 100644
--- a/skimage/restoration/inpaint.py
+++ b/skimage/restoration/inpaint.py
@@ -2,9 +2,9 @@ from __future__ import print_function, division
 
 import numpy as np
 import skimage
-from .._shared.utils import assert_nD
 from scipy import sparse
 from scipy.sparse.linalg import spsolve
+from scipy.ndimage.filters import laplace
 
 
 def inpaint_biharmonic(img, mask, multichannel=False):
@@ -12,18 +12,18 @@ def inpaint_biharmonic(img, mask, multichannel=False):
 
     Parameters
     ----------
-    img : 2D{+color} np.array
+    img : nD{+color channel} np.ndarray
         Input image.
-    mask : 2D np.array
+    mask : nD np.ndarray
         Array of pixels to be inpainted. Have to be the same size as one
         of the 'img' channels. Unknown pixels has to be represented with 1, 
         known pixels - with 0.
     multichannel : boolean, optional
-        Defines if the last `img` dimension is a color dimension.
+        If True, the last `img` dimension is considered as a color channel.
 
     Returns
     -------
-    out : 2D{+color} np.array
+    out : nD{+color channel} np.array
         Input image with masked pixels inpainted.
 
     Example
@@ -41,102 +41,58 @@ def inpaint_biharmonic(img, mask, multichannel=False):
     .. [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
-
-    Realization is based on:
-    .. [2]  John D'Errico,
-        http://www.mathworks.com/matlabcentral/fileexchange/4551-inpaint-nans,
-            method 3
     """
     
     def _inpaint(img, mask):
         out = np.copy(img)
-        out_h, out_w = out.shape
-        out_l = out.size
 
-        def _in_bounds(idx):
-            if len(idx) == 1:
-                return 0 <= idx <= out_l - 1
-            else:
-                return (0 <= idx[0] <= out_h - 1) and (0 <= idx[1] <= out_w - 1)
-        
+        # Initialize sparse matrices
+        matrix_unknown = sparse.lil_matrix((np.sum(mask), out.size))
+        matrix_known = sparse.lil_matrix((np.sum(mask), out.size))
+
+        def _get_neighborhood(idx, radii):
+            bounds_lo = (idx - radii).clip(min=0)
+            bounds_hi = (idx + np.add(radii, 1)).clip(max=out.shape)
+            return bounds_lo, bounds_hi
+
         # Find indexes of masked points in flatten array
-        mask_mn = np.array(np.where(mask)).T
         mask_i = np.ravel_multi_index(np.where(mask), mask.shape)
 
-        # Initialize sparse matrix
-        # TODO: only points required for computation might be considered
-        matrix_unknown = sparse.lil_matrix((np.sum(mask), out.size), dtype=np.int32)
-        matrix_known = sparse.lil_matrix((np.sum(mask), out.size), dtype=np.int32)
+        # Find masked points and prepare them to be easily enumerate over
+        mask_pts = np.array(np.where(mask)).T
 
-        # INFO: kernels can be reworked using scipy.signal.convolve2d
-        #       and np.array([0, 1, 0], [1, -4, 1], [0, 1, 0])
+        # Iterate over masked points
+        for mask_pt_n, mask_pt_idx in enumerate(mask_pts):
+            # Get bounded neighborhood of selected radii
+            b_lo, b_hi = _get_neighborhood(mask_pt_idx, radii=np.array([2]))
 
-        # 1 stage. Find points 2 or more pixels far from edges
-        kernel = [1, 2, -8, 2, 1, -8, 20, -8, 1, 2, -8, 2, 1]
-        offset = [-2 * out_w, -out_w - 1, -out_w, -out_w + 1,
-                  -2, -1, 0, 1, 2, out_w - 1, out_w, out_w + 1, 2 * out_w]
+            # 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))
 
-        for idx, (i, (m, n)) in enumerate(zip(mask_i, mask_mn)):
-            if 2 <= m <= out_h - 3 and 2 <= n <= out_w - 3:
-                for k, o in zip(kernel, offset):
-                    if i + o in mask_i:
-                        matrix_unknown[idx, i + o] = k
-                    else:
-                        matrix_known[idx, i + o] = k
+            # 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)
 
-        # 2 stage. Find points 1 pixel far from edges
-        kernel = [1, 1, -4, 1, 1]
-        offset = [-out_w, -1, 0, 1, out_w]
-
-        for idx, (i, (m, n)) in enumerate(zip(mask_i, mask_mn)):
-            if (m in [1, out_h - 2] and 1 <= n <= out_h - 2) or \
-               (n in [1, out_w - 2] and 1 <= m <= out_w - 2):
-                for k, o in zip(kernel, offset):
-                    if i + o in mask_i:
-                        matrix_unknown[idx, i + o] = k
-                    else:
-                        matrix_known[idx, i + o] = k
-
-        # 3 stage. Find points on the edges
-        kernel = [1, 1, -3, 1, 1]
-        offset = [-out_w, -1, 0, 1, out_w]
-        offset_mn = [(-1, 0), (0, -1), (0, 0), (0, 1), (1, 0)]
-
-        for idx, (i, (m, n)) in enumerate(zip(mask_i, mask_mn)):
-            if (m in [0, out_h - 1] and 1 <= n <= out_w - 2) or \
-               (n in [0, out_w - 1] and 1 <= m <= out_h - 2):
-                for k, o_mn in zip(kernel, offset_mn):
-                    if _in_bounds((m + o_mn[0], n + o_mn[1])):
-                        o = offset[offset_mn.index(o_mn)]
-                        if i + o in mask_i:
-                            matrix_unknown[idx, i + o] = k
-                        else:
-                            matrix_known[idx, i + o] = k
-
-        # 4 stage. Find corner points
-        kernel = [1, 1, -2, 1, 1]
-        offset = [-out_w, -1, 0, 1, out_w]
-        offset_mn = [(-1, 0), (0, -1), (0, 0), (0, 1), (1, 0)]
-
-        for idx, (i, (m, n)) in enumerate(zip(mask_i, mask_mn)):
-            if m in [0, out_h - 1] and n in [0, out_w - 1]:
-                for k, o_mn in zip(kernel, offset_mn):
-                    if _in_bounds((m + o_mn[0], n + o_mn[1])):
-                        o = offset[offset_mn.index(o_mn)]
-                        if i + o in mask_i:
-                            matrix_unknown[idx, i + o] = k
-                        else:
-                            matrix_known[idx, i + o] = k
+                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 columns
+        # 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 over defect points
+        # Solve linear system for masked points
         matrix_unknown = matrix_unknown[:, mask_i]
         matrix_unknown = sparse.csr_matrix(matrix_unknown)
         result = spsolve(matrix_unknown, rhs)
@@ -148,15 +104,12 @@ def inpaint_biharmonic(img, mask, multichannel=False):
 
         result = result.ravel()
 
-        # Put calculated points into the image
-        for idx, (m, n) in enumerate(mask_mn):
-            out[m, n] = result[idx]
+        # 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
 
-    assert_nD(img, 2 + multichannel)
-    assert_nD(mask, 2)
-
     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')
diff --git a/skimage/restoration/tests/test_inpaint.py b/skimage/restoration/tests/test_inpaint.py
index 1fe3e766..abe89df3 100644
--- a/skimage/restoration/tests/test_inpaint.py
+++ b/skimage/restoration/tests/test_inpaint.py
@@ -14,27 +14,11 @@ def test_inpaint_biharmonic():
     mask[0, 4:] = 1
     out = inpaint.inpaint_biharmonic(img, mask)
     ref = np.array(
-        [[0., 0.0625, 0.25, 0.5625, 0.56947314],
-         [0., 0.0625, 0.25, 0.47029959, 0.57644628],
-         [0., 0.0625, 0.24664256, 0.49225207, 0.68956612],
-         [0., 0.0625, 0.25, 0.5625, 1.],
-         [0., 0.0625, 0.25, 0.5625, 1.]]
-    )
-    assert_allclose(ref, out)
-
-
-def test_inpaint_edges():
-    img = np.tile(np.square(np.linspace(0, 1, 5)), (5, 1))
-    mask = np.zeros_like(img)
-    mask[[0, -1], :] = 1
-    mask[:, [0, -1]] = 1
-    out = inpaint.inpaint_biharmonic(img, mask)
-    ref = np.array(
-        [[0.12199519, 0.15599245, 0.28348214, 0.44445398, 0.48737981],
-         [0.08799794, 0.0625, 0.25, 0.5625, 0.53030563],
-         [0.07949863, 0.0625, 0.25, 0.5625, 0.54103709],
-         [0.08799794, 0.0625, 0.25, 0.5625, 0.53030563],
-         [0.12199519, 0.15599245, 0.28348214, 0.44445398, 0.48737981]]
+        [[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)
 

From 893d93749116aa7605191b7ae03bd4cf11a4e323 Mon Sep 17 00:00:00 2001
From: Egor Panfilov <multicolor.mood@gmail.com>
Date: Fri, 15 Jan 2016 13:59:32 +0300
Subject: [PATCH 11/13] ENH: Moved nested functions, refactor code

---
 skimage/restoration/inpaint.py | 164 +++++++++++++++++----------------
 1 file changed, 83 insertions(+), 81 deletions(-)

diff --git a/skimage/restoration/inpaint.py b/skimage/restoration/inpaint.py
index f0c417b3..1b318b97 100644
--- a/skimage/restoration/inpaint.py
+++ b/skimage/restoration/inpaint.py
@@ -1,4 +1,4 @@
-from __future__ import print_function, division
+from __future__ import division
 
 import numpy as np
 import skimage
@@ -7,23 +7,89 @@ 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 : nD{+color channel} np.ndarray
+    img : (M, N[, ..., P][, C]) ndarray
         Input image.
-    mask : nD np.ndarray
-        Array of pixels to be inpainted. Have to be the same size as one
-        of the 'img' channels. Unknown pixels has to be represented with 1, 
+    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.
+        If True, the last `img` dimension is considered as a color channel,
+        otherwise as spatial.
 
     Returns
     -------
-    out : nD{+color channel} np.array
+    out : (M, N[, ..., P][, C] ndarray
         Input image with masked pixels inpainted.
 
     Example
@@ -43,73 +109,6 @@ def inpaint_biharmonic(img, mask, multichannel=False):
             http://www.ima.umn.edu/~damelin/biharmonic
     """
     
-    def _inpaint(img, mask):
-        out = np.copy(img)
-
-        # Initialize sparse matrices
-        matrix_unknown = sparse.lil_matrix((np.sum(mask), out.size))
-        matrix_known = sparse.lil_matrix((np.sum(mask), out.size))
-
-        def _get_neighborhood(idx, radii):
-            bounds_lo = (idx - radii).clip(min=0)
-            bounds_hi = (idx + np.add(radii, 1)).clip(max=out.shape)
-            return bounds_lo, bounds_hi
-
-        # 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 radii
-            b_lo, b_hi = _get_neighborhood(mask_pt_idx, radii=np.array([2]))
-
-            # 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
-        # TODO: consider images in [-1:1] scale
-        result[np.where(result < 0)] = 0
-        result[np.where(result > 1)] = 1
-
-        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
-
     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')
@@ -119,16 +118,19 @@ def inpaint_biharmonic(img, mask, multichannel=False):
 
     img = skimage.img_as_float(img)
     mask = mask.astype(np.bool)
-    
-    if not multichannel:
-        img = img.reshape(img.shape + (1,))
 
-    out = np.zeros_like(img)
-    
+    if not multichannel:
+        img = img[..., np.newaxis]
+
+    out = np.copy(img)
+
     for i in range(img.shape[-1]):
-        out[..., i] = _inpaint(img[..., i], mask)
+        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.reshape(out.shape[:-1])
+        out = out[..., 0]
 
     return out

From 98fec5211ae7b4ae2d17ce04bd008693e1e29f46 Mon Sep 17 00:00:00 2001
From: Egor Panfilov <multicolor.mood@gmail.com>
Date: Fri, 15 Jan 2016 20:04:08 +0300
Subject: [PATCH 12/13] FIX: Returned handle for empty input arrays

---
 skimage/restoration/inpaint.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/skimage/restoration/inpaint.py b/skimage/restoration/inpaint.py
index 1b318b97..bf696e04 100644
--- a/skimage/restoration/inpaint.py
+++ b/skimage/restoration/inpaint.py
@@ -77,9 +77,9 @@ def inpaint_biharmonic(img, mask, multichannel=False):
 
     Parameters
     ----------
-    img : (M, N[, ..., P][, C]) ndarray
+    img : (M[, N[, ..., P]][, C]) ndarray
         Input image.
-    mask : (M, N[, ..., P]) ndarray
+    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.
@@ -89,7 +89,7 @@ def inpaint_biharmonic(img, mask, multichannel=False):
 
     Returns
     -------
-    out : (M, N[, ..., P][, C] ndarray
+    out : (M[, N[, ..., P]][, C]) ndarray
         Input image with masked pixels inpainted.
 
     Example
@@ -108,6 +108,9 @@ def inpaint_biharmonic(img, mask, multichannel=False):
             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:

From 3ae2f1236bc149a972a745c6f6955c2864577011 Mon Sep 17 00:00:00 2001
From: Egor Panfilov <multicolor.mood@gmail.com>
Date: Sat, 16 Jan 2016 20:02:24 +0300
Subject: [PATCH 13/13] TST, DOC: Added test for 3D spatial data, cleaned up
 gallery example

---
 doc/examples/filters/plot_inpaint.py      |  6 ++---
 skimage/restoration/tests/test_inpaint.py | 29 ++++++++++++++++++++++-
 2 files changed, 31 insertions(+), 4 deletions(-)

diff --git a/doc/examples/filters/plot_inpaint.py b/doc/examples/filters/plot_inpaint.py
index ce661843..104694ec 100644
--- a/doc/examples/filters/plot_inpaint.py
+++ b/doc/examples/filters/plot_inpaint.py
@@ -44,15 +44,15 @@ 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, cmap=plt.cm.gray, interpolation='nearest')
+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, interpolation='nearest')
+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, cmap=plt.cm.gray, interpolation='nearest')
+axes[2].imshow(image_result)
 axes[2].set_xticks([]), axes[2].set_yticks([])
 
 plt.show()
diff --git a/skimage/restoration/tests/test_inpaint.py b/skimage/restoration/tests/test_inpaint.py
index abe89df3..04e8c8c2 100644
--- a/skimage/restoration/tests/test_inpaint.py
+++ b/skimage/restoration/tests/test_inpaint.py
@@ -6,12 +6,13 @@ from numpy.testing import (run_module_suite, assert_allclose,
 from skimage.restoration import inpaint
 
 
-def test_inpaint_biharmonic():
+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],
@@ -23,6 +24,32 @@ def test_inpaint_biharmonic():
     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)