An Improved Robust Fractal Image Compression Based on M-Estimator

: In this paper, a robust fractal image compression method based on M-estimator is presented. The proposed method applies the M-estimator to the parameter estimation in the fractal encoding procedure using Huber and Tukey’s robust statistics. The M-estimation reduces the inﬂuence of the outliers and makes the fractal encoding algorithm robust to the noisy image. Meanwhile, the quadtree partitioning approach has been used in the proposed methods to improve the efﬁciency of the encoding algorithm, and some unnecessary computations are eliminated in the parameter estimation procedures. The experimental results demonstrate that the proposed method is insensitive to the outliers in the noisy corrupted image. The comparative data shows that the proposed method is superior in both the encoding time and the quality of retrieved images over other robust fractal compression algorithms. The proposed algorithm is useful for multimedia and image archiving, low-cost consumption applications and progressive image transmission of live images, and in reducing computing time for fractal image compression.


Introduction
The idea of fractal image compression (FIC) was originally introduced by Barnsley et al. in 1987 [1], and the first automatic algorithm for FIC was developed by Jacquin in 1992 [2]. Due to the high compression ratio and simple decompression method, the FIC algorithm attracts many researchers' interests [3][4][5][6]. However, there are two drawbacks of fractal encoding: its complex computation and its retrieved image quality. Many fractal-coding algorithms have been developed to improve on these drawbacks [7][8][9][10][11][12]. Firstly, based on the basic concepts and methods of fractal image compression algorithms [13][14][15][16][17][18][19], a particle swarm optimization (PSO) method by utilizing the visual information of the edge property is proposed, and this method speeds up the encoder 125 times faster with only 0.89 dB decay of image quality in comparison to the full search method [20][21][22][23][24]. Then, several fractal image compression algorithms based on modified gray-level transform with fitting plane, characteristics of fractal and partitioned iterated function system, and wavelet transform with diamond search are proposed, which have better results in terms of compression speed and quality of results. Thereafter, high-efficiency video coding provided significant performance improvement in the compression ratio [25][26][27][28][29]. At the same time, various optimization algorithms [30][31][32][33], robustness [34][35][36][37], and feature extraction research [38][39][40][41] provide new ideas for the field of image compression.
However, little attention has been given to the robustness of such fractal-based methods until Ghazel et al. [42] proposed a fractal-based method to enhance and restore a noisy image. The method is based on the Additive White Gaussian Noise (AWGN), whose amplitude distribution follows a Gaussian distribution, while the power spectral density is uniformly distributed. Indeed, one of the original motivations for this study was the observation that a noisy image is somewhat denoised when it is fractally coded. Lu et al. [43] 2 of 12 presented an enhanced fractal predictive denoising algorithm for denoising the AWGN images by using a quadratic gray-level function. A new similarity measure for fractal encoding has been introduced by Jeng et al. [44], the so-called, Huber Fractal Image Compression (HFIC) schema. It brings the linear Huber regression technique from the robust statistics into the encoding procedure of the fractal image compression. Lin [45] proposed a similar technique of linear robust regression with least absolute derivation (LAD), least trimmed squares (LTS), and Wilcoxon to FIC.
The linear regression techniques are the attempts toward the design of robust fractal image compression [46][47][48]. However, there is a primary drawback to the current schema: the complexity of computations [49][50][51][52]. It needs almost ten thousand seconds for a full search HFIC algorithm to encode a standard gray-level Lena image. In our analysis, there are many calculations performed that are not required. For example, the fractal graphics compression technology based on a genetic algorithm [53,54] considers and utilizes the self-transformation characteristics of the image, but its disadvantage is that it requires an extensive search and is time-consuming. The algorithm proposed in this paper can reduce the compression time by effectively eliminating many unnecessary calculations. This approach is shown to significantly improve the time of compression procedures without affecting the quality of the reconstructed image. The proposed algorithm is a robust Mestimator-based fractal image compression algorithm (RMFIC). It is a full-search fractal encoding algorithm using the quadtree partitioning method and it uses the principle of linear robust statistics. The goal of the work is to improve the efficiency and the retrieved image quality of existing robust fractal encoding methods, such as HFIC and LAD-FIC. We chose two M-estimators for simulation; Huber's M-estimator [55] and Tukey's bisquare M-estimator [56].
The rest of the paper is organized as follows. Section 2 reviews the principle of original fractal image compression. Section 3 reviews the concepts of M-estimator and further describes Huber's case and Tukey's bisquare case. Section 4 describes the proposed RMFIC in detail. The experimental results are demonstrated in Section 5. Last, Section 6 presents the conclusion.

