Image Completion with Hybrid Interpolation in Tensor Representation

The issue of image completion has been developed considerably over the last two decades, and many computational strategies have been proposed to fill-in missing regions in an incomplete image. When the incomplete image contains many small-sized irregular missing areas, a good alternative seems to be the matrix or tensor decomposition algorithms that yield low-rank approximations. However, this approach uses heuristic rank adaptation techniques, especially for images with many details. To tackle the obstacles of low-rank completion methods, we propose to model the incomplete images with overlapping blocks of Tucker decomposition representations where the factor matrices are determined by a hybrid version of the Gaussian radial basis function and polynomial interpolation. The experiments, carried out for various image completion and resolution up-scaling problems, demonstrate that our approach considerably outperforms the baseline and state-of-the-art low-rank completion methods.


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][10][11] or their hybrid versions [12] seem to be more appropriate. The hybrid strategies have been studied in many other research papers [13][14][15][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][18][19][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][25][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][33][34][35][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][40][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][49][50][51] on the use of various interpolation methods for solving image completion problems, other recent works [52][53][54][55][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.

Preliminary
Mathematical notations and preliminaries of tensor decomposition models are adopted from [31]. Tensors are multi-way arrays, denoted by calligraphic letters (e.g., Y). Let Y ∈ R I 1 ×...×I N be the N-way array. The elements of Y are denoted as y i 1 ,...,i N , where ∀n : 1 ≤ i n ≤ I n , 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 j contains the j-th column of Y. The symbol || · || F denotes the Frobenius norm of a matrix; || · || denotes the 2-nd norm. The sets of real numbers and natural numbers are represented by R and N , respectively. The symbols x and x stand for the floor and ceiling functions of x, respectively.
Let M = [m i 1 ,...,i N ] ∈ R I 1 ×...×I N be the N-th way observed tensor with missing entries, and Ω = [ω i 1 ,...,i N ] ∈ R I 1 ×...×I N be a binary tensor that indicates the locations of available entries in M. If m i 1 ,...,i N is observed, then ω i 1 ,...,i N = 1; otherwise, ω i 1 ,...,i N = 0. The locations of positive entries inΩ = 1 − Ω correspond to missing entries in M. The number of observed entries is equal to be a zero-filled incomplete tensor.
be a subtensor of M that contains only the entries pointed to by Ω.

Image Completion with Tucker Decomposition
For the N-th order tensor Y ∈ R I 1 ×I 2 ×...×I N , the Tucker decomposition [58] with the ranks J = (J 1 , J 2 , . . . , J N ) can be formulated as J n ] = [u i n ,j n ] ∈ R I n ×J n for n = 1, . . . , N is the factor matrix capturing the features across the n-th mode of Y. The operator × n stands for the standard tensor-matrix contraction along the n-th mode, which is defined as The optimization problem for low-rank image completion with the Tucker decomposition can be formulated as follows: where Y is given by model (3), and Φ(·) is a penalty function that imposes the desired constraints onto the core tensor G and the factor matrices otherwise, no changes. Assuming ∀n : J n < I n , the tensor Y in (3) has a low Tucker rank. Problem (5) can be solved by performing iterative updates with the Tucker decomposition in each step. The Tucker decomposition can be computed in many ways, depending on the constraints imposed on the estimated factors. If the nonnegativity constraints are used (as specified), any nonnegative least square (NNLS) solver can be applied in the alternating optimization scheme. Neglecting the computational cost of using an NNLS solver and the cost of computing the core G in (3), the total computational complexity for approximating the solution to (5) in K iterations can be roughly estimated as O(K ∑ N n=1 J n ∏ N p=1 I p ).

