Image Completion in Embedded Space Using Multistage Tensor Ring Decomposition

Tensor Completion is an important problem in big data processing. Usually, data acquired from different aspects of a multimodal phenomenon or different sensors are incomplete due to different reasons such as noise, low sampling rate or human mistake. In this situation, recovering the missing or uncertain elements of the incomplete dataset is an important step for efficient data processing. In this paper, a new completion approach using Tensor Ring (TR) decomposition in the embedded space has been proposed. In the proposed approach, the incomplete data tensor is first transformed into a higher order tensor using the block Hankelization method. Then the higher order tensor is completed using TR decomposition with rank incremental and multistage strategy. Simulation results show the effectiveness of the proposed approach compared to the state of the art completion algorithms, especially for very high missing ratios and noisy cases.


INTRODUCTION
Joint analysis of datasets recorded by different sensors is an effective approach for investigating a physical phenomenon. For this propose, the acquired datasets can be recorded in different matrices and jointly analyzed (Ozerov and Févotte, 2009).
By increasing the volume of the recorded data and also the number of the recorded signals, matrices are not any more efficient for analyzing and representing the datasets. For overcoming this limitation, the acquired datasets are represented in tensor formats. Tensors are higher order arrays used in different signal processing applications such as Blind Source Separation (BSS) (Cichocki et al., 2009;Bousse et al., 2016), image in-painting/completion (Liu et al., 2012;Yokota et al., 2016;Wang et al., 2017), time series analysis (Kouchaki et al., 2014;Sedighin et al., 2020) and machine learning (Rabanser et al., 2017;Sidiropoulos et al., 2017).
The acquired datasets can be incomplete or uncertain due to different reasons such as noise, outliers, hardware failure or human error (Liu et al., 2012;Yokota et al., 2016;Wang et al., 2017). Even, sometimes, the resolution of the recorded data is not sufficient due to hardware limitations. In these situations, tensor completion can be very helpful for efficient data analysis.
Tensor completion is an approach for recovering missing or uncertain elements of a data tensor using its available elements (Liu et al., 2012;Yokota et al., 2016;Wang et al., 2017). At the first look, the problem seems to be ill-posed and difficult, but indeed, it can be done using some assumptions for the original data tensor (Gandy et al., 2011;Liu et al., 2012;Bengua et al., 2017).
Generally, tensor completion approaches can be divided into two groups: low rank based approaches (Liu et al., 2012;Bengua et al., 2017) and tensor decomposition based approaches (Yokota et al., 2016, Yokota et al., 2018. However, in some completion methods such as (Yuan et al., 2018b), both of the approaches have been combined.
In low rank based completion approaches, it is assumed that the original tensor is low rank. Under this assumption, the completion is to find the missing elements of a tensor in a way that the completed tensor has low rank approximation (Signoretto et al., 2010;Gandy et al., 2011;Liu et al., 2012;Bengua et al., 2017;Yuan et al., 2018b). In contrast to matrices, there is no unique definition for the tensor rank. One of the basic definitions for the tensor rank is the Canonical Polyadic Decomposition (CPD) rank, which is the number of rank-one tensors after CP decomposition of the original tensor. Estimating the CPD rank of a tensor is difficult, therefore, other alternative definitions such as weighted summation of the ranks of different unfoldings of the data tensor have been proposed (Liu et al., 2012). Algorithms such as Simple Low Rank Tensor Completion (SiLRTC) or High accuracy Low Rank Tensor Completion (HaLRTC) (Liu et al., 2012;Bengua et al., 2017) can be considered as examples for rank minimization based completion algorithms.
In tensor decomposition based approaches, the incomplete tensor is decomposed into several lower order or smaller factor matrices and/or core tensors. These latent factors usually preserve the structural information of the original tensor. Therefore, the completed tensor can be reconstructed using the estimated latent factors. Existing different approaches for tensor decomposition, result in different decomposition based completion approaches. In (Yokota et al., 2016), a completion approach based on CPD has been proposed. In (Yokota et al., 2018), Tucker decomposition has been exploited for tensor completion.
Besides the mentioned tensor decomposition based approaches, recently, new completion methods based on Tensor Train (TT) (Grasedyck et al., 2015;Ko et al., 2018;Yuan et al., 2019) and Tensor Ring (TR) (Wang et al., 2017;Yuan et al., 2018a;Yuan et al., 2018b) decompositions have also been proposed. TT decomposition is an approach for decomposing a tensor into a series of inter-connected third order core tensors (Oseledets, 2011), as (see Figure 1).
where X ∈ R I1×I2×/×IN is the N-th order tensor to be decomposed and G (n) ∈ R Rn−1×In×Rn is the n-th core tensor. The vector [R 0 , R 1 , . . . , R N ] is the TT rank vector and R 0 R N 1 (Oseledets, 2011). The (i 1 , . . . , i N )-th element of X, i.e., x i 1 ,...,i N , can be estimated as where G (n) i n is the i n -th lateral slice of the n-th core tensor. In TR decomposition, similar to TT decomposition, a tensor is decomposed into a chain of inter-connected third order core tensors with a notation similar to (Eq. 1), with, R 0 R N ≥ 1 . The elements of a tensor with TR structure can be computed as where "tr" denotes the trace operator.
Recently, a promising approach, called Hankelization, has been developed for improving the quality of tensor completion approaches (Yokota et al., 2018;Sedighin et al., 2020). Previously, Hankelization has been used as an initial step for a time series analysis algorithm, called Singular Spectrum Analysis (SSA) (Golyandina et al., 2013;Rahmani et al., 2016;Hassani et al., 2017;Kalantari et al., 2018). Generally speaking, Hankelization is an approach for transferring a lower order dataset to a higher order one by exploiting Hankel matrix structure. In this approach, the initial data is reshaped into a tensor whose slices have Hankel or block Hankel structure. Recall that in a matrix with Hankel structure, all of the elements in each skew diagonal are the same. Hankelization provides possibility of exploiting local correlations of the elements. Due to the Hankel structure, it is expected that the resulting higher order tensor provided by the Hankelization is low rank. This low-rankness in addition to a rank incremental strategy for Tucker model has been used for tensor completion in (Yokota et al., 2018). In (Asante-Mensah et al., 2020), TR decomposition has been used for the completion of Hankelized images. Also in (Sedighin et al., 2020), a new approach for time series completion based on a two-stage Hankelization has been proposed. In this approach, a one way data, i.e., the time series, is first reformatted into a Hankel matrix, and then the resulting Hankel matrix is block Hankelized into a 6-th order tensor using block (patch) Hankelization. In block Hankelization approach, in contrast to classic Hankelization methods, blocks of elements are Hankelized instead of individual elements (Sedighin et al., 2020).
In this paper, we mainly focus on image completion. Using the block Hankelization approach of (Sedighin et al., 2020), the original incomplete tensor (image) is first transformed into a higher order 7-D tensor. Then a TR decomposition has been applied for image completion. Applying TR decomposition for the completion of Hankelized images, or better to say, in the embedded space, has been investigated in few papers such as (Asante-Mensah et al., 2020) and (Sedighin et al., 2021). Similar to many TT and TR completion approaches, determining proper ranks for completion is an important issue, so in this paper, the rank incremental strategy, developed by us in (Sedighin et al., 2021), is used for determining the TR ranks. Moreover, a multistage strategy is developed for further improving the quality of reconstructed images. Multistage strategy has been previously proposed by us in (Sedighin et al., 2020) for time series completion. To best of our knowledge, this is the first time of using multistage strategy for image completion along with adaptive TR decomposition. This is the main difference of the proposed approach with the simulation presented in (Sedighin et al., 2021). The contributions of this paper can be summarized as.
• Applying block Hankelization on each frontal slice of the incomplete image which results in a 7-th order tensor. The block (patch) Hankelization has been applied on each slice of the color image separately, which results in three 6-th order tensors. Then these 6-th order tensors have been combined with each other to produce a 7-th order tensor • Applying TR decomposition with automatic rank incremental strategy proposed in (Sedighin et al., 2021). • Developing a multistage strategy for improving the quality of the final reconstructed image.
Extensive simulation results confirmed the quality of the proposed algorithm in image completion, especially for high missing ratios and noisy cases.
The rest of this paper is organized as follows: Notations and preliminaries are presented in Section 2, different TT and TR completion algorithms are briefly reviewed in Section 3. Section 4 is devoted to the Block Hankelization and the proposed algorithm is presented in Section 5. Finally, results and discussion are presented in Section 6 and Section 7, respectively.

NOTATIONS AND PRELIMINARIES
Notations used in this paper are adopted from . It is assumed that vectors, matrices and tensors contain real valued elements. An N-th order tensor is denoted by a bold underlined capital letter as X ∈ R I 1 ×I 2 ×/×I N , where I n is the number of elements in the n-th mode. An I 1 × I 2 matrix is denoted by a bold capital letter as X ∈ R I1×I2 and vectors are denoted by bold lower case letters as x ∈ R I1 . The Matricization or unfolding of a tensor is reshaping that tensor into a matrix with the same elements. The mode-n matricization of a tensor is defined as X (n) ∈ R I n ×I 1 ···I n−1 I n+1 ···I N in which Mode-{n} canonical matricization of a tensor is also defined as X <n> ∈ R I 1 I 2 ...I n ×I n+1 ...I N in which X <n> (i 1 · · · i n , i n+1 · · · i N ) x i1,i2,...,iN . Furthermore, X [n] ∈ R I n ×I n+1 /I N I 1 /I n−1 is an unfolding of a tensor in a way that Hadamard (component wise) product of two tensors is denoted by ;. Mode-n product of a tensor A and a matrix B is denoted by A × n B. vec(X) and rank(X) are the vectorization and rank of a matrix, respectively.

TENSOR TRAIN AND TENSOR RING BASED COMPLETION ALGORITHMS
Tensor Train decomposition has been exploited in many completion algorithms. In (Yuan et al., 2019), two completion algorithms, TT Weighted OPTimization (TT-WOPT) and TT Stochastic Gradient Descent (TT-SGD) have been proposed. The two algorithms are based on the minimization of the following weighted cost function using a gradient descent method: where Ω is a binary mask tensor whose elements are 1 for the observed and 0 for the missing elements, X is the observed tensor and X(θ) is the estimated tensor with TT structure as and θ G (1) , . . . , G (N) . A TT based completion algorithm using system identification has been proposed in (Ko et al., 2018). In this approach, the indices of the observed elements are the inputs of a system and the observed values are the outputs. Alternating Least Squares (ALS) and Alternating Directions Fitting (ADF) based approaches have also been proposed in (Grasedyck et al., 2015).
TT rank minimization has been used in many of tensor completion algorithms. Indeed, these algorithms are basically categorized as the rank minimization based approaches, however, since they are using TT rank which is closely related to TT decomposition, we have briefly reviewed them in this section.
In these algorithms, the cost function is defined as (Bengua et al., 2017) where α n is the weight of the n-th term. The above cost function is a weighted summation of the ranks of the canonical unfoldings of the input tensor, called TT rank. Algorithms, such as SiLRTC-TT and TMac-TT have been developed based on the minimization of the TT rank and parallel matrix factorization (Bengua et al., 2017). Furthermore, minimization of the TT rank in addition to a sparsity assumption of the mode-n matricizations of the incomplete tensor has been proposed in (Yang et al., 2018). TR decomposition has also been used in many completion algorithms. TR-WOPT has been proposed in (Yuan et al., 2018a) and is based on minimizing the cost function similar to (Eq. 4), where X(θ) has a TR structure. In (Wang et al., 2017), an approach called TR-ALS (Alternating Least Squares) has also been developed for tensor completion.
TR-LRF (TR Low Rank Factors) algorithm has been proposed in (Yuan et al., 2018b). This algorithm is based on a combination of rank minimization and TR decomposition approaches, in which the rank minimization has been applied on different matricizations of the core tensors resulting from TR decomposition of the incomplete tensor.
Moreover, in , a completion algorithm called MTRD has been introduced based on a more balanced matricization of a tensor with TR structure. A new matricization has been defined as X <n,d> ∈ R I a I a+1 ···I n ×I n+1 I n+2 ···I a−1 , where and a is defined as In (Yu et al., 2019), another completion approach, called TRNNM (TR Nuclear Norm Minimization), based on minimizing the nuclear norm of X <n,d> has been proposed. Other completion approaches based on balanced matricizations of tensors with TR structure have also been proposed in (Huang et al., 2020a,b).

BLOCK HANKELIZATION
As mentioned earlier, Hankelization is an effective approach in signal processing for exploiting local correlations of pixels or elements. In many of time series completion or forecasting algorithms, Hankelization has been used for transforming a lower order signal into a higher order matrix or tensor ). An illustration for Hankelizing a time series has been shown in Figure 2.
In (Yokota et al., 2018), Hankelization is performed by multiplying a special matrix, called duplication matrix, by each of the modes of the original tensor followed by a tensorization step.
Block (patch) Hankelization using blocks of elements (instead of single elements) has been proposed in (Sedighin et al., 2020). The block Hankelization is also performed by multiplying matrices, called block duplication matrices (shown in Figure 3), by different modes of the original tensor (Sedighin et al., 2020). Each block duplication matrix is of size S H,k ∈ R PT k (Ik /P−Tk +1)×Ik , where S H,k is the block duplication matrix which is multiplied by the k-th mode of the original tensor, P is the block size, T k is the k-th window size and I k is the size of the k-th mode of the original tensor which should be dividable by P.
The resulting matrix after multiplication of matrix H ∈ R I1×I2 by block duplication matrices S H,1 ∈ R PT1(I1/P−T1+1)×I1 and The matrix H b is then folded (tensorized) into a 6-th order tensor of size P × P × T 1 × D 1 × T 2 × D 2 , where P × P is the block size, T 1 and T 2 are the window sizes for the first and second modes, respectively, D 1 I 1 /P − T 1 + 1 and D 2 I 2 /P − T 2 + 1. Block Hankelization has been illustrated in Figure 4. Note that in this figure, each colored block is a P × P matrix [for more details please see (Sedighin et al., 2020)].

