Fractal Image Coding Based on a Fitting Surface

. A no-search fractal image coding method based on a fitting surface is proposed. In our research, an improved gray-level transform with a fitting surface is introduced. One advantage of this method is that the fitting surface is used for both the range and domain blocks and one set of parameters can be saved. Another advantage is that the fitting surface can approximate the range and domain blocks better than the previous fitting planes; this can result in smaller block matching errors and better decoded image quality. Since the no-search and quadtree techniques are adopted, smaller matching errors also imply less number of blocks matching which results in a faster encoding process. Moreover, by combining all the fitting surfaces, a fitting surface image (FSI) is also proposed to speed up the fractal decoding. Experiments show that our proposed method can yield superior performance over the other three methods. Relative to range-averaged image, FSI can provide faster fractal decoding process. Finally, by combining the proposed fractal coding method with JPEG, a hybrid coding method is designed which can provide higher PSNR than JPEG while maintaining the same Bpp.


Introduction
The basic idea of fractal image coding was first proposed by Barnsley [1].Its kernel issue is to find an iterated function system whose fixed point can approximate the input image well.Then, Jacquin [2] introduced the first practical fractal coding scheme.After many years of development, it has been successfully used in many image processing applications such as image compression [3], image denoising [4,5], image retrieval [6,7], image magnification [8,9], and image watermarking [10][11][12].
Although the fractal coding method has the advantage of potential high compression ratio, resolution independence, and fast decoding, its main drawback is very time-consuming in the encoding process.In order to reduce the computational complexity of fractal encoding, converting the global search to local search is an effective way to solve this problem.It mainly consists of classification techniques and feature vector techniques.For the former ones, the range blocks and domain blocks are first divided into different categories and the block matching is only carried out within the same category [13,14].For the latter ones, the block matching process is carried out in the feature space.Due to the lower dimension of features and more effective searching strategies, such as kd-tree method, the fractal encoding process can be finished in a short time [15][16][17].Typically, Furao and Hasegawa [18] introduced the first no-search fractal encoding algorithm which can provide faster encoding process and higher compression ratio at the expense of poor quality of the decoded image.Afterwards, Wang et al. [19,20] proposed the adaptive plane and fitting plane to improve the quality of the decoded image, respectively.In particular, the latter one can achieve higher compression ratio, better decoded image quality, and shorter encoding time than the previous similar methods.In our research, we propose an improved gray-level transformation whose fitting surface is used to fit both the range blocks and domain blocks.Moreover, the fitting surface in our method can approximate the range block better than Wang's fitting plane.Due to the no-search and quadtree method adopted, the proposed method can also maintain a fast encoding speed.Moreover, an initial image for fast fractal decoding is also proposed.Generally, accelerating the fractal decoding can be realized by two ways: one is to adopt improved iteration strategies [21][22][23], such as the one buffer decoding method [23].The other one is to select an initial image which can approximate the original image well.Moon et al. [24] proposed an approximated range-averaged image as the initial image and the rangeaveraged image is considered as the ideal initial image.In our research, we propose a novel initial image which can provide faster decoding speed.Since the proposed fitting surface itself is more similar to the corresponding range block, all the fitting surfaces can constitute a better fitting surface image (FSI) as the initial image in fractal decoding process without extra computations.Experimental results show that the proposed no-search fractal encoding method can provide better performance than Jacquin [2], Tong [25], and Wang's methods [20].Moreover, compared with the range-averaged image, FSI can also provide faster decoding speed.Finally, by combining our proposed fractal encoding method with JPEG, we design a hybrid coding method which can provide higher PSNR than JPEG while maintaining the same Bpp.
This paper is organized as follows.In Section 2, we describe the conventional fractal image coding method and the corresponding improved gray-level transforms.In Section 3, an improved gray-level transform with a second order fitting surface is proposed.In Section 4, a novel initial image is introduced to accelerate the fractal decoding process.In Section 5, the experimental results, such as Bpp, encoding time, and decoded image quality, are given and analyzed in detail.The final conclusion is presented in Section 6.