RBF Interpolation
The RBF interpolation [59] is a commonly used mesh-free method for approximating unstructured and high-dimensional data with high-order interpolants. Given the set of I distinct data points (x i , y i ) for i = 1, . . . , I, ∀i : x i ∈ R P , y i ∈ R, the aim is to find the approximating function y(x), which is referred to as the interpolant, satisfying the condition: ∀i : y(x i ) = y i . For the RBFs ψ : [0, ∞) → R, the interpolant takes the form where {w i } are real-value weighting coefficients. For I data points, we have The weighting coefficients {w i } can be computed from the system of linear equations: where y = [y 1 , . . . , (8) can be solved with any linear least-square solver.

Proposed Algorithm
For interpolation of N-way incomplete tensor M, the points x i and x j in (7) are expressed by index values: (6) can also be regarded in a wider sense, and hence the norm l 2 can be replaced with the distance function D(x i , x j ). For any distance function, the following conditions are satisfied: z). In the N-dimensional space, any data point y i 1 ,...,i N can be modelled with the interpolant where τ > 0 is a scaling factor. Let D([i 1 , i 2 , . . . , i N ], [j 1 , j 2 , . . . , j N ]) be an additively separable distance function, i.e., where d (n) (i n , j n ) is a one-variable function that expresses the distance metrics for the variables (i n , j n ) in the n-th mode. We also assume that the RBF ψ is multiplicatively separable. That is Considering separability conditions (10),(11), model (9) can be reformulated as follows: where f Following the standard tensor-matrix contraction rule given in (4), formula (12) can be presented in the following form: Remark 1. The model given in (13) has a form identical to model (3), where W is the core tensor, and {F (n) } are factor matrices. Hence, the RBF interpolation model in N-dimensional space, in which separability conditions (10), (11) are satisfied, boils down to the Tucker decomposition model, where the factor matrices {F (n) } are previously determined by rational functions.
If d (n) (i n , j n ) = i j n −1 n is expressed by an exponential function, τ = 1, and ψ(ξ) = −ξ, then f , and model (12) takes the form of the standard multivariate polynomial regression: where numbers (J 1 , . . . , J N ) determine the degrees of the polynomial with respect to each mode, i.e.,∀n : J n = D n + 1, where D n is a degree of the polynomial along the n-th mode. Thus the multivariate polynomial regression can be also presented in the form of the Tucker decomposition model. The RBF ψ can take various forms, such as Gaussian, polyharmonic splines, multiquadrics, inverse multiquadrics, and inverse quadratics [59]. The Gaussian RBF (GRBF), expressed by ψ(ξ) = exp{ξ}, is commonly used for interpolation; however, it cannot be used for constructing an interpolant with polynomial precision. For example, the GRBF does not approximate a linear function y(x) well. Hence, the polynomial regression in (14) is more suitable for linear or slowly varying functions, but it often yields underestimated interpolants if the polynomial has too low of a degree. An increase in the degree leads to an ill-conditioned regression problem and unbiased least squares estimates. To tackle this problem, both interpolation approaches can be combined, which leads to the following interpolation model: Model (15) can be equivalently presented in the tensor-matrix contraction form: where C = [c r 1 ,...,r N ] ∈ R R 1 ×...×R N is the core tensor of the polynomial regression model, and ∀n : System (16) has ∏ N n=1 J n + ∏ N n=1 R n variables, and only ∏ N n=1 I n equations, resulting in ambiguity of its solution under the assumption that ∀n : J n = I n . To relax this problem, a side-condition is imposed: Vectorizing models (16), (17) and performing the straightforward mathematical operations, we get n=1 J n , and c = vec(C) ∈ R ∏ N n=1 R n . The symbol ⊗ denotes the Kronecker product, and vec(·) means the vectorization operator.
To apply model (16) for image completion, letŶ be obtained according to Definition 1. System (18) can be applied to model non-zero entries inŶ, using the following transformations:ŷ = vec(Ŷ ) = Svec(Y ),F = SFS,P = SP, andŵ = Sw, where S = diag{vec(Ω)} ∈ R ∏ N p=1 I p ×∏ N p=1 I p is a diagonal matrix with binary values on the main diagonal. The matricesF andP have |Ω| zero-value rows, andF also has the same number of zero-value columns. After removing the zero-value rows and columns from the transformed system, the observed entries in M can be expressed by the model where y(ω) = vec(Y (Ω)) ∈ R |Ω| is a vectorized version of observed entries in M, F(ω, ω) ∈ R |Ω|×|Ω| is a submatrix ofF by selecting non-zero rows and columns, P(ω, :) ∈ R |Ω|×∏ N n=1 R n is obtained from P by removing all zero-value rows, and w(ω) = vec(W (Ω)) ∈ R |Ω| , taking into account rule (2). Note that |Ω| << ∏ N n=1 I n , if the number of missing entries in M is relatively large. The matrix F(ω, ω) is positive-definite because it is generated by a GRBF. The matrix P(ω, :) might be ill-conditioned if the polynomial functions of higher degrees are used, but the second term in (16) is used to better approximate linear relationships, and hence there is no need to use higher degrees. In this approach, ∀n : R n = 3, which leads to second-degree polynomials. The system matrix in (19) is therefore symmetric and positive-definite, and any least-square (LS) solver can be used to compute the vectors w(ω) and c from (19), given y(ω).
System (18) can also be used to compute the missing entries in M. Having the estimates for w(ω) and c, the missing entries can be obtained by solving the following system of linear equations: where F(ω, ω) ∈ R |Ω|×|Ω| is obtained from F by removing the rows and selecting the columns that are indexed by Ω, and the vectorȳ = y(ω) contains the estimates for the missing entries in M. Hence, the completed image is expressed by The system matrix in (19) has the order of |Ω| + ∏ N n=1 R n . Applying Gaussian elimination to (19), the computational complexity amounts to O (|Ω| + ∏ N n=1 R n ) 3 , and it is relatively large if the number of missing entries in M is small. It is therefore necessary to reduce the computational complexity for the LS problem associated with (19). To tackle this issue, let input tensorŶ = Ŷ (s 1 ,...,s N ) ∈ R I 1 ×...×I N be partitioned into small overlapping subtensors {Ŷ (s 1 ,...,s N ) }, where ∀n : s n = 1, . . . , S n and S n is the number of partitions ofŶ along its n-th mode. The total number of subtensors is equal to S Y = ∏ N n=1 S n . Each subtensor can be expressed byŶ (s 1 ,...,s N ) = y Γ(s 1 ),...,Γ(s N ) ∈ R L 1 ×...×L N , where ∀n : L n = I n S n ≤ I n . The set Γ(s n ) contains the indices of the entries inŶ which belong toŶ (s 1 ,...,s N ) along its n-th mode, and it can be expressed by ∀n : Γ(s n ) = {γ(s n ) + 1, . . . , γ(s n ) + L n } ∈ N L n for s n < S n {γ(S n ), . . . , I n } for s n = S n where γ(s n ) = (s n − 1)(L n − η n ) . The parameter η n = θ n L n 100 determines the number of overlapping pixels along the n-th mode, and θ n ∈ [0, 99] expresses the percentage of the overlap along the n-th mode.
The proposed methodology for image completion should be applied separately to each subtensor. The missing pixels are completed using a very limited volume of the input tensor but mostly restricted to their nearest neighborhood. It is thus a strategy that resembles PDE-based image inpainting, but it is much more flexible for highly dissipated known pixels and allows us to reduce the computational cost dramatically. Moreover, the factor matrices {F (n) } and {P (n) } in (16) are expressed by radial functions, and hence can be precomputed before using the subtensor partitioning procedure. In RBF-based interpolation methods, the samples or pixels that are close to the boundary of the region of interest are usually not well approximated. However, due to overlapping, boundary perturbation effects in our approach are considerably relaxed.
The pseudo-code of the proposed image completion algorithm is presented in Algorithm 1. Compute factor matrices F (n) using metrics d, and P (n) using set R.

