**1. Introduction**

Image completion aims to synthesize missing regions in an incomplete image or a video sequence with the content-aware information captured from accessible or unperturbed regions. The missing regions may result from removing unwanted sub-areas, occlusion, or unobserved or considerably damaged random pixels.

Image completion is one of the fundamental research topics in the area of computer vision and graphics technology, motivated by widespread applications in various applied sciences [1]. It is used for the restoration of old photographs, paintings, and films by removing scratches, dust spots, occlusions, or other user-marked objects, such as annotations, subtitles, stamps, logos, etc. In telecommunications, image completion techniques may fix error concealment problems or recover missing data-blocks lost during transmission or video compression [2,3]. The recent works also emphasize their usefulness in remote sensing imaging to remove occlusions, such as clouds and "dead" pixels [4,5].

The topic of image completion has been extensively studied for at least two decades, resulting in the development of several computational approaches to address this problem. The pioneering work in fully automated digital image inpainting dates back to 2000 when Bertalmio et al. [6] introduced a milestone algorithm that was able to automatically fill in missing regions from their neighboring information without providing user-defined sample images. This algorithm is based on partial differential equations (PDEs) that determine isophotes (brightness-level lines) along which the information on the surrounding structure is propagated upwards. It is a fully automatic algorithm, but it inpaints efficiently only in narrow lines, does not recover textures, and produces blurred results. Another approach to image inpainting is to synthesize the texture. This concept was introduced

by Efros and Leung [7] in 1999, but their algorithm requires a sample texture image to recover the texture in the missing region. Bertalmio et al. [8] also developed a hybrid version of both inpainting techniques, in which an incomplete image is decomposed into its texture and structure components, one reconstructed by texture synthesis, and another by PDE-based image inpainting. When a region to be completed is relatively large, the exemplar-based image synthesis algorithms [9–11] or their hybrid versions [12] seem to be more appropriate. The hybrid strategies have been studied in many other research papers [13–16], and currently, image completion based on simultaneous fill-in of texture and structure is a fundamental approach, especially for recovering large missing regions. Various neural network architectures, e.g., convolutional neural network and generative adversarial network, can also perform hybrid image completion [17–20]. A survey of image completion methods can be found in [1,5,21,22].

The above-mentioned image completion methods are efficient even for very large missing regions; however, they cannot be applied for incomplete images with many uniformly distributed missing pixels (e.g., about 90 %). No texture information can be learned from any subregions of such an image. There is also no continuous boundary of missing regions, and hence the neighboring information cannot be propagated towards the centers of missing regions. In such cases, different methods must be applied. Assuming that an incomplete image has a low-rank structure and satisfies the incoherence condition [23], usually represented by clusters of similar patches, the image completion boils down to a rank minimization problem. Since it is an NP-hard problem, many computational strategies have been proposed to approximate it by a convex relaxation, usually involving matrix or tensor decomposition methods. One of them assumes a convex approximation of a rank minimization problem with a nuclear-norm minimization problem, which can be solved easily using singular value decomposition (SVD) [24–26]. Due to the orthogonality condition, low-rank representations yielded by SVD contain negative values, which is not profitable for representing a large class of images. Other low-rank models (not necessarily restricted to SVD) can be used to tackle this problem [27].

One of the commonly used methods for extracting low-rank part-based representations from nonnegative matrices is nonnegative matrix factorization (NMF) [28]. It has already found many relevant applications in image analysis, and can also be used for solving image completion problems [29,30]. In this approach, the missing regions are sequentially updated with an NMF-based low-rank approximation of an observed image, which resembles the phenomena of propagating the neighboring information towards the missing regions in the PDE-based inpainting methods. However, due to the non-convexity of alternating optimization, NMF is sensitive to its initialization, especially for insufficiently sparse data. The ambiguity effects can be considerably relaxed if tensor decomposition models are used [31]. Moreover, multi-linear decomposition models prevent cross-modal interactions, which is particularly useful for image representations. Such a low-rank image completion methodology is mostly based on the concept of tensor completion issues that have been extensively studied [32–36]. There are many tensor decomposition models that are used for image completion tasks, including the fundamental ones, such as CANDECOMP/PARAFAC (CP) [37,38] and Tucker decomposition [39–41], as well as tensor networks, such as tensor ring [42], tensor train [43,44], hierarchical Tucker decomposition [45], and other tensor decomposition models [46].

Tensor completion methods are intrinsically addressed for processing color images (3D multi-way arrays), but they can also be applied to gray-scale images using various tensorization or scanning operators, e.g., the ket augmentation [47]. They are also very flexible in incorporating penalty terms or constraints; however, their efficiency strongly depends on the image to be completed. In practice, a low-rank representation is always a certain approximation of the underlying image, controlled by the rank of a tensor decomposition model. The problem of rank selection is regarded in terms of a trade-off between under- and over-fitting, and many approaches exist to tackle it. For example, Yokota et al. [38] proposed to increase the rank with recursive updates. In the early recursive steps, a low-rank structure is recovered, and then it is gradually updated to a higher-rank structure that contains more details. The rank can also be controlled by the decreasing rank

procedure [46] or using the singular value thresholding strategy [43]. Thus, one of the drawbacks of this kind of image completion method is the problem of selecting the optimal rank of decomposition, and it is usually resolved by heuristic procedures. Moreover, the tensor decomposition-based image-completion methods usually involve a high computational cost if all factor matrices or core tensors are updated in each iterative step.