Review of Fractal Image Coding
According to the Collage theorem, for an input image F, the fractal image encoding process is to construct an iteration function system  which can be used to reconstruct the original image approximately.The reconstruction procedures can be described as follows: where  indicates the th iteration.F 0 is the initial image.F denotes the fixed point of  which is an approximated image of F.
In the practical implementation of the fractal encoding process, the input image is firstly partitioned into nonoverlapping range blocks.The domain blocks are obtained by sliding a square window throughout the original image.Commonly, the size of the domain blocks is twice the size of the range blocks.For each range block, the corresponding best matching domain block can be found in a domain block pool.The domain block pool is built by contracting all the domain blocks to the same size of the range blocks.Furthermore, the domain block pool is extended by eight isometric transformations which contain the identity, rotation through 90 ∘ , rotation through 180 ∘ , rotation through 270 ∘ , reflection about the middle vertical axis, reflection about the middle horizontal axis, reflection about the first diagonal, and reflection about the second diagonal.Figure 1 illustrates the above eight isometric transformations in order.In order to reduce the matching error, each contracted domain block should be transformed by the following gray-level transform: where  and  denote the scaling coefficient and offset coefficient, respectively.For each range block R, the best matching domain block can be obtained by minimizing the following equation: where ‖•‖ is the two-norm.I denotes a matrix whose elements are all ones.Error is the minimum matching error for each range block.With the least squares method,  and  can be computed as follows: where ⟨•, •⟩ is the inner product. and  are the mean values of the contracted domain block and the range block, respectively.The fractal encoding process can be finished by storing all the transformation parameters as the fractal code.At the fractal decoding phase, the transformation parameters are recursively applied to any starting image and an iteration image sequence is constructed.The image sequence converges after about ten iterations.By selecting the commonly used black image as the initial image, Figure 2 illustrates the first six iteration images for Barbara image.
Different from (2), Tong and Pi [25] proposed a modified gray-level transform as follows: The corresponding minimum matching error function can be described as where   is the scaling coefficient. and  are also the mean values of the contracted domain block and the range block, respectively.The main advantage of this method is that we only need to search the optimal coefficient   and it is more efficient than the previous method in (2).In order to reduce the minimum matching error further, Wang et al. [20] proposed an improved gray-level transform as follows: The corresponding minimum matching error function can be described as where P R and P D are the fitting planes of the range block R and the contracted domain block D, respectively.The fitting plane can be described as follows: where  is the width or height of the range blocks. and  are the variables to represent the position of each pixel. is the pixel gray value., , and  are the parameters to be estimated.

An Improved Gray-Level Transform Using a Fitting Surface
In this section, a second order surface instead of a first order one in ( 9) is proposed which can be represented as follows: where For arbitrary range block R, we need to find the coefficients , , , , and  that can make the surface best fit the range block R. We can solve this problem by the least squares method: where  is pixel gray value in (10).Set the derivatives of ( 12) with respect to , , , , and  to zero, respectively.By solving the system of equations, we can obtain the optimal coefficients as follows: Figure 3(a) illustrates an original surface which can be considered as a range block.Figures 3(b) and 3(c) are Wang's fitting plane and the proposed fitting surface, respectively.We can observe that the proposed fitting surface approximates the original surface better than the fitting plane.Substantially, the mean squared errors of Wang's fitting plane and the proposed fitting surface are 300 and 200, respectively.
Furthermore, we propose an improved gray-level transform by making use of the surface stated in (10) as follows: The corresponding minimum matching error function can be denoted as where S R is the fitting surface of the range block R and both the range block R and the domain block D adopt S R as their fitting surface.On one hand, we do not need to use another surface to represent the domain block D and one set of parameters can be saved to improve the compression ratio.On the other hand, we suppose that, for arbitrary range block R, the best matching domain block D is very similar to the range block R and can also provide minor matching error with the same fitting surface S R .In order to check the effectiveness of the proposed gray-level transform, in Figure 4(a), we select one hundred 8 × 8 range blocks randomly in Barbara image.Figure 4(b) illustrates their matching errors and we can see that, for most of the range blocks, our proposed method can provide smaller matching errors than Wang's method [20].In addition, the minimum matching error function in (15) can be also described as follows: From ( 16), we know that, since the proposed fitting surface S R can approximate the range block R better than the fitting plane, smaller matching errors can be achieved in our method.This is consistent with Figure 4(b).

A Novel Initial Image for Fast Decoding
At the decoding phase of our proposed method, in each level of the quadtree structure established in the encoding process, all the range blocks will be reconstructed by the same way as (14).One iteration is finished by combining all the reconstructed range blocks with different sizes.S R can be computed off line by (10) and all the fitting surfaces will be used repeatedly.Since ( 14) is invariable in the decoding process,  once the iteration strategy is specified, the convergence speed of the decoding process depends heavily on the initial image.The initial image is more similar to the input image; the reconstructed image in each iteration will approximate the decoded image better.In our research, instead of the initial black image which is commonly used, we proposed an initial image represented as follows: where S  R is the th fitting surface, ⋃(•) denotes the seamless splicing operation, and an initial image called fitting surface image (FSI) can be constituted by all the S R s with the same order as corresponding range blocks.Since all the S R s can be computed off line, we can get FSI without extra computations.By Moon's opinion [24], the range-averaged image is considered as the ideal initial image.The rangeaveraged image is combined with the blocks whose pixel intensities are all the mean values of the corresponding range blocks.Compared with the fitting surface we proposed, the blocks in range-averaged image can be only regarded as the simplest fitting surface.Thus, the proposed fitting surface in (10) can approximate arbitrary range block better than mean value block, and analogously the proposed FSI will be more similar to the original image than the range-averaged image.If the size of range blocks is set to be 8 × 8, Figures 5(a), 5(b), and 5(c) illustrate the original image, range-averaged image, and FSI for Barbara image, respectively.We can observe that, compared with the range-averaged image, FSI can preserve more details and approximate the original image better.