Baseline Fractal Image Compression
The fractal coding process can be described as follows: the original image is partitioned into non-overlapping blocks called range blocks (R) which are usually 8 × 8 pixels each. There are also domain blocks (D) which are twice the size of range blocks and overlap such that a new domain block starts every pixel. For each range block (r j ), a domain block (D t ) that best approximates this range block under the predetermined contractive transformations (ω i ), is searched for. In detail, these transformations contract the domain blocks spatially to 8 × 8 pixels and then transform the gray levels by a combination of rescaling, offset adjustment, reflection, and rotation. The general fractal affine transformation can be expressed as where γ i is the down-sampling operation. We denote d = γ i (D t ), and ϕ i can be expressed as where (t x , t y ) is the location of domain block D t in original image, s is the contract scaling, and o is the brightness offset. The 2 × 2 sub-matrix a b c d represents the isometrical transformations in Equation (3).
Hence, the general affined transformation can also be written as Equation (2) calculates the reconstructed (approximated) range block as the transformation of d, denoted as d.
The parameters are selected so that the transformed domain block (D t ) is a good approximation of the range block r j by minimizing the matching error where, d k = τ(d), k is the index of the isometrical transformation. The parameter s and o can be obtained by the least-squares method as The fractal code consists of t * x , t * y , k * , s * and o * . In practice, t * x , t * y , k * , s * and o * can be quantized by 8, 8, 3, 5, and 7 bits, respectively.
At decoder, the transformation for all range blocks can be iteratively applied to any initial image which then converges to the fractal approximation of the image.

Robust M-Estimator
The method of least squares is a prototypical M-estimator [57] since the estimator is defined as a minimum of the sum of squares of the residuals. However, the ordinary least square method is very sensitive to outlying observations. There are several approaches to deal with this problem. One approach is robust regression, which is to employ a fitting criterion that is not as vulnerable as least squares to unusual data. Statisticians have developed many robust methods based on different techniques. The most common general method of robust regression is M-estimation, introduced by Huber. Consider the general regression model for the i-th of n observations, where Here, ε i is an unobserved random variable, which is added to the linear relationship between the dependent variable y i and the vector of regressors x i .
The general M-estimator is defined to minimize the objective function where the function ρ(·) is a symmetric positive definite function with a unique minimum at zero and is chosen to be increasing slower than quadratically. Let ψ = ρ be the derivation of ρ. Differentiating the objective function with respect to the coefficients β, and setting the partial derivatives to 0, produces k + 1 a system of estimating equations for the coefficients Define the weight function w(x) = ψ(x)/x, and let w i = w(ε i ). Then the estimating equations can be written as This equation is a weighted least-squares problem, aiming at minimizing IRLS (iteratively reweighted least-squares) algorithm can be used. At the iteration t, the new weighted-least-squares estimates are computed by is the current weight matrix.

Huber's M-Estimator
Huber first introduced a robust M-estimator [58], named Huber's M-estimator. He chose the ρ-function to be formed by Equation (14) ρ H (e) = e 2 /2, for |e| ≤ k k|e| − e 2 /2, for |e| > k Here k is the tuning constant. In particular, k = 1.345σ for the Huber case [34], where σ is the standard deviation of the errors. A common approach [58] is to takê σ = median{abs(e − median(e))}/0.6745 Huber's M-estimator is a mixed l 2 and l 1 minimization problem, but LS is only a l 2 minimization problem (since in LS case, ρ(e) is always equal to ε 2 /2). When there are outliers in the measurements, generally the corresponding residuals are relatively large. Since ρ H (e i ) of Huber's M-estimator with n beyond 7 is smaller than that of LS, Huber's Mestimator is less sensitive to relatively large residuals. Hence it can suppress the influences of outliers.
From Equation (14), we obtain Hence, in the Huber case, the weight function, which we denote as w H (e), is defined by Equation (17) w H (e) = ψ(e)/e = 1, for |e| ≤ k e/k, for |e| > k Considering the general regression model (1), Let ε = e/σ and w i = w(ε i ) = w(e/σ), where σ is defined by Equation (15). The main steps of IRLS for the estimation of the parameter using the Huber M-estimator are described as Algorithm 1.