Results
This section presents an experimental study that was carried out to demonstrate the performance of the proposed algorithms. The tests were performed for a few image completion problems using the following RGB images: Barbara, Lena, Peppers, and Monarch, which are presented in Figure 1. All of them have a resolution of 512 × 512 pixels.

Setup
The incomplete images were obtained by removing some entries from tensors representing the original images. The following test cases were analyzed: • A: 90 % uniformly distributed random missing tensor fibers in its third mode (color), which corresponds to 90 % missing pixels ("dead pixels"), • B: 95 % uniformly distributed random missing tensor entries ("disturbed pixels"), • C: 200 uniformly distributed random missing circles-created in the same way as in the first case, but the disturbances are circles with a random radius not exceeding 10 pixels, • D: resolution up-scaling-an original image was down-sampled twice by removing the pixels according to a regular grid mask with edges equal to 1 pixel. The aim was to recover the missing pixels on the edges.
We compared the proposed methods with the following: FAN (filtering by adaptive normalization) and EFAN (efficient filtering by adaptive normalization) [60], SmPC-QV (smooth PARAFAC tensor completion with quadratic variation) [38], LRTV (low-rank total-variation) [61], TMAC-inc (low-rank tensor completion by parallel matrix factorization with the rank-increasing) [62], C-SALSA (constrained split augmented Lagrangian shrinkage algorithm) [63], fALS (filtered alternating least-squares) [64], and KA-TT (ket augmentation tensor train) [65]. FAN and EFAN are based on adaptive Gaussian low-pass filtration. SmPC-QV performs low-rank tensor completion with smoothness-penalized CP decomposition and gradually increasing CP rank. LRTV accomplishes low-rank matrix completion using total variation regularization. TMAC-inc also belongs to a family of low-rank tensor completion, and in this approach, an incomplete tensor is unfolded with respect to all modes, and the resulting matrices are completed by applying low-rank matrix factorizations together with the adaptive rank-adjusting strategy. C-SALSA performs image completion using the variable splitting approach to solve an LS image reconstruction problem with a strong nonsmooth regularizer. fALS and KA-TT combine low-pass filtration with standard Tucker decomposition and tensor train models, respectively. The proposed algorithm is referred to as Tensorial Interpolation for Image Completion (TI-IC), and it is presented in Algorithm 1. It combines two strategies: RBF-interpolation with an exponential function, and multivariate polynomial regression. To emphasize the importance of both terms in the model (16), we also present the results obtained separately for each of them, using the same partitioning strategy in each case. The TI-IC algorithm with only the exponential term is referred to as TI-IC(Exp). When only the polynomial regression is used, TI-IC will be denoted as TI-IC(Poly).
TI-IC is flexible with respect to the choice of the distance function d (n) (·, ·), degrees of the interpolation polynomials, partitioning, and overlapping rates. Since {i n , j n } lie on a line, d (n) (i n , j n ) = |i n − j n | seems to be the best choice. Factor matrices {P (n) } were determined by quadratic polynomials, hence ∀n : R n = 3. Higher-order polynomials result in ill-conditioning of the system matrix in (19) and do not noticeably improve the performance. The partitioning and overlapping rates were set experimentally to [S 1 , S 2 , S 3 ] = [32,32,1] and [θ 1 , θ 2 , θ 3 ] = [33.33, 33.33, 0]. As the resolution of M is 512 × 512, the overlapping amounts to 5 pixels across the first and second mode for each subtensorŶ ∈ R 16×16×3 . For larger subtensors, computational time increased considerably, and we did not observe a noticeable improvement in the quality of recovered images. For smaller subtensors, the performance decreased. The scaling factor in the exponential RBFs was also determined experimentally, and to compute F (n) we set τ = 3 for TI-IC, and τ = 5 for TI-IC(Exp).
In the iterative algorithms, the maximum number of iterations was set to 1000, and the threshold for the residual error was equal to 10 −12 . The maximum rank was limited to 50.
The algorithms were implemented in MATLAB 2016a and run on the distributed cluster server in the Wroclaw Centre for Networking and Supercomputing (WCSS) (https://www.wcss.pl/en/) using PLGRID (http://www.plgrid.pl/en) queues and parallel workers. The resources were limited to 10 cores (ncpus) and 32 GB RAM (mem). The workers can be employed to run each algorithm for various initializations in parallel, or they can be used to process subtensorsŶ (s 1 ,...,s N ) in Algorithm 1 in parallel, in such a way that each subtensor is processed by one CPU core. The block partitioning procedure was implemented with the blockproc function in MATLAB 2016a, which has an option to run the computation across the available workers. We analyzed both options, i.e., when it was enabled and when it was disabled.

Image Completion
The recovered images were validated quantitatively using the signal-to-interference ratio (SIR) measure [31], defined as SIR = 20 log 10 The SIR values were averaged over the colormaps. The speeds of the algorithms were compared by measuring the averaged runtime of each test.  Figure 3-5, respectively. Due to the random initialization of some baseline algorithms, all the tests were repeated 100 times, and the SIR samples are presented in Figure 6 in the form of box-plots, separately for each test. The mean runtime of the evaluated algorithms and the corresponding standard deviations for each test case are listed in Table 1

SIR [dB]
Test C Test D

Discussion
The experiments were carried out for typical but challenging image completion problems. In test A, we knew only 10% of the pixels in the "Barbara" image, and the aim was to recover the 90% missing pixels. The results illustrated in Figure 2, 6A show that good quality reconstructions were obtained when the EFAN, SmPC-QV, fALS, KA-TT, TI-IC(Exp), and TI-IC algorithms were used, but the image recovered with TI-IC had the highest SIR score. TI-IC was also quite fast in this test (see Table 1). It lost in the runtime category only to FAN and EFAN, but its parallel version TI-IC * ) was more than 250 times faster than SmPC-QV. The latter performs the CP decomposition, but the difference in computational speed comes from the fact that in our method, the factor matrices are precomputed, and only the core tensor in the Tucker decomposition is estimated using the data. The results obtained in Test B are presented in Figure 3, 6B. They confirm the conclusions drawn from Test A, but it should be noted that TI-IC strengthened its leading position in the SIR performance. Moreover, its runtime was shorter than that in the previous test because only 5% of the entries were known, and hence the system matrix in (19) was smaller. Test C compared algorithms for the completion of many small-scale missing regions (holes), distributed across the image. The results presented in Figure 4, 6C show that EFAN and SmPC-QV failed to provide satisfactory reconstructions in this test, but TI-IC outperformed the other algorithms considerably. Obviously, a lower number of missing pixels in the image to be completed results in a noticeable increase in the runtime, but it was still below the runtime of low-rank tensor completion methods, such as SmPC-QV, TMac-inc, fALS, and LRTV. In Test D, only 50 % of pixels were unknown, but not all the tested algorithms handled this case well. In this test, TI-IC also yielded the best reconstruction (see Figure 5, 6D) but only slightly better than that obtained with TI-IC(exp). Hence, the low-degree polynomial regression in Test D does not affect the result considerably, and due to the computation time, it can be neglected. In other tests, both approaches (RBF interpolation and polynomial regression), combined appropriately, were essential to yielding high-quality results.