Experiments
The 512 × 512 Barbara image is used in our experiment.The scaling coefficient   is quantized by two bits and they are 0, 0.25, 0.5, and 0.75.The parameters , , , , and  for the fitting surface are quantized by 4, 4, 4, 4, and 6 bits, respectively.The experiments are carried out on a Pentium Dual-Core 2.93 GHz PC and programmed using MATLAB software.In order to estimate the performance of our proposed method, we will compare our algorithm with Jacquin [2], Tong [25], and Wang's methods [20] by the encoding time, the quality of the decoded image, and Bpp.The decoded image quality is measured by peak signal to noise ratio (PSNR): PSNR = 10 log 10 ( 255 2 where  and  are the height and width of the input image.f and f * represent the original image and the decoded image, respectively.

Compare the Proposed Method with the Other Three
Methods.In this part, the no-search and quadtree method will be incorporated into the fractal coding method.The quadtree structure we use contains four levels and the sizes of range blocks are 16 × 16, 8 × 8, 4 × 4, and 2 × 2, respectively.In the encoding process, the input image is firstly divided into nonoverlapping range blocks of 16 × 16 pixels.We encode the range blocks one by one with the fractal encoding method and compare the matching error with a given error threshold   ,  = 1, 2, 3, which changes by  +1 = 2 ×   as the quadtree level increases.Whether the range blocks should be divided into four subblocks is determined by whether the matching error is larger than the threshold.The above procedures will be repeated level by level until the size of range blocks is 2 × 2. We summarize the above procedures as follows: Step 1. Set an error threshold  1 and divide the input image into range blocks of 16 × 16 pixels.Assign extra two bits for each range block which are used to label the quadtree level it belongs to.
Step 2. For each range block, as illustrated in Figure 6, appoint the domain block as the best matching domain block directly which has the same center with the range block, compute the corresponding coefficients , , , , and  with (13), search the best scaling coefficient   , and obtain the minimum matching error by (15).If the minimum matching error is less than   , store the corresponding coefficients , , , , , and   as the fractal code and label them.If not, divide the range block into four subblocks uniformly and insert them into the next level.
Step 3. Repeat Step 2 until the size of range blocks is 2 × 2.
Encode them with the gray-level transform in (5) and store all the coefficients as the fractal code.We compare the proposed method with the other three methods under the same environment and the main differences are the gray-level transforms in (2), ( 5), (7), and (14).The encoding speed is mainly determined by the total number of blocks encoded in four quadtree levels.Smaller number of blocks matching can result in shorter encoding time.From the last column of Table 1, it is clearly seen that, for each  1 , the proposed method needs to encode the least number of range blocks and this results in the shortest encoding time which is shown in the last row of Table 2.The number of encoded blocks increases from the proposed method to Wang's method to Jacquin and Tong's methods.Jacquin and Tong's methods need to encode the most range blocks and this will result in longer encoding time which is also shown in the eleventh and twelfth rows of Table 2.We also can observe that Jacquin and Tong's methods perform almost the same number of blocks matching.Furthermore, in each block matching, the parameters in Jacquin's method are computed by ( 4) which leads to more computations than Tong's method.Thus, the encoding time of Jacquin's method is longer than Tong's method.
From Figures 3 and 4, we know that the fitting surface S R can approximate the corresponding range block R better; the corresponding matching error (R, D) can achieve a smaller value than Wang's method.Consequently, for a given threshold  1 , more range blocks will be encoded in the quadtree level with larger block size.From Table 1, we can see that, from Jacquin and Tong's methods to Wang's method to our method, the number of blocks with larger size increases and the number of blocks with smaller size decreases.The number of 16 × 16 blocks in our method is always larger than the other three methods.On the one hand, smaller matching errors in the encoding process can result in higher quality of the decoded image which is listed in the sixth row of Table 2. On the other hand, larger number of blocks with larger size will result in improving Bpp. Figure 7(a) also illustrates the PSNR versus Bpp for the four methods; we can see that the proposed method can obtain the highest decoded image quality all the time with different Bpps.If Bpp is larger than 1.2, the PSNR of Wang's method is slightly smaller than Tong's method.If not, Wang's method can obtain higher PSNR than Tong's method.Figure 7(b) illustrates the PSNR versus encoding time for the four methods.We can see that, under the same encoding time, our proposed method can provide the highest PSNRs than the other three methods.Figure 8 illustrates the decoded Barbara image with  1 = 11.Figure 9 illustrates the quadtree structure for Barbara image with  1 = 11.In summary, compared with the other three methods, the proposed method can provide higher PSNR at the same Bpp and higher PSNR at the same encoding time.

