A study of local smoothness-informed convolutional neural network models for image inpainting

CNN has been designed since the 1980s and was later diversely developed, for example, the well-known CNN developed by Fukushima [13] and Zhang et al. [14]. We utilize CNN to the gray-level image inpainting and design the CNN with the same sized output to the input for every layer to avoid the data loss. The CNN accepts one damaged image and outputs one inpainted image; the image is convoluted by the kernels and activated from input to output. The kernels stride over the image by one step to complete the convolution for the whole image and the Rectified Linear Units [15] is used as the activation operator. The main computational details and architecture of the CNN are as follow.

Consider a CNN (denoted as g(x, y; w, b)) with L layers consist of kernels with odd size, the gray value at location (x, y) of the i^th output of the (l + 1)^th layer reads (2) gil+1(x,y;wil+1,bil+1)=σ(∑k=1Kl∑p=−ss∑q=−ssgkl(x+p,y+q;wil,bil)wikl+1(c+p,c+q)+bil+1) s=S−12,c=S+12, \matrix{ {g_i^{l + 1}(x,y;w_i^{l + 1},b_i^{l + 1}) = \sigma \left( {\sum\limits_{k = 1}^{{K^l}} \sum\limits_{p = - s}^s \sum\limits_{q = - s}^s g_k^l(x + p,y + q;w_i^l,b_i^l)w_{ik}^{l + 1}(c + p,c + q) + b_i^{l + 1}} \right)} \hfill \cr {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;s = {{S - 1} \over 2},c = {{S + 1} \over 2},} \hfill \cr } where σ is an activation operator, K^l is the amount of the outputs of the l^th layer, wikl+1 w_{ik}^{l + 1} denotes the k^th kernel corresponding to the i^th output of the (l + 1)^th layer, bil+1 b_i^{l + 1} denotes the i^th component of the bias vector of the (l + 1)^th layer, g^l is the input of the l + 1 layer, l = 0, ⋯, L − 1, S is the kernel size, s is the convolutional stride depart from the center of the kernel, c is the center of the kernel.

Training the CNN aims to learn the kernels and biases using the loss function presented in Section 3.2.

Let g(x, y) as the unknown function representing the inpainted image. Inpainting aims to determine g(x, y). The main idea of the deep learning methods for image inpainting is to generate an estimation with the available data of Ω₂ and train the neural network (finding the kernels and biases) via the loss minimization process. The basic loss function is the least square total loss (3) Ld=∫Ωh(g−u)2dxdy. {L_d} = \int_\Omega h{(g - u)^2}dxdy. However, this loss function lacks the connectivity conservation between Ω₁ and its surroundings. Without connectivity, the estimation in Ω₁ has weaker smoothness, which usually results in odd or irregular estimations. TV is wildly used for edge maintaining or low-order smoothness, which was earlier designed for noise removal by Rudin et al. [4]. The TV-based energy functional for inpainting (see, e.g., [2, 16]) with f = hu is given by (4) Lcs=∫Ω(h2(g−u)2+λ|∇g|)dxdy, {L_{cs}} = \int_\Omega \left( {{h \over 2}{{(g - u)}^2} + \lambda |\nabla g|} \right)dxdy, where λ > 0 is a weighted coefficient, ∇ is the gradient operator. We may choose a very small λ if noise removal is not a major task.