To relax the problem with the right rank selection while keeping computational costs at a very low level, we propose a very simple alternative approach to low-rank image completion. Our strategy is based on the Tucker decomposition model in which the full-rank factor matrices are previously estimated with a hybrid connection of two interpolation methods. Since the factor matrices are precomputed, only the core tensor must be estimated. Despite the full-rank assumption, it involves a relatively low computational cost because the core tensor is sparse with non-zero entries determined by the available pixels. Motivated by several works [48–51] on the use of various interpolation methods for solving image completion problems, other recent works [52–56] on image processing aspects, and the concept of tensor product spline surfaces [57], we show the relationship of the Tucker decomposition with factorizable radial basis function (RBF) interpolation and use it to compute the factor matrices. RBF interpolation is a mesh-free method, which is very profitable for recovering irregularly distributed missing pixels, but it may incorrectly approximate linear structure. Hence, we combine it with low-degree polynomial interpolation, and both interpolation methods can be expressed by the Tucker decomposition model. Adopting the idea of PDE-based inpainting, we propose to compute the interpolants using only restricted surrounding subareas, instead of all accessible pixels. Hence, the whole image is divided into overlapping blocks, and the Tucker decomposition model is applied to each block separately. As many interpolation methods do not approximate the boundary entries well, the overlapping is necessary to avoid discontinuity effects. The proposed methodology is applied to various image completion problems, where the incomplete images are obtained from true images by removing many random pixels or many small holes. One of the experiments is performed for resolution up-scaling, where a low-resolution image is up-scaled to high resolution using the proposed algorithm. All the experiments demonstrate that the proposed method considerably outperforms the baseline and state-of-the-art low-rank image-completion methods in terms of the performance, and it is much faster than the methods based on tensor decompositions.

The remainder of this paper is organized as follows. Section 2 reviews some preliminary knowledge about fundamental tensor operations, the Tucker decomposition model, low-rank tensor completion, and RBF interpolation. The proposed algorithm is described in Section 3. The numerical experiments performed on various image completion problems are given and discussed in Section 4. The last section contains the conclusions.

#### **2. Preliminary**

Mathematical notations and preliminaries of tensor decomposition models are adopted from [31]. Tensors are multi-way arrays, denoted by calligraphic letters (e.g., <sup>Y</sup>). Let Y ∈ <sup>R</sup>*I*1×...×*IN* be the *N*-way array. The elements of Y are denoted as *yi*1,...,*iN* , where ∀*n* : 1 ≤ *in* ≤ *In*, *n* = 1, ... , *N*. Boldface uppercase letters (e.g., *Y* ) denote matrices; lowercase boldface ones stand for vectors (e.g., *y*); non-bold letters are scalars. The vector *y<sup>j</sup>* contains the *j*-th column of *Y*. The symbol || · ||*<sup>F</sup>* denotes the Frobenius norm of a matrix; || · || denotes the 2-nd norm. The sets of real numbers and natural numbers are represented by <sup>R</sup> and <sup>N</sup> , respectively. The symbols *<sup>x</sup>* and *<sup>x</sup>* stand for the floor and ceiling functions of *x*, respectively.

Let <sup>M</sup> = [*mi*1,...,*iN* ] <sup>∈</sup> <sup>R</sup>*I*1×...×*IN* be the *<sup>N</sup>*-th way observed tensor with missing entries, and <sup>Ω</sup> = [*ωi*1,...,*iN* ] <sup>∈</sup> <sup>R</sup>*I*1×...×*IN* be a binary tensor that indicates the locations of available entries in M. If *mi*1,...,*iN* is observed, then *ωi*1,...,*iN* = 1; otherwise, *ωi*1,...,*iN* = 0. The locations of positive entries in <sup>Ω</sup>¯ <sup>=</sup> <sup>1</sup> <sup>−</sup> <sup>Ω</sup> correspond to missing entries in <sup>M</sup>. The number of observed entries is equal to |Ω| = {-(*i*1,..., *iN*) : *ωi*1,...,*iN* > 0}.

*Appl. Sci.* **2020**, *10*, 797

**Definition 1.** *Let*

$$\mathcal{J} = [\mathfrak{g}\_{i\_1,\ldots,i\_N}], \quad \text{where} \quad \mathfrak{H}\_{i\_1,\ldots,i\_N} = \begin{cases} m\_{i\_1,\ldots,i\_N} & \text{if } & \omega\_{i\_1,\ldots,i\_N} = 1 \\ 0 & \text{otherwise} \end{cases} \tag{1}$$

*be a zero-filled incomplete tensor.*

#### **Definition 2.** *Let*

$$\mathcal{J}(\Omega) = [\mathcal{g}\_{i\_1,\ldots,i\_N}]\_\prime \quad \text{where} \quad \mathcal{g}\_{i\_1,\ldots,i\_N} = \begin{cases} m\_{i\_1,\ldots,i\_N} & \text{if } \qquad \omega\_{i\_1,\ldots,i\_N} = 1\\ \mathcal{O} & \text{otherwise} \end{cases} \tag{2}$$

*be a subtensor of* M *that contains only the entries pointed to by* Ω*.*

Thus: |Y(Ω)<sup>|</sup> <sup>=</sup> <sup>|</sup>Ω|. Let *<sup>ω</sup>* <sup>=</sup> vec(Ω) <sup>∈</sup> <sup>R</sup>∏*<sup>N</sup> <sup>n</sup>*=<sup>1</sup> *In* be a vectorized version of Ω.