Conclusions
In this study, we showed the relationship between the models of RBF interpolation and Tucker decomposition (Remark 1). We combined the exponential RBF interpolation and polynomial regression in one model and experimentally demonstrated that such a hybrid method achieved the highest SIR scores in all the tests. The proposed algorithm (TI-IC) can be applied to a wide spectrum of image-completion problems. The incomplete images can contain many single missing entries or missing pixels distributed across the image, a large number of small-scale regions (holes), or regularly shaped missing regions, such as in resolution up-scaling problems. The TI-IC algorithm is also computationally efficient. It provides reconstructions of the highest quality and in a much shorter time than the tested low-rank tensor image-completion methods. Its runtime depends on the number of missing entries in an input tensor, and it is shorter if more entries are unknown. The computational complexity of the proposed method can be controlled by the block partitioning strategy, as proven in Remark 2. Assuming the overlapping in this partitioning strategy, we avoided visible disturbances around boundary entries of the blocks, which is an intrinsic effect of RBF interpolation methods. Furthermore, due to the use of RBFs, the overlapping blocks can be processed in parallel computer architectures, and our experiments demonstrated that the use of a parallel pool of workers in MATLAB considerably shortened the runtime of the proposed algorithm.
Summing up, the proposed algorithm outperforms all the tested image completion methods for a wide spectrum of tests. Its computational runtime is also satisfactory and considerably shorter than that for the low-rank tensor decompositions. The proposed algorithm can also be efficiently implemented on parallel computer architectures.

Conflicts of Interest:
The authors declare no conflict of interest.