Algorithm 1 IRLS using Huber M-estimator
Input Initial value of estimatesβ (0) , random variableε i (0) , and weight w i 2: whileβ didn't converge to a convergence value do

5:
Compute the objective value E based on the objective function by Equation (10) 6: end while

Tukey's Bisquare M-Estimator Figures, Tables and Schemes
The Tukey's bisquare function can suppress the outlier even further. In Tukey's bisquare case, ρ-function is defined by Equation (18) as follows: where, the tuning constant t = 4.6851σ for Tukey's Bisquare case in particular [34], where σ is the standard deviation of the errors, and can be obtained by Equation (15).
The objective function is defined in Equation (10). Note that the Huber objective function increase without bound as the residual e departs from 0. In contrast, Tukey's bisquare objective function levels eventually level off (for |e| > k). The weights for the Huber estimator decline when |e| > k, and the weights for Tukey's bisquare are 0 for |e| > k.
Here W TB denotes the weight function of Tukey's bisquare case.
Considering the general regression model (1), Let ε = e/σ and w i = w(ε i ) = w(e/σ). The main steps of IRLS for the estimation of the parameter using Tukey's bisquare Mestimator are described as Algorithm 2.

5:
Compute the objective value E based on the objective function by Equation (10) 6: end while

Robust M-Estimator for FIC
To apply the general model (9) to the fractal coding process, a 2-D image block with L × L size can be rearranged as a 1-D vector, denoted as d 1 , d 2 , . . . , d L 2 . The model can be represented by Appl. Sci. 2022, 12, 7533 6 of 12 Compared with the general regression model, o and s are taken as the parameters to be estimated, denoted β = [o s] T as the parameter vector and the total number of regression data is equal to L 2 . Letβ = [ôŝ] T to be the estimated parameter vector. The residuals can be defined as follow.
The M estimateβ minimizes the objective function of Equation (10), where n = L 2 , ε = e/σ and σ is the tuning constant. The optimized scaling parameterŝ and offset parameterô can be obtained by the IRLS algorithm. We defined weights for residuals as follows: where, the w(·) function is defined by Equations (17) and (19), respectively, standing for Huber and Tukey's statistics. Thus, the estimated parameterβ can be calculated via Equation (13) aŝ Thus, we can obtain the estimated parameter in the IRLS iterated process directly bŷ In the HFIC algorithm proposed by Jeng et al. [19], the Huber M-estimator is employed in the scaling parameter and offset parameter estimation. The principle of HFIC is as follows: for each isometrical transformation, find the optimized scaling and offset parameters by minimizing the objective function of Equation (10). After all eight isometrical transformations are applied, record the optimized parameters with the minimized objective.
Note that HFIC is very time-consuming since there will be eight times of iterated IRLS algorithm for one pair of range block r and domain block d. The experimental results also confirm that it needs ten thousand seconds to finish a full-search Huber fractal coding. Here what we improved is eliminating the calculation during the obtaining of the optimized parameters to speed up the denoising FIC algorithm.
Although decreasing the domain's searching scope can speed up the compression process, it will have a negative effect on the quality of the reconstruction image. To show the performance of a robust estimator, in this study, we focus on the full-search algorithm. Namely, for a given range block r, all domain blocks will be applied to the fractal affine transformation to find the optimized one with the minimized matching error.
Quadtree partitioning is a method used to represent an image as a hierarchical quadtree data structure in which the root of the tree is the initial image and each node contains four subnodes [4]. A node represents a square image portion and its four subnodes correspond to the four quadrants of the square image portion. Quadtree-based FIC is the FIC algorithm used in quadtree partitioning. Since quadtree-based FIC is widely used in fractal compression algorithms, and it is already proven to be superior in performance over the fixed partitioning in the FIC algorithms, we chose quadtree-based FIC as the main frame of the proposed robust algorithm.
As in the analysis above, the most time-consuming parts of the original HFIC are the eight times iterated calculations on the parameter estimation. In the proposed algorithm, we decrease it to only once with almost no sacrifice on the quality. For each range block (r) and down-sampled domain block (d), in the preprocess stage, the four quadrants of the given block are ordered based on the brightness average value (r ). Then it is rotated to the one with the brightest quadrant in the upper left corner. After this preprocess, the matching error is compared between the range block (r ) and domain block (d ). Here r and its quadrants are all with the brightest quadrant in the upper left, as is domain block d . Hence, in the proposed algorithm, the IRLS algorithm is applied during the matching between r and d , instead of r and d. Hence, the eight times of parameters estimation for each r and d is changed to be only once, which is between r and d .
Algorithm 3 describes the main principle of the proposed decode algorithm.