Fast Fractal Decoding Based on FSI.
In the following part, we will estimate the performance of FSI.In order to describe the convergence speed conveniently, we define the following RMSE to describe the difference between each iteration and the final decoded image: where f  and f Decoded are the th iteration image and decoded image, respectively.Based on the fractal code, with a number of iterations, we first get the final decoded images as the reference image, and then we select the range-averaged image and FSI as the initial image, respectively.The same reconstruction method, described by (14), is adopted for both above two initial images.Figure 10(a) illustrates the RMSEs of each iteration with respect to the decoded image for Barbara image.Figure 10(b) illustrates the RMSEs with respect to the decoding time.We can see that FSI can provide smaller RSME all the time and achieve faster fractal decoding speed than the range-averaged image.

Hybrid Coding Method Based on the Proposed Fractal
Coding Method.In this part, we first compare the proposed no-search fractal coding method with JPEG by compression ratio (CR) and PSNR.The error threshold  1 in fractal coding is set to be 11 and the experimental results are listed in Table 3.We can see that, with almost the same decoded image quality, JPEG can provide higher CR than the fractal coding method for the whole image, but, for the blocks of 16 × 16 pixels, the fractal coding method can provide higher CR while maintaining almost the same decoded image quality as JPEG.Thus, we suppose that if the blocks of 16 × 16 pixels are encoded by the proposed fractal algorithm and the remaining blocks are encoded by JPEG, better performance can be obtained.The encoding process of the above hybrid coding method can be described as follows.
Step 1. Assign extra one bit to label the range block encoded by the fractal coding method.
Step 2. Partition the input image into range blocks of 16 × 16 pixels and encode them one by one.For arbitrary range block, if the collage error is smaller than  1 , encode the range block by the proposed fractal encoding method and label it; if not, partition the range block into four subblocks of 8 × 8 pixels.
Step 3. Encode all the remaining blocks of 8 × 8 pixels by JPEG encoder.
The corresponding decoding process can be described as follows.
Step 1.According to marking vector, decode the blocks of 16 × 16 pixels and 8 × 8 pixels with the fractal decoder and the JPEG decoder, respectively.
Step 2. Repeat Step 1 about ten times and the decoding process will converge to the final decoded image.
From the above analysis we know that, for the range blocks of 16 × 16 pixels, the fractal coding method can obtain    good decoded image quality and high CR, simultaneously, but JPEG can only obtain either good decoded image quality or high CR.Thus, the decoded image quality can be improved by the local image encoded by the fractal coding method.Figure 11 illustrates the PSNR versus Bpp for the hybrid method and JPEG.It can be clearly seen that, with the same Bpp, the hybrid method can provide better PSNR than JPEG.

Conclusions
In order to reduce the matching error between the range block and its matching domain block, an improved gray-level  transform with a fitting surface is proposed in this paper.Firstly, the proposed fitting surface can approximate the range block better.Consequently, higher quality of the decoded image can be achieved with respect to the other three similar methods.Moreover, we adopt the same fitting surface for the range block and its matching domain block.This can save many bits and result in improving Bpp.In addition, due to smaller matching errors, more range blocks with larger size can be encoded and this can reduce the total number of blocks matching; that is, shorten the encoding time.Experimental results show that the proposed method can yield superior performance over the other three similar methods.Moreover, FSI can also accelerate the fractal decoding process.Finally, by combining the proposed fractal encoding method with JPEG, PSNR can be achieved higher than JPEG while maintaining the same Bpp.

Figure 2 :
Figure 2: The first six iterations for Barbara image in fractal decoding.

Figure 3 :Figure 4 :
Figure 3: Comparison between the fitting plane and the fitting surface.(a) Original surface.(b) Fitting plane.(c) Fitting surface.

Figure 5 :
Figure 5: Comparison between Moon's range-averaged image and FSI with respect to the original image.(a) Original image.(b) Rangeaveraged image.(c) FSI.

Figure 6 :
Figure 6: Relative position of the range block and its corresponding domain block.

Figure 7 :
Figure 7: Comparison between the proposed method and the other three methods for Barbara image.(a) PSNR versus Bpp.(b) PSNR versus encoding time.

Figure 10 :
Figure 10: Comparison of the fractal decoding for the range-averaged image and FSI.(a) RMSE versus iterations.(b) RMSE versus decoding time.

Figure 11 :
Figure 11: PSNR versus Bpp for the hybrid coding method and JPEG.

Table 1 :
Number of blocks in different quadtree levels for the four methods.

Table 2 :
Performance comparison between the proposed method and the other three methods.

Table 3 :
Performance comparison between the proposed fractal coding method and JPEG while  1 = 11.