PROPOSED ALGORITHM
As mentioned earlier, the main idea of this paper is to apply TR decomposition on the Hankelized incomplete image. In many TR completion algorithms, the cost function used for completion can be written as where Ω is the mask tensor, whose elements are 1 for the observed and 0 for the missing elements of the input tensor, X is the incomplete input tensor, X is the completed tensor with TR representation and θ G (1) , G (2) , . . . , G (N) . The cost function in the embedded space, i.e., after Hankelization of the datasets, can be written as where Ω H and H are the resulting tensors after Hankelization of Ω and X, respectively and H(θ) is the TR estimation of H. FIGURE 3 | Illustration of a block duplication matrix. The lower order tensor is transformed into a higher order tensor by multiplying this matrix into its different modes following by a folding step.
Frontiers in Artificial Intelligence | www.frontiersin.org August 2021 | Volume 4 | Article 687176 Similar to (Yokota et al., 2018) and by using a Majorization Minimization approach, minimization of (Eq. 10), can be achieved by minimizing the following auxiliary function where H(θ k ) is the estimated tensor in the TR format in the k-th iteration and 1 is a tensor with the size equal to Ω H whose all elements are 1. The cost function (Eq. 11), can be re-written as (Yokota et al., 2018;Sedighin et al., 2020) whereH Ω H ;H + 1 − Ω H ; H(θ k ) . It can be inferred from (Eq. 12) that the minimization of (Eq. 10) is equivalent to the TR decomposition of the input tensor, where in each iteration, the observed elements are kept fixed. Considering (Eq. 12), TR estimation is done in an iterative manner, where the number of iterations is denoted by I int . In each iteration, the TR decomposition ofH, i.e., H(θ) has been estimated and then the updatedH is computed. As mentioned in the introduction, determining the ranks for TT and TR decompositions is a challenging task. In papers such FIGURE 4 | Illustration of the block Hankelization. In the first step, the initial matrix is multiplied by two block duplication matrices, and in the second step, the resulting matrix is folded into a 6-th order P × P × T 1 × D 1 × T 2 × D 2 tensor, where each colored block is a P × P matrix. as (Wang et al., 2017;Yuan et al., 2018b), the TR ranks are determined in advance and as inputs for the algorithm. These fixed rank approaches usually fail in reconstruction of the images with high missing ratios. Also, determining proper input ranks is difficult. In papers such as (Yokota et al., 2018;Sedighin et al., 2020), the input ranks are not determined in advance, but are increased gradually during iterations. Experiments show that these adaptive rank approaches are more successful in image completion when missing ratio is high. In this paper, for controlling the TR ranks, i.e., the size of the core tensors, rank incremental approach, presented in detail in (Sedighin et al., 2021), has been exploited. In (Sedighin et al., 2021), sensitivities of the estimation error to each of the core tensors are estimated and the size of the core tensors whose sensitivities are larger than a threshold are increased. However, in contrast to (Sedighin et al., 2021), in this paper, both of the ranks of the selected core tensor are increased. The step for increasing the ranks in each iteration can be 1 or any other natural number. It is also possible not to increase the ranks in some iterations, i.e., changing the rank increasing step to 0. The latter case is mainly applicable for images with high missing ratios. The estimated H is then re-transformed into the original image space by de-Hankelization, which results in X. De-Hankelization is a procedure for transferring back a Hankelized dataset into its original format. De-Hankelization of a block Hankelized tensor is done by averaging the corresponding frontal slices (slices with the same color in Figure 4) and then averaging the blocks with the same color