Algorithm 3 The Robust M-estimator for FIC algorithm (RMFIC)
Input the original image, the tolerance value T, Output the fractal codes 1: Pre-decimate the image into domain image and partition the original image into blocks with size of L × L 2: Preprocess the blocks. Do flip or rotation to make the blocks with its brightest quad-rant in its upper left corner. Record the index k of the transformation for each block and its quadrants 3: while there is a block didn't encoded do 4: Fetch a range block r j 5: Obtain the optimal scaling parameter s and offset parameter o by IRLS algorithm and record its corresponding fractal codes 6: Calculate the matching error E i by Equation (10).

7:
if E i < T or the block size of the d i is already small enough then 8: go to 13 9: else 10: Partitioning the current block into 4 sub-blocks (quadrants) by size of L i /2 × L i /2 11: go to 5 12: end if 13: Record the fractal codes with the minimum matching error. 14: end while Two M-estimators are chosen, viz, Huber and Tukey's bisquare statistics. The respective IRLS algorithms are described in Sections 3.1 and 3.2.

Experimental Results
In order to evaluate the performance of this proposed algorithm (RMFIC), a computer simulation is performed. The standard gray-level image of Lena and Boats with 512 × 512 pixels size was used. The coding block size in the RMFIC is set to be 32 × 32, 16 × 16, or 8 × 8 using the quadtree partitioning method, and the tolerance value T = 8. For the proposed algorithm based on Huber's M-estimator (RMFIC-H), we adopt the object function of Equation (10) as the matching error function. The initial values of the estimated parameters are obtained by Least Squares. The tuning constant parameter k = 1.345σ, and in the study case, σ is defined to be fixed value 5. For the proposed algorithm based on Tukey's bisquare M-estimator (RMFIC-TB), we adopt MSE (mean-square-error) in Equation (25), as the matching error function. The initial values of the estimated parameters are obtained by Least Squares regression. The tuning constant parameter t = 4.6851σ, and σ is defined by Equation (15). We use a blank image as the initial image for decoding and the iteration number is 20, which is sufficient to reconstruct the image. The simulation is done under the Microsoft Visual C++ 2008 on Windows XP with SP3, Intel Core2 Duo 2.53 GHz CPU, and 2G RAM platform.
where f i is the grey value of the pixel i in the original image andf i is the gray level of the pixel i in the retrieved image. N denotes the total number of pixels.
The distortion between the original image f and the retrieved imagef is measured by the peak signal-to-noise ratio (PSNR). The definition of PSNR is PSNR( f ,f ) = 10 log 10 255 To show the quality of the reconstruction of the image, we corrupt the Lena image and the Boats image with salt and pepper noise with noise level as high as 10%. Figure 1 illustrates the retrieved images from the corrupted Lena image by full search FIC with normal least square regression (FIC-LS) and proposed RMFIC using Huber's and Tukey's robust statistics (RMFIC-H and RMFIC-TB). Figure 1a Figure 2 is the test case on Boats image which is corrupted by salt and pepper noise with a 10% noise level. Figure 2a,b show the original Boats image, and corrupted Boats image. Figure 2c shows the retrieved image from full search FIC with normal least squares regression with 22.68db by PSNR. Figure 2d,e illustrate the retrieved image from proposed methods (RMFIC) with corresponding qualities of, 25.57 db and 24.56 db by PSNR for Boats corrupted image. From the experimental results, we can form the conclusion that the RMFIC algorithm using both the Huber case and Tukey's bisquare case show the ability of robustness for FIC. Namely, it is much more insensitive to outliers, compared to LS-based full search FIC. Note that there are several tiny blocks not reconstructed correctly from the retrieved image from RMFIC-TB by vision. The existence of these tiny blocks is due to the definition of Tukey's bisquare M-estimator. From the weight Equation (19), we can see that if |e| > k, the weight value is zero. That means if the outlier is "detected", the bisquare estimator can get rid of this value from the data set. However, in the same situation, Huber's estimator has set the weight value to 1. Thus, the Tukey's bisquare M-estimator sacrifices some data. When encoding the small blocks by fractal, there is a case where pixels are taken as outliers and the weight is taken as zero. That is why these tiny incorrect retrieved blocks exist. Further optimization for this case is required in the future. Table 1 shows the simulation comparative results on Lena's image with salt and pepper noise of 5%, 10%, 15%, 20%, and 30% noise levels, respectively. Both these algorithms are with full search schema. As we discussed in Section 4, the most time-consuming part of the original HFIC is the eight times of iterated IRLS algorithm for one pair of range block r and domain block d. To encode the noisy Lena image in Figure 1b with the full search HFIC method, requires 9704 s. By eliminating the unnecessary calculations on the original HFIC algorithm, the encoding time of the proposed RMFIC-H (using Huber M-estimator) was reduced almost 87.5%, and its relative reconstructed image quality is improved since the quadtree partitioning method is adopted. Even when the noise level is up to 30%, the proposed RMFIC-H method still shows a better quality by PSNR. RMFIC-TB (using Tukey's bisquare M-estimator) also shows considerable robustness, compared to the original FIC algorithm with normal least squares regression. Compared to the existing robust fractal-based algorithms, such as HFIC and LAD-FIC, RMFIC-TB shows a similar performance on the reconstructed image from the noisy image. However, the encoding time of RMFIC-TB is not optimized enough, and further work is still required in the future. as zero. That is why these tiny incorrect retrieved blocks exist. Further optimization for this case is required in the future.  as zero. That is why these tiny incorrect retrieved blocks exist. Further optimization for this case is required in the future.  For the bell-shaped noises, such as zero-mean Gaussian noise and Laplace noise, the original HFIC algorithm does not show significant robustness under them. We take a simple examination for zero-mean Gaussian noise and the Lena image of 10% noise level to test, by the baseline FIC (LS). And the reconstructed image can obtain PSNR of 26.48, by RMFIC-H algorithm, the PSNR is 26.40 and by RMFIC-TB, this value is 26.43.

Conclusions
In this paper, an improved robust fractal encoding method based on M-estimator has been proposed (RMFIC). This proposed algorithm applies the maximum likelihood estimation technique to the quadtree-based fractal image compression procedure. Improvement is made by reducing the unnecessary computations for the parameter estimation procedure. Both Huber's and Tukey's bisquare M-estimators are used in this work. The proposed methods show good robustness against salt and pepper noise. Comparative experimental results are given for the proposed method and the existing robust fractal encoding method. The encoding time of RMFIC-H is reduced by almost 87.5%, compared with the original full search HFIC, and its retrieved image quality is better than other robust algorithms such as HFIC and LAD-FIC. The RMFIC-TB also shows considerable robustness.
However, both M-estimators do not make a significant improvement in image quality for Gaussian noise, so further research is still required in the future. In addition, the proposed method only considered the normal raw inputting image, and then added the noise. However, the problem of using the raw corrupted blurred image as the original image has gradually attracted the attention of scholars and has become an important research topic of image compression recently. This will be another topic that the proposed method will further study.