The variational problem to find a solution g is to minimize the energy functional (4). Ultimately, determining g is to solve the associated Euler-Lagrange equation (5) h(g−u)−λ∇⋅(∇g|∇g|2+ε)=0 h(g - u) - \lambda \nabla \cdot \left( {{{\nabla g} \over {\sqrt {|\nabla g{|^2} + \varepsilon } }}} \right) = 0 with ɛ = 0. Here ɛ > 0 is a small parameter for numerical implementation to avoid dividing by zero. We set ɛ according to the flatness of the detectable domain since the estimation depends on the given information, and ɛ is given by (6) ε=e−1A∫Ω2|∇g|2n2dxdy, \varepsilon = {e^{ - {1 \over A}\int_{{\Omega _2}} {{|\nabla g{|^2}} \over {{n^2}}}dxdy}}, where A is the total area of Ω, n > 0 is used to adjust the gradient magnitude. The given parameter ɛ ∈ (0, 1] depends inversely on the flatness to conserve the TV role which follows the idea of edge detection [17]. For comparison in implementation, we actually minimize the shifted TV of (4), which is |∇g|2+ε. \sqrt {|\nabla g{|^2} + \varepsilon } . We note that the parameter ɛ is to avoid dividing by zero in minimizing (4). Traditionally, the Euler-Lagrange equation is solved in a numerical iterative algorithm. Let (7) ∂g(x,y;t)∂t=h(g(x,y;t)−u)−λ∇⋅(∇g(x,y;t)|∇g(x,y;t)|2+ε) {{\partial g(x,y;t)} \over {\partial t}} = h(g(x,y;t) - u) - \lambda \nabla \cdot \left( {{{\nabla g(x,y;t)} \over {\sqrt {|\nabla g(x,y;t{{)|}^2} + \varepsilon } }}} \right) denotes the spatial variant on temporal level, t is an artificially introduced temporal variable for numerical realization. The final solution g satisfies ∂g(x,y;t)∂t→0 {{\partial g(x,y;t)} \over {\partial t}} \to 0 . However, for the convenience of using the optimization algorithm we mainly apply to train the CNN in this paper (see next subsection), we convert the Euler-Lagrange equation to the minimization problem of the following loss function (8) Leu=∫Ω|h(g−u)−λ∇⋅(∇g|∇g|2+ε)|2dxdy. {L_{eu}} = \int_\Omega {\left| {h(g - u) - \lambda \nabla \cdot \left( {{{\nabla g} \over {\sqrt {|\nabla g{|^2} + \varepsilon } }}} \right)} \right|^2}dxdy. Consequently, solving the variational problem for the energy functional (4) shifts to the minimization problem of the loss function (8). Note that there might be slight differences between their solutions, due to the addition of small ɛ.

The TV on the whole domain Ω was initially for noise removal. For image inpainting without noise removal to the partial unknown domain, it mainly to enhance the connectivity between the corrupted domain and the detectable domain. In this paper, we explore a local H¹ or H² regularization for the connectivity, which features one of the necessary geometrical details of the images. The regularization is to condition the low-order derivatives of the sub-domain of detectable domain Ω₂ that mostly connected to Ω₁ on the given function f (inspired by [10]) to preserve the various levels of smoothness to enhance connection between Ω₁ and its surroundings. With this point, g subjects to (9) Dg(x,y)=Du(x,y), (x,y)∈Ω3, Dg(x,y) = Du(x,y),\quad (x,y) \in {\Omega _3}, where D is a first or second order differential operator, Ω₃ ⊂ Ω₂ is around Ω₁. Assume f ∈ ℂ²(Ω₃). Theoretically, the first order differential operator ∇ (which we call local H¹ regularization) is able to enhance the smoothness of g approaching to f, and the second order differential operator Δ (which we call local H² regularization) conducts a smoother connection.

Denote V_g(x, y) = Dg(x, y), V_u(x, y) = Du(x, y). Combining the loss function L_d, determining g is to minimize the following loss (10) LD=∫Ωh(g−u)2dxdy+β∫Ω3|Vg−Vu|2dxdy, {L_D} = \int_\Omega h{(g - u)^2}dxdy + \beta \int_{{\Omega _3}} |{V_g} - {V_u}{|^2}dxdy, β is a penalty coefficient.

In this section, we present the computational details of deep learning models mentioned in previous section for image inpainting. The integral is discretely represented in a summation and derivatives are approximated by finite differences. Avoiding data loss, we pad zero sides to the boundaries of the image matrix.

Denote f = f (x, y) as the input image, w as the weighted matrixes representing the kernels, b as the bias vectors, g = g(x, y; w, b) as the estimation of the CNN, in this subsection (x, y) ∈ [0, W) × [0, W) ⊂ ℕ × ℕ as the coordinate of any pixel or discrete point of the image and W as the length of the boundary of the image.

Computing the gradient for the TV, we tend to use the Sobel operator [18] {S_x, S_y} to approximate the first order partial derivatives g_x and g_y, where (11) Sx=[−101−202−101], Sy=[121000−1−2−1]. {S_x} = \left[ {\matrix{ { - 1} & 0 & 1 \cr { - 2} & 0 & 2 \cr { - 1} & 0 & 1 \cr } } \right],\quad {S_y} = \left[ {\matrix{ 1 & 2 & 1 \cr 0 & 0 & 0 \cr { - 1} & { - 2} & { - 1} \cr } } \right]. With this operator, ∀(x, y) ∈ [0, W) × [0, W), (12) gx(x,y)≈(Sx*g)(x,y),gy(x,y)≈(Sy*g)(x,y), {g_x}(x,y) \approx ({S_x}*g)(x,y),{g_y}(x,y) \approx ({S_y}*g)(x,y), where * is the convolution operator symbol.