Algorithm 1 | Pseudo code of the proposed algorithm
INPUT: Incomplete image X, mask tensor Ω, vector of block sizes p [P 1 , P 2 , . . . , P L ] where L is number of stages, t [T 1 , T 2] , rank vectors r 1 , . . . , r N , error level ϵ, maximum value of the rank R max , internal iteration number I int , the rank increasing step (inc) and the threshold for selecting the core tensors for rank incremental in each iteration (tol). OUTPUT: Completed image X. 1: Initialize the missing elements of X by zero. 2: for l 1 : L do 3: Block Hankelize the input incomplete image (X) and the mask tensor (Ω) by block size P l × P l and window size t [T 1 , T 2] which results in H and Ω H . 4: PutH H 5: while max(r l ) < R max (or the normalized approximation error is higher than the error level ϵ) do 6: for j 1 : I int do 7: Compute the TR decomposition ofH, i.e., H with rank vector r l . 8:H Ω H ; H + (1 − Ω H ) ; H 9: end for 10: Increase the elements of the rank vector r l using the approach of (Sedighin et al., 2021). 11: Compute X by de-HankelizingH (in noisy cases by de-Hankelizing H).

13:
Apply smoothing by replacing each estimated element (for Ω 0) by the average of its four neighbors in the frontal slice and keeping the observed elements fixed. 14: end while 15: Put X X.

16: end for
Frontiers in Artificial Intelligence | www.frontiersin.org August 2021 | Volume 4 | Article 687176 and alphabet (Sedighin et al., 2020). An illustration for de-Hankelization has been shown in Figure 5. The estimated X is then modified as X Ω ; X + (1 − Ω) ; X and smoothness, by replacing each estimated element of X (for Ω 0) by the average of its four neighbors in the frontal slice, is applied. The above procedure is repeated until the maximum value of the TR rank vector achieves the limit or the normalized approximation error in the embedded space i.e., (or the change in the normalized approximation error between two consecutive iterations) becomes less than a pre-determined error level ϵ.
In this paper, we have employed a multistage strategy to further improve the quality of final results. This is the main difference of this paper with the simulation presented in (Sedighin et al., 2021), which was the result of a single stage algorithm. The main idea and the structure of the multistage strategy have been illustrated in Figure 6. In this approach, the incomplete image is block Hankelized with block size P 1 × P 1 and completed by the algorithm. Then the output of the algorithm is again processed as the new input of the algorithm with block size P 2 × P 2 , where P 2 < P 1 . The initial TR rank vector for the first stage is set to all 1 vector. For the next stages, the initial values of the rank vectors are set higher than 1 but lower than the final ranks of the previous stage. The procedure can be repeated for several stages until the desired accuracy is achieved. The multistage strategy has a benefit of providing a good initialization for the algorithm and improving the quality of the final result (especially for very high missing ratios and noisy cases), in comparison to a single stage algorithm.
The pseudo code of the proposed algorithm has been listed in Algorithm 1.