Computing the second order partial derivatives, we consider the convolutions twice by the operators S_x and S_y, so the divergence is approximated as (13) ∇⋅(∇g|∇g|2+ε)(x,y)≈Sx*(gx|∇g|2+ε)(x,y)+Sy*(gy|∇g|2+ε)(x,y) \nabla \cdot \left( {{{\nabla g} \over {\sqrt {|\nabla g{|^2} + \varepsilon } }}} \right)(x,y) \approx {S_x}*\left( {{{{g_x}} \over {\sqrt {|\nabla g{|^2} + \varepsilon } }}} \right)(x,y) + {S_y}*\left( {{{{g_y}} \over {\sqrt {|\nabla g{|^2} + \varepsilon } }}} \right)(x,y)

A H¹ or H² regularization in Ω₃ is respectively denoted as (14) |Vg−Vu|2=|∇g−∇u|2 |{V_g} - {V_u}{|^2} = |\nabla g - \nabla u{|^2} or (15) |Vg−Vu|2=|Δg−Δu|2. |{V_g} - {V_u}{|^2} = |\Delta g - \Delta u{|^2}.

Near the boundary of Ω₁, it is not convenient to compute derivatives with the operators S_x, S_y, so we approximate the derivatives ∀(x, y) ∈ Ω₃ as follow (16) gx(x,y)≈(D1*g)(x,y),gy(x,y)≈(D1T*g)(x,y)ux(x,y)≈(D1*u)(x,y),uy(x,y)≈(D1T*u)(x,y), \matrix{ {{g_x}(x,y) \approx ({D_1}*g)(x,y),{g_y}(x,y) \approx (D_1^T*g)(x,y)} \hfill \cr {{u_x}(x,y) \approx ({D_1}*u)(x,y),{u_y}(x,y) \approx (D_1^T*u)(x,y),} \hfill \cr } where (17) D1=[−1,1]. {D_1} = [ - 1,1]. For the rectangular damaged domain considered in this paper, the second order partial derivatives are approximated by second-order central difference operator, or applying D₁ or D1T D_1^T twice.

Let d = [0, W) × [0, W), d₂ ⊂ d as the discrete Ω₂, d₃ ⊂ d₂ as the discrete Ω₃. With the instructive knowledge presented in the above subsection, determining g is to solve the optimal problems minw,bL˜d,minw,bL˜eu,minw,bL˜cs,minw,bL˜D \mathop {\min }\limits_{w,b} {\tilde L_d},\mathop {\min }\limits_{w,b} {\tilde L_{eu}},\mathop {\min }\limits_{w,b} {\tilde L_{cs}},\mathop {\min }\limits_{w,b} {\tilde L_D} separately, where (18) L˜d=1W2∑(x,y)∈dh(g−u)2L˜eu=1W2∑(x,y)∈d[h(g−u)−λ∇⋅(∇g|∇g|2+ε)]2L˜cs=1W2∑(x,y)∈d[h2(g−u)2+λgx2+gy2+ε]L˜D=1W2∑(x,y)∈dh(g−u)2+βW2∑(x,y)∈d3|Vg−Vu|2. \matrix{ {{{\tilde L}_d} = {1 \over {{W^2}}}\sum\limits_{(x,y) \in d} h{{(g - u)}^2}} \hfill \cr {{{\tilde L}_{eu}} = {1 \over {{W^2}}}\sum\limits_{(x,y) \in d} {{\left[ {h(g - u) - \lambda \nabla \cdot \left( {{{\nabla g} \over {\sqrt {|\nabla g{|^2} + \varepsilon } }}} \right)} \right]}^2}} \hfill \cr {{{\tilde L}_{cs}} = {1 \over {{W^2}}}\sum\limits_{(x,y) \in d} \left[ {{h \over 2}{{(g - u)}^2} + \lambda \sqrt {g_x^2 + g_y^2 + \varepsilon } } \right]} \hfill \cr {{{\tilde L}_D} = {1 \over {{W^2}}}\sum\limits_{(x,y) \in d} h{{(g - u)}^2} + {\beta \over {{W^2}}}\sum\limits_{(x,y) \in {d_3}} |{V_g} - {V_u}{|^2}.} \hfill \cr }

We adopt the Adam [19] algorithm to update the parameters. Computational process is stated with loss function L˜D {\tilde L_D} in Algorithm 1.