RESULTS
In this section, the effectiveness of the proposed multistage approach has been investigated and compared with the state of the art algorithms. The proposed algorithm has been compared with the HaLRTC algorithm (Liu et al., 2012) which is a rank minimization based completion algorithm and also with MDT (Yokota et al., 2018), TR-ALS (Wang TABLE 2 | Comparison of the performance of the algorithms for the reconstruction of images with 90% missing ratio. PSNR, SSIM and normalized approximation error corresponding to each image have been written in brackets.

Missed image
HaLRTC TT-WOPT TR-ALS TR-LRF MDT Proposed  (Yuan et al., 2018a) and TR-LRF (Yuan et al., 2018b), which are tensor decomposition based approaches. We have used standard (typical) color images for our simulations.
In the first simulation, the approaches have been compared for the completion of images with high missing ratios, i.e., 90%, 95% and 99% missing pixels. The parameters regarding to each algorithm have been listed in Table 1. The results have been presented in Tables 2-4, respectively. Peak Signal to Noise Ratio (PSNR), Structural SIMilarity (SSIM) and normalized approximation error ||T− X|| F ||T|| F , where T is the original image, have been written beneath each figure in brackets. For the proposed approach, the two stage strategy has been exploited. The block size of the first stage was set to P 1 8 and for the second stage was set to P 2 4. The window size of the both stages was t [2, 2]. Rank incremental step for 90% and 95% missing ratios were 1 for both of the stages. For 99% missing ratio, in the first stage, the ranks were increased by 1 in each 4 iterations, and for the second stage the rank increasing step was 1 in each iteration. The size of the reconstructed images in Table 2 and Table 3 are 128 × 128 × 3 and in Table 4    been set to [15,15,3]. In TR-ALS, the initial rank for 90% missing ratio was set to [3, 3, 3] and for 95% missing ratio was set to [25,25,25] (Table 1). As the results show, for high missing ratios, the proposed approach provides higher performance in comparison to the fixed rank approaches such as TT-WOPT, TR-ALS and TR-LRF. This higher performance is more significant for 99% missing ratio where TR-ALS completely fails. Indeed, the TR-ALS algorithm could not provide any output for images with 99% missing ratio. This shows the importance of adaptive rank selection for image completion. The higher performance of the algorithm is also clear in comparison with the low rank completion based algorithm such as HaLRTC. The results show that the performance of HaLRTC dramatically deteriorates when the missing ratio increases. The proposed algorithm has also slightly higher performance in comparison to the MDT approach, which is a tensor decomposition based approach with rank incremental strategy. The run time of our algorithm is larger than the other considered approaches. This is due to the multistage structure of our algorithm. Additionally, our algorithm in each stage has inner iterations for computing TR decomposition of the input and outer iterations for increasing the ranks. These two kinds of iterations are repeated for each stage. Moreover, in our approach 3-rd order images are transformed into 7-th order tensors. Analyzing this higher order tensor also increases the run time of the algorithm. These together increase the run time of our approach, however this longer run time results in the better performance of our approach. In the next simulation, the algorithms have been compared for the completion of structurally missing images. The images with blocks and slices missing elements were completed by the algorithms and the results have been presented in Table 5. The image size of this simulation was 128 × 128 × 3 and the algorithms parameters were similar to Table 1 for 90% missing ratio. The completed images along with the PSNR's, SSIM's and normalized approximation errors, show the superiority of the proposed approach in comparison to the other completion methods.
In the next simulation, the effectiveness of the multistage strategy has been investigated. For this purpose, one stage  algorithms with P 8 and P 4 have been compared with a two stage algorithm with P 1 8 and P 2 4 for the completion of an image with 99% missing ratio. The window size was set to t [2, 2] and the initial rank vector for the second stage was set to [10, 10, . . . , 10]. The results have been illustrated in Table 6. As the results show, the quality of the reconstructed image resulted from the two stage algorithm is higher than the images resulted from single stage algorithms. This is because the first stage of the algorithm can provide a good initialization for the second stage, which can result in a better completion performance comparing to single stage algorithms. Finally, for investigating the effect of noise on the performance of the proposed algorithm, several incomplete noisy images have been completed by the proposed and the MDT approaches. The size of the images in this simulation was 128 × 128 × 3. The missing ratio of the images was 90% and the remaining pixels have been contaminated by noise with normal distribution and standard deviation σ. The two methods have been compared for different σ's. The parameters of the two algorithms were similar to Table 1 for 90% missing ratio. The ranks of the second stage of the proposed algorithm have been set equal to 10 for σ 0.1 and 3 − 6 for other σ's. For the MDT approach, noise levels have been selected as 10 −2 , 4 × 10 −2 , 10 −1 , 1.6 × 10 −1 , 2.5 × 10 −1 and 3.6 × 10 −1 for σ 0.1, 0.2, 0.3, 0.4, 0.5 and 0.6, respectively. For each σ, seven incomplete noisy images have been completed by each algorithm and the resulting PSNR's have been averaged and the standard deviation has been computed. Note that in noisy cases, the line 12 of the Algorithm 1 moves to after line 13 and the output of the line 13 will be the output of each stage (i.e., before replacing the observed elements). The results have been presented in Table 7. As the results show, the qualities of the both approaches decrease by increasing σ. However, even in these noisy cases, the quality of the proposed approach is higher than MDT.

DISCUSSION
A new approach for image completion in the embedded space by block Hankelization and by exploiting TR decomposition has been proposed and extensively tested. In this approach, the incomplete image has been transformed into a higher order 7-th order tensor using block Hankelizaion and then a TR decomposition with rank incremental strategy and smoothness have been applied for completion. Moreover, a multistage strategy, which has been previously applied by us for time series completion, has been used for increasing the quality of the final completed image. Simulation results and comparisons with the state of the art algorithms indicated the advantage of the proposed algorithm.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.