- Email: [email protected]
Genetic algorithm with a hybrid select mechanism for fractal image compression
Digital Signal Processing 20 (2010) 1150–1161
Contents lists available at ScienceDirect
Digital Signal Processing www.elsevier.com/locate/dsp
Genetic algorithm with a hybrid select mechanism for fractal image compression Ming-Sheng Wu a , Yih-Lon Lin b,∗ a b
Department of Electrical Engineering, Cheng Shiu University, Kaohsiung, Taiwan Department of Information Engineering, I-Shou University, Kaohsiung, Taiwan
a r t i c l e
i n f o
a b s t r a c t
Article history: Available online 4 January 2010
In this paper, a genetic algorithm with a hybrid select mechanism is proposed to speed up the fractal encoder. First, all of the image blocks including domain blocks and range blocks are classified into three classes: smooth; horizontal/vertical edge; and diagonal/subdiagonal edge, according to their discrete cosine transformation (DCT) coefficients. Then, during the GA evolution, the population of every generation is separated into two clans: a superior clan and an inferior clan, according to whether the chromosome type is the same as that of the range block to be encoded or not. The hybrid select mechanism proposed by us is used to select appropriate parents from the two clans in order to reduce the number of MSE computations and maintain the retrieved image quality. Experimental results show that, since the number of MSE computations in the proposed GA method is about half of the traditional GA method, the encoding time for the proposed GA method is less than that of the traditional GA method. For retrieved image quality, the proposed GA method is almost the same as the traditional GA method or only has a little decay. Moreover, in comparison with the full search method, the encoding speed of the proposed GA method is some 130 times faster than that of the full search method, whereas the retrieved image quality is still relatively acceptable. © 2009 Elsevier Inc. All rights reserved.
Keywords: Fractal image compression Discrete cosine transformation Classification method Genetic algorithm Hybrid select mechanism
1. Introduction Fractal image compression is potentially a great coding scheme since it features a high compression ratio and good retrieved image quality. The original idea was proposed by Barnsley [1–3] according to the iterated function system (IFS). In 1992, the first fractal image coding algorithm was realized by Jacquin [4]. The underlying premise of the coding scheme is based on the partitioned iteration function system (PIFS) which utilizes the self-similarity characteristic in a natural image to achieve compression. In addition to image compression, the application of fractals is extremely widespread, in such things as image restoration [5], medical image classification [6], character recognition [7], and watermarking [8,9]. It is therefore little wonder this technique attracts the attention of many researchers. The drawback to fractal image compression is that the encoding process is very time-consuming since, for each range block, a large number of computations of similarity measures must be done to find the best matched domain block. Hence the main research topic of the field is focused on how to reduce the encoding time. In the past, the techniques of partition and classification were used to reduce the number of computations of similarity measures to speedup the encoder [10– 12]. Then the idea of spatial correlation was also added to the fractal encoding algorithm to reduce the search space of the optimal solution. Truong et al. [13] limit the search space for the current range block on the neighborhoods of the
*
Corresponding author. Address for correspondence: No. 1, Sec. 1, Syuecheng Rd., Dashu Township, Kaohsiung County 840, Taiwan. Fax: +886 7 6578944. E-mail addresses: [email protected], [email protected] (Y.-L. Lin).
1051-2004/$ – see front matter doi:10.1016/j.dsp.2009.12.009
©
2009 Elsevier Inc. All rights reserved.
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
1151
matched domain blocks of the neighboring range blocks by utilizing the spatial correlations between neighboring blocks in both the domain pool and the range pool. In comparison with the full search algorithm, this algorithm achieves a 2.6 times speedup ratio yet preserves almost the same retrieved image quality. Furao and Hasegawa [14] present a no search fractal encoding algorithm. In this algorithm, they select a domain block which has the same center with the range block to be encoded as the matched one. This is possible since the two blocks usually have the same texture in a natural image. Utilizing such a characteristic and merging quad-tree technique, the authors have achieved the fastest encoding velocity to date, but unfortunately both the retrieved image quality and the compression ratio suffer as a result. Recent efforts to push the genetic algorithm (GA) to a numerical optimum solution have become highly popular among researchers. GA was developed by Holland (1975) [15] over the course of the 1960s and 1970s and popularized by Goldberg [16]. The genetic algorithm is a global search technique that mimics natural selection and natural genetics. GA is capable of solving many large, complex problems. Experience shows such problems are hard to solve, since their search spaces generally have very rough landscapes riddled with many local optima. Many researchers have used the advantages of GA to achieve some significant results [17–21]. Since the landscape of natural images is also very rough, GA is well suited to the search for the best match in fractal image compression. For example, Mitra [22] proposed a search strategy using GA in which the chromosome consists of the PIFS of all range blocks. In this method, the fitness of a chromosome is measured by the distance between the attractor and the given image. Although a faster speed can be achieved, the retrieved image quality is not good enough because the fitness is defined in terms of the whole image. Wu et al. [23] used a schema genetic algorithm (SGA) to speed up the fractal encoder. In SGA, the genetic operators are adapted according to the schema theorem in the evolutionary process performed on the range blocks. The encoding speed of such a method is about 100 times faster than that of the full search method and can preserve the retrieved image quality. In this paper, a genetic algorithm with a hybrid select mechanism is proposed to speed up fractal encoding. All of the image blocks are classified into three classes: smooth type; diagonal/sub-diagonal edge type; and horizontal/vertical edge type. Two DCT coefficients – the lowest vertical coefficient and the lowest horizontal coefficient – are used to determine which type the block belongs to. Then, for each range block to be encoded, GA is used to find the best matched solution. During the GA evolutionary process, all the chromosomes are classified as one of two clans: a superior clan or an inferior clan. In the superior clan, the type of chromosome is the same as that of the range block. In the inferior clan, it is not. If a chromosome belongs to the superior clan, it is a good candidate for the best match. Hence its fitness value will be calculated. The chromosomes belong to the superior clan are ranked according to their fitness values and placed at the front of the chromosomes belong to the inferior clan. In contrast, if a chromosome belongs to the inferior clan, it is not a good candidate for the best match. Hence its fitness value need not be calculated. The chromosomes belong to inferior clan will be assigned the same rank. Finally, two selection mechanisms – uniform random selection and ranking selection – are merged in the GA mating process to select appropriate parents in order to find a suitable domain block and speedup the encoder. By using the evolutionary technique, we find that the encoding speed can be improved, since the number of fitness computations is reduced. Moreover, since good chromosomes are ranked at the front, good genes can be carried into the next generation during the GA evolutionary process. In this way, the quality of the retrieved image can also be maintained. The rest of the paper is organized as follows. We introduce the conventional fractal image coding scheme in Section 2. Section 3 describes the proposed genetic algorithm in this paper. Some experimental simulations are performed in Section 4 to verify the improvement of our proposed algorithm. Finally, a conclusion is made in Section 5. 2. Fractal image encoding Fractal image compression is based on the local self-similarity property and PIFS. The related implementation process is stated as follows. Suppose the original gray level image f is of size m × m. To simplify the discussion, let m be 256. The domain pool is defined as the set of all possible blocks of size 16 × 16 of the image f , which makes up (256 − 16 + 1) × (256 − 16 + 1) = 58,081 blocks. The range pool is defined to be the set of all non-overlapping blocks of size 8 × 8 of the image f , which makes up (256/8) × (256/8) = 1024 blocks. For each block v from the range pool, the fractal transformation is constructed by searching all of the elements in the domain pool to find the block with the greatest similarity. Let u denote a subsampled domain block which is of the same size as v. The similarity of u and v is measured using mean square error (MSE), which is defined by
MSE =
7 7 1
64
2
u (i , j ) − v (i , j )
(1)
j =0 i =0
The fractal transformation allows Dihedral transformations of the domain blocks, i.e., the 8 orientations of the blocks generated by rotating the blocks counterclockwise at angles 0◦ , 90◦ , 180◦ and 270◦ and flipping with respect to the line y = x, respectively, as shown in Table 1. Thus for a given block from the range pool, there are 58,081 × 8 = 464,648 MSE computations to obtain the block with the greatest similarity from the domain pool. Thus, in total, one needs 1024 × 464,648 = 475,799,552 MSE computations to encode the whole image using this full search compression method.
1152
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
Table 1 The 8 transformations in the Dihedral group. Rotate 0◦
1 0
0 1
0 1
1 0
Rotate 90◦
0 1 −1 0
Rotate 180◦
−1 0
0 −1
Flip with the line x = y from above 1 0 0 −1
0
−1
−1 0
Rotate 270◦
0 −1 1 0
−1 0 0
1
The fractal transformation also allows the adjustment of the contrast p and the brightness q on the sub-sample block u. Thus the similarity is to minimize the quantity d = p · uk + q − v , where uk , 0 k 7 are the 8 orientations of u. By calculus, p and q can be computed directly by
p=
[ N u , v − u , 1 v , 1 ] , [ N u , u − u , 1 2 ]
q=
1 N
v , 1 − p u , 1
= [1 1 . . . 1] T . where N is the number of pixels of the range block and 1 Finally, the position x and y of the domain block, the contrast p, the brightness q, and the orientation k constitute the fractal code of the given range block v. For a 256 × 256 image, both x and y require 8 and 8 bits, respectively, to represent the position of the domain block. For contrast p, brightness q, and the orientation k, 5, 7 and 3 bits are required, respectively. Hence one needs 8 + 8 + 5 + 7 + 3 = 31 bits in total to encode a range block. Finally, as v runs over all 1024 range blocks in the range pool, the encoding process is completed. To decode, one first makes up the 1024 affine transformations from the compression codes. Next, one chooses any initial image and performs the 1024 affine transformations on the image to obtain a new image. The transformation proceeds recursively according to the Contractive Mapping Fixed-Point Theorem and Collage Theorem until the sequence of images converges to the encoded image. The stopping criterion of the recursion is designed in accordance with the user’s application and the final image is the retrieved image of fractal coding. 3. Genetic algorithm for fractal image compression In this section, a genetic algorithm with hybrid select mechanism proposed by us will be introduced. First, by using two DCT coefficients: the lowest vertical coefficient F (1, 0) and the lowest horizontal coefficient F (0, 1), all of the image blocks including domain blocks and range blocks are classified into three classes: smooth type, diagonal/sub-diagonal edge type, and horizontal/vertical edge type according to their texture. Then when the GA is used to encode each range block, a hybrid select mechanism composed of two different kinds of selection mechanisms are used to select the appropriate parents to speed up the encoder. Below, the classification method for the image blocks is introduced first. Second, the proposed GA algorithm designed to reduce encoding time will also be stated in detail. 3.1. The classification algorithm for image blocks The main goal of the classification algorithm is to classify all of the image blocks including domain blocks and range blocks into three classes. By using DCT transformation, an image block can be transferred from the spatial domain to the frequency domain, in which the DCT coefficients located in the upper-left represent the low frequency information of the image block and reflect the rough contour of the image block. In contrast, the DCT coefficients located in the lower-right represent the high frequency information of the image block and reflect the fine texture of the image block. Hence an image block can be determined as belonging to which type, only by using its lower frequency DCT coefficients. Based on the viewpoint stated above, we design a classifier which uses the lowest vertical coefficient F (1, 0) and the lowest horizontal coefficient F (0, 1) to execute the action of classification. Assume f is a given image block of size 8 × 8, F (1, 0) and F (0, 1) can be computed from
√ F (1, 0) =
2
8
√ F (0, 1) =
7
7
2
8
f (i , j ) cos θi
(2)
f (i , j ) cos θ j
(3)
j =0 i =0 7
7
j =0 i =0
where θi = (2i + 1)π /16, i = 0, 1, . . . , 7. The magnitude of F (1, 0) reflects the intensity variation between the left half and right half of image block f and the magnitude of F (0, 1) reflects the intensity variation between the upper half and lower half of image block f . For an image block f , if both | F (1, 0)| and | F (0, 1)| are small, it means that the block f has no clear texture. Hence the landscape on
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
1153
Fig. 1. The classification of H type blocks and non-H type blocks.
block f tends to smoothness. Block f will be classified as smooth type and denoted as S type. Second, if | F (1, 0)| is large and | F (0, 1)| is small, f will have a vertical edge. But if | F (1, 0)| is small and | F (0, 1)| is large, f will have a horizontal edge. According to Dihedral characteristics, for the two cases, f is classified as horizontal/vertical edge type and denoted as H type. Finally, if both | F (1, 0)| and | F (0, 1)| are large, it means that f will have a diagonal or sub-diagonal edge. Such f will be classified as diagonal/sub-diagonal type and denoted as D type. To determine the different types of image blocks, these blocks first execute the DCT transformation to find their F (1, 0) and F (0, 1) values. Then the classification algorithm is used to determine which type these image blocks belong to, according to their respective F (1, 0) and F (0, 1) values. First, H type blocks are separated out. For each image block, we take the absolute values of its F (1, 0) and F (0, 1) to constitute a DCT absolute coefficient pair (| F (1, 0)|, | F (0, 1)|). Then the DCT absolute coefficient pair is undergoes the normalization process. After normalization, the normalized DCT absolute coefficient pair (| F (1, 0)| N , | F (0, 1)| N ) is obtained and is labeled on a Cartesian coordinate system in which | F (1, 0)| N is the x-coordinate and | F (0, 1)| N is the y-coordinate. The result is shown in Fig. 1. Here | F (1, 0)| N and | F (0, 1)| N denote normalized | F (1, 0)| and| F (0, 1)|, respectively. We can find from Fig. 1 that the positions of all the normalized image blocks are located on the unit circle of the first quadrant. According to the texture characteristics of the image block, if the normalized DCT absolute coefficient pair of an image block f is near the point (1, 0), it will have a clear vertical edge since | F (1, 0)| N is large and | F (0, 1)| N is small. Hence the image block f will be classified as H type. For the same reason, if the normalized DCT absolute coefficient pair of an image block f is near the point (0, 1), it will have a clear horizontal edge since | F (1, 0)| N is small and | F (0, 1)| N is large. Such an image block f will also be classified as H type. Two threshold ◦ and T ◦ are used to separate H type blocks from non-H type blocks. Therefore, an image block is classified angles: T H1 H2 ◦ or between T ◦ and 90◦ . Once the H as H type if its position in Fig. 1 is located on the unit circle between 0◦ and T H1 H2 type blocks are determined, the S type blocks are also isolated from the remaining non-H type blocks. For the remaining image blocks, according to their DCT absolute coefficient pairs (| F (1, 0)|, | F (0, 1)|) which have not been normalized, we find an image block whose distance from the origin of the coordinate is greater than the others and label its coefficient pair as (| F (1, 0)|max , | F (0, 1)|max ). Then the DCT absolute coefficient pair (| F (1, 0)|, | F (0, 1)|) of the remaining image blocks is divided by ((| F (1, 0)|max )2 + (| F (0, 1)|max )2 )1/2 . Through the scaling process, the scaled DCT absolute coefficient pair (| F (1, 0)| R , | F (0, 1)| R ) is obtained and is re-labeled on a Cartesian coordinate system in which | F (1, 0)| R is the x-coordinate and | F (0, 1)| R is the y-coordinate. The result is shown in Fig. 2. In Fig. 2, (| F (1, 0)|max, R , | F (0, 1)|max, R ) is the image block whose distance from the origin of the coordinates is greatest. We find that, through the scaling process, the positions of the ◦ and T ◦ . If the position of an image block f is near rest of the image blocks are located in the unit circle between T H1 H2 the origin of the coordinate, it will have a smooth landscape since both | F (1, 0)| R and | F (0, 1)| R are small. Such an image block f will be classified as S type. A threshold circle whose radius is T S is used to separate the S type and non-S type blocks. Hence an image block f is classified as S type if its position is located in the threshold circle. Finally, the non-S type blocks will be classified as D type. For these non-S type blocks, they have a clear diagonal or sub-diagonal edge since their | F (1, 0)| R and | F (0, 1)| R values are large. The classification steps for the image blocks are discussed in detail below: 1. For all the image blocks, calculate the F (1, 0) and F (0, 1) from (2) and (3), respectively, and constitute the DCT absolute coefficient pair (| F (1, 0)|, | F (0, 1)|). 2. Normalize the DCT absolute coefficient pair (| F (1, 0)|, | F (0, 1)|) for all the image blocks to obtain the normalized DCT absolute coefficient pair(| F (1, 0)| N , | F (0, 1)| N ). Label the position of each image block according to its (| F (1, 0)| N , | F (0, 1)| N ) on the Cartesian coordinate system. The result is shown in Fig. 1.
1154
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
Fig. 2. The classification of S type blocks and D type blocks.
Fig. 3. Chromosome formation.
3. Determine the H type blocks from Fig. 1. An image block is classified as H type if its position is located on the unit ◦ or between T ◦ and 90◦ . circle between 0◦ and T H1 H2 4. For the remainder in Step 3, find the pair (| F (1, 0)|max , | F (0, 1)|max ) according to their DCT absolute coefficient pair (| F (1, 0)|, | F (0, 1)|) obtained in Step 1. Divide the (| F (1, 0)|, | F (0, 1)|) of the remaining blocks by ((| F (1, 0)|max )2 + (| F (0, 1)|max )2 )1/2 to obtain the scaled DCT absolute coefficient pair (| F (1, 0)| R , | F (0, 1)| R ). Label their position on the Cartesian coordinate system. The result is shown in Fig. 2. 5. Determine the S type blocks from Fig. 2. An image block is classified as S type if its position is located in a circle whose radius is T S . 6. Classify the rest of the blocks in Step 5 as D type. Once the classification process is complete, for each range block, a genetic algorithm with a hybrid select mechanism proposed by us will be used on the classified image blocks in order to search for the best matched domain block. 3.2. A genetic algorithm with a hybrid select mechanism for FIC Here, a genetic algorithm with hybrid select mechanism is introduced in detail. For each range block to be encoded, the GA algorithm executes GA evolution on the classified domain blocks obtained in Section 3.1 to find the best matched solution. Its main goal is to reduce the amount of computation of similarity measure between range block and sub-sample domain block in order to speed up the encoding speed while maintaining the retrieved image quality. The related evolutionary strategies with respect to the GA are detailed as follows. (1) Chromosome: Since the fractal encoding scheme utilizes the PIFS to encode every range block, one takes the absolute position (x, y ) of a domain block and the Dihedral transformation k to constitute a chromosome. A chromosome is 19 bits in length as shown in Fig. 3, in which x, y and k are the genes of 8, 8 and 3 bits, respectively. (2) Fitness function: The distance between range block and sub-sampled domain block is measured by MSE. The fitness value is defined as the reciprocal of MSE. (3) Initial population: Chromosomes are initialized randomly. (4) Selection: Here, a hybrid selection mechanism which is composed of uniform random selection and ranking selection is adopted to find appropriate parents and reduce the amount of MSE computation. All of the chromosomes are divided into two clans:
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
1155
Fig. 4. Uniform crossover.
a superior clan and an inferior clan. If the type of chromosome is the same as that of the range block to be encoded, it is a good candidate for the best match. It is proof that the gene in the chromosome is good. Hence it should have a higher probability to be selected as a parent in order to generate better offspring. Such a chromosome will be classified as belonging to the superior clan. In contrast, if the type of chromosome is not the same as that of the range block to be encoded, it is not a good candidate for the best match. Its gene is not good and hence the probability it is selected as a parent should be lower. Such a chromosome will be classified as belonging to the inferior clan. For those chromosomes belonging to the superior clan, their genes are good. We hope that these genes can be preserved to next generation to generate better offspring. Hence we let their fitness be calculated and be ranked from high to low according to their fitness values. But for those chromosomes belonging to the inferior clan, their genes are not good. These genes will be eliminated during the GA evolution. Therefore, for reducing the number of MSE computations, their fitness will not be calculated. We assign them the same rank and rank them directly behind the superior clan. For these chromosomes, their rank is equal to the number of chromosomes in the superior clan plus 1. Once all the chromosomes have been ranked, the hybrid select mechanism will be executed in order to select appropriate parents from the ranked chromosomes. First, the ranking selection is carried out to select parents from all the ranked chromosomes. If a chromosome belonging to the superior clan is selected, it will be immediately taken as the parent. But if the select mechanism determines the parent as coming from the inferior clan, the uniform random selection will be invoked in order to select a parent from the inferior clan. For example, suppose the size of the population is 7 and the number of chromosomes belonging to the superior clan is 3. First, the ranking select mechanism is used to find a parent from the population. For the 3 chromosomes in the superior clan, the selected probabilities according to their rank are 7/28, 6/28, and 5/28, in that order. The probability that the parent comes from the inferior clan is 4/28 + 3/28 + 2/28 + 1/28 = 10/28. By use of the random generator, if the random number lies between 0 and 7/28, the chromosome whose rank is 1 will be selected. For the same reason, if the random number lies between 7/28 and 13/28 or between 13/28 and 18/28, the chromosome whose rank is 2 or 3 will be selected, respectively. But if the random number lies between 18/28 and 1, the ranking select mechanism determines that the parent comes from the inferior clan. Hence the uniformly random selection will be invoked to select a parent from the 4 chromosomes in the inferior clan. On average, the selected probability for each chromosome in the inferior clan is (10/28)/4 = 5/56. The fact the parent comes from the inferior clan maintains the diversity of the population. (5) Crossover: The uniform crossover operation, as shown in Fig. 4, is adopted with probability P c to determine whether the operation will be performed or not. The temporary offspring are generated by executing the crossover operation on the two parents. (6) Mutation: The mutation operation with probability P m is applied to the temporary offspring in order to generate the offspring for the next generation. Here the goal of the mutation is to maintain the diversity of the population and to avoid pre-maturity. (7) Stopping criterion: The goal of stopping criterion is to stop the GA evolutionary process. Here, a pre-set number of iterations is adopted to stop the evolution. Connecting the above operational strategies, a genetic algorithm with a hybrid select mechanism is proposed. The detailed evolutionary steps regarding the GA are given as follows: Initially, set the population size, the crossover mask, crossover rate P c , mutation mask, and mutation rate P m . 1. Generate the initial population of chromosomes randomly. 2. Divide the chromosomes into two clans: superior clan and inferior clan according to their types. Calculate the fitness values of the chromosomes belonging to the superior clan. 3. Rank all the chromosomes. First, for the chromosomes belonging to the superior clan, rank them according to their fitness values. Second, let the chromosomes belonging to the inferior clan have the same rank, and their rank is the number of chromosomes in the superior clan plus 1.
1156
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
Table 2 ◦ and T ◦ . The rank condition in part for the determination of T H1 H2 Rank
min( D 1,0 , D 0,1 )
Phase angle
337 338 339 340 341 342 343 344 345 346
0.1921 0.1922 0.1925 0.1934 0.1968 0.1970 0.1980 0.1995 0.2000 0.2008
78.98◦ 11.03◦ 11.05◦ 78.90◦ 11.29◦ 11.30◦ 11.36◦ 11.45◦ 11.48◦ 11.53◦
4. If a pre-set number of iterations is reached, then stop and record the fractal code. Otherwise, go to the next step. 5. Execute the hybrid select mechanism to select parents. First execute the ranking selection. If the selected chromosome belongs to the superior clan, take it immediately as the parent. In contrast, if the selected chromosome belongs to the inferior clan, invoke the uniformly random selection to select a parent from the inferior clan. 6. Perform the uniform crossover to generate a temporary offspring. 7. Perform the mutation operation on the temporary offspring to generate the offspring for the next generation. 8. Go to Step 2. 4. Experimental results The images Lena, Baboon, F16, and Pepper are tested to demonstrate the speed-up rate and retrieved quality of the proposed GA method. For a given image of size 256 × 256, let the size of range block be 8 × 8. The software simulations are performed using BCB v. 6.0 on a Celeron Core2 Duo 2.0 GHz, 2 G RAM, windows XP PC. The quality of the retrieved image g in comparison to the original image f is measured by peak signal to noise ratio (PSNR) given by
PSNR = 10 × log
2552
1 m×n
n−1 m−1 j =0
i =0
( f (i , j ) − g (i , j ))2
where m × n is the image size. 4.1. The determination of threshold values for the classification method All the image blocks will first be classified before the GA evolutionary process is executed. Hence three threshold values: ◦ , T ◦ and T must be obtained in order to execute the classification algorithm. It is obvious that the number of domain T H1 S H2 blocks is much larger than that of range blocks. If we use the domain pool to determine the threshold values, it will produce a lot of computation overhead. Thus, to reduce the computation load, we use the range pool to determine the threshold values instead of the domain pool and partition the range pool into three classes of almost equal numbers of members. First, ◦ and T ◦ in Fig. 1 are obtained to separate H type and non-H type blocks from all the image two threshold angles: T H1 H2 blocks. For each normalized range block, since the more that are near point (1, 0) or point (0, 1), the clearer the vertical or horizontal texture is, we calculate the value min( D 1,0 , D 0,1 ) for all the range blocks and rank them according to value from small to large. Here D 1,0 = ((| F (1, 0)| N − 1)2 + (| F (0, 1)| N − 0)2 )1/2 represents the distance between the normalized range block and point (1, 0) and D 0,1 = ((| F (1, 0)| N − 0)2 + (| F (0, 1)| N − 1)2 )1/2 represents the distance between the normalized range block and point (0, 1). The min( D 1,0 , D 0,1 ) represents the minimum value of both D 1,0 and D 0,1 . Table 2 shows the rank condition in part and the test image is Lena. In the last row of Table 2, the rank of the range block is 346. Its min( D 1,0 , D 0,1 ) is 0.2008 and the phase angle of the line connecting the range block and origin is 11.53◦ . For an image of size 256 × 256, with the size of range block 8 × 8, the number of range blocks is (256/8) × (256/8) = 1024. We let the first ◦ and T ◦ in Fig. 1, T ◦ is assigned to the third of range blocks, about 342, be classified as H type. For the assignment of T H1 H2 H1 middle point between the last rank of the H type blocks and the 1st rank of non-H type blocks. Furthermore, we let both ◦ ◦ . Hence from ◦ the L 1 to x-coordinate and L 2 to y-coordinate have the same span. In other words, we keep T H2 = 90 − T H1 ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ Table 2, we obtain the result T H1 is (11.30 + 11.36 )/2 = 11.33 and T H2 is 90 − 11.33 = 78.67 . If we let the range ◦ is ((90◦ − 78.90◦ ) + 11.29◦ )/2 = 11.20◦ and T ◦ is 90◦ − 11.20◦ = 78.80◦ , blocks of the first 340 be classified as H type, T H1 H2 since the 340th range block is near the y-coordinate. Second, the threshold radius T S is determined from the remaining 1024 − 342 = 682 range blocks. For the same reason, for the remaining range blocks, since the closer to the origin (0, 0), the smoother the characteristics, we calculate the value (| F (1, 0)|2R + | F (0, 1)|2R )1/2 for all the remaining range blocks and rank them according to their value from small to large. Table 3 shows the rank condition in part. In the last row of Table 3, the rank of the range block is 346, from the remaining 682, and the distance from the origin is 0.0822. We take the first half of the range blocks as S type and take the second half of the range blocks as D type. Hence both the number of S type blocks and the number of D type blocks are 682/2 = 341. For the same reason, T S is assigned to the middle point between the last
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
1157
Table 3 The rank condition in part for the determination of T S . Rank
Distance with origin
337 338 339 340 341 342 343 344 345 346
0.0781 0.0793 0.0795 0.0797 0.0803 0.0805 0.0808 0.0813 0.0821 0.0822
Table 4 Three threshold values and the classification result of the domain blocks by using the proposed classification method.
Lena Pepper F16 Baboon
◦ T H1
◦ T H2
TS
The number of H type blocks
The number of S type blocks
The number of D type blocks
11.33◦ 12.49◦ 12.23◦ 13.01◦
78.67◦ 77.51◦ 77.77◦ 76.99◦
0.0804 0.0986 0.0459 0.1504
19,223 19,276 18,782 20,350
11,852 12,818 14,924 14,759
27,006 25,987 24,375 22,972
Fig. 5. The number of chromosomes in the superior clan and inferior clan versus the number of generations for Lena – without elitism.
rank of the S type blocks and the 1st rank of the D type blocks. Therefore T S is (0.0803 + 0.0805)/2 = 0.0804. Table 4 lists the three threshold values and the classification result of the domain blocks by using the proposed classification method for ◦ , T ◦ , and T are 11.33◦ , 78.67◦ , and 0.0804, respectively, and the number of Lena, Pepper, F16, and Baboon. For Lena, T H1 S H2 domain blocks of H type, S type, and D type are 19,223, 11,852, and 27,006, respectively. 4.2. Comparison of proposed GA method and traditional GA method – without elitism Some simulations are done to demonstrate the speed-up ratio and retrieved quality of the proposed GA method in comparison to the traditional GA method. For GA operations, the proposed GA method uses a hybrid select mechanism proposed by us and the traditional GA method uses a ranking select mechanism. The uniform crossover is adapted for both GA methods. Elitism is ignored. At the same time, the GA parameters of the two methods are identical, in which P c is 0.6 and P m is 0.01. Fig. 5 shows the relationship between the number of chromosomes in the superior clan and the number of chromosomes in the inferior clan versus the number of generations for the proposed GA method. The test image is Lena and the size of the population is 250. The two curves in Fig. 5 represent the number of chromosomes in the superior clan and the inferior clan in every generation, respectively, in which the values taken are an average of 1024 range blocks in an image. At the beginning, the number of chromosomes in the superior clan is less than the number of chromosomes in the inferior clan. The number of initial chromosomes for the superior clan and inferior clan are 83.8 and 166.2, respectively, on average. During the GA evolution, the number of chromosomes from the superior clan increases gradually and the number of chromosomes from the inferior clan decreases gradually. At the 8th generation, the number of chromosomes in the superior clan exceeds
1158
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
Table 5 The comparison of proposed GA method and tradition GA method together with full search method – without elitism. Method
Time (sec)
MSE computations
PSNR (dB)
Lena
Full search Proposed GA Traditional GA
2589.53 35.48 47.64
475,799,552 2,831,391 5,120,000
28.91 26.25 26.23
Pepper
Full search Proposed GA Traditional GA
2587.91 32.68 47.93
475,799,552 2,361,127 5,120,000
29.84 26.65 27.10
F16
Full search Proposed GA Traditional GA
2582.92 34.25 47.70
475,799,552 2,653,875 5,120,000
25.24 23.01 23.32
Baboon
Full search Proposed GA Traditional GA
2591.36 33.45 48.32
475,799,552 2,501,905 5,120,000
20.15 19.34 19.47
the number of chromosomes in the inferior clan. Finally, at the 20th generation, the number of chromosomes from the superior clan and inferior clan have changed to 182.9 and 67.1, respectively, on average. We find from Fig. 5 that since the fitness for chromosomes in the inferior clan need not be calculated, we can save about two-thirds of the MSE computations at the GA evolutionary initial stage. Even at a later stage of the GA evolutionary process, we can still save about one-quarter of the MSE computations. On average, in comparison with the traditional GA method, about 45% of MSE computations can be saved. Table 5 lists some comparative results of our proposed GA method, the traditional GA method and the full search method for Lena, Pepper, F16, and Baboon. The size of the population is 250 and the pre-set number of iterations is 20. For a proper comparison, the two GA methods are performed over 10 runs and the results are averaged. For Lena, comparing our proposed GA method with the traditional GA method, the PSNRs of the two methods are almost the same, but for the proposed GA method, the number of MSE computations is only 55% that of the traditional GA method. In terms of the speed up rate, the encoding time of the proposed GA method can save about 25% off that of the traditional GA method, in which the encoding time of our method already includes the time spent on the DCT calculation and classification of all the image blocks. Moreover, for Pepper, F16, and Baboon, the encoding time of the proposed GA method saves about 32%, 28%, and 31%, respectively, from that of the traditional GA method, and the penalty is only 0.13–0.45 dB decay on retrieved image quality. Finally, in comparison with the full search method, the encoding speed of the proposed GA method is about 75 times faster than that of the full search method. 4.3. Comparison of proposed GA method and traditional GA method – with elitism In this section, elitism is considered. The relevant GA operators and parameters are the same as in Section 4.2. Fig. 6 shows both the number of chromosomes in the superior clan and the number of chromosomes in the inferior clan versus the number of generations for the proposed GA method. The test image is also Lena and the size of the population is 250. At the beginning, the number of initial chromosomes for the superior clan and inferior clan are 83.3 and 166.7, respectively, on average. But within a short time, the number of chromosomes in the superior clan rapidly exceeds those in the inferior clan, thus constituting the majority of chromosomes in the population. Finally, at the 20th generation, the number of chromosomes in the superior clan and inferior clan have changed to 213.5 and 36.5, respectively, on average. In theory, such a situation is disadvantageous to our method, since the fitness of chromosomes in the superior clan needs to be calculated. But as a result of elitism, at each generation, half the chromosomes are coming directly from the former generation and almost all belong to the superior clan. For these chromosomes, their fitness is known and need not be re-calculated. Hence for our method, the number of MSE computations can be substantially decreased. Furthermore, for the remaining 125 chromosomes generated by way of the GA evolutionary process, about half belong to the inferior clan and their fitness need not be calculated for our method. But for the traditional GA method, the fitness of all 125 chromosomes need to be calculated. Overall, in comparison with the traditional GA method, about 45% of MSE computations can be saved. The result is almost the same as one that ignores elitism. Table 6 also lists the comparative results of the proposed GA method, traditional GA method and full search method for Lena, Pepper, F16, and Baboon. The size of the population is 250 and the pre-set number of iterations is 20. For proper comparison, the two GA methods are performed over 10 runs and the results are averaged. For the number of MSE computations, the proposed GA method saves about 46–54% from the traditional GA method. In terms of the speed up rate, the encoding time of the proposed GA method is some 24–30% less than that of the traditional GA method, while the penalty is only 0.2–0.43 dB decay on retrieved image quality. Remember too that the encoding time of our method already includes the time spent on the DCT calculation and classification of all the image blocks. Finally, in comparison with the full search method, the encoding speed of our proposed GA method is some 130 times faster than that of the full search method, while the retrieved image quality remains acceptable.
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
1159
Fig. 6. The number of chromosomes in the superior clan and inferior clan versus the number of generations for Lena – with elitism. Table 6 The comparison of proposed GA method and tradition GA method together with full search method – with elitism. Method
Time (sec)
MSE computations
PSNR (dB)
Lena
Full search Proposed GA Traditional GA
2589.53 19.47 25.50
475,799,552 1,463,472 2,688,000
28.91 27.03 27.46
Pepper
Full search Proposed GA Traditional GA
2587.91 18.05 25.70
475,799,552 1,240,140 2,688,000
29.84 27.90 28.29
F16
Full search Proposed GA Traditional GA
2582.92 18.59 25.64
475,799,552 1,329,282 2,688,000
25.24 23.86 24.11
Baboon
Full search Proposed GA Traditional GA
2591.36 18.55 25.83
475,799,552 1,319,639 2,688,000
20.15 19.52 19.72
Results for the retrieved images by our proposed GA, traditional GA, and full search methods are shown in Fig. 7. The retrieved images of the two GA methods are obtained by decoding two compression codes selected respectively from the 10 runs of the two GA methods in Table 6. Fig. 7(a) is the original Lena image and Fig. 7(b) is the initial image used for fractal decode. Figs. 7(c) and 7(d) are the retrieved images using the full search method and traditional GA method, respectively. Fig. 7(e) is the retrieved image using the proposed GA method. 5. Conclusion In this paper, a genetic algorithm with hybrid select mechanism is proposed to speed up the encoding speed of a fractal image. First, all the image blocks are classified into three classes: S type, H type, and D type according to their DCT coefficients. Then, during the GA evolution, the population at every generation is separated into two clans: a superior clan and an inferior clan according to whether the type of chromosome is the same as that of the range block to be encoded or not. The hybrid select mechanism proposed by us is used to select appropriate parents from the two clans in order to reduce the amount of MSE computations while maintaining the retrieved image quality. First, an experiment using GA without elitism is performed. For the number of MSE computations, our proposed GA method requires 46–55% less than that of the traditional GA method. For the speed up rate, the proposed GA method can save about 25–32% from the encoding time of the traditional GA method. Second, for the experiment using GA with elitism, the number of MSE computations of the proposed GA method is some 46–54% fewer than that of the traditional GA method. For the speed up rate, the proposed GA method can save about 24–30% from the encoding time of the traditional GA method. For retrieved image quality, the experiment also shows that, whether the GA utilizes elitism or not, the results of our proposed GA method is almost identical to the traditional GA method, or has only a little decay. Finally, a comparison between our proposed GA method and a full search method is carried out. For GA without elitism, the encoding speed of our method is about 75 times faster than that of the full search method. Even more significantly, for GA with elitism, the encoding speed is some 130 times that of the full search method, while the retrieved image quality remains acceptable. There are, in addition, other soft computing methods such as evolutionary strategy, genetic programming, and particle swarm optimization etc. can be use to speedup the fractal encoder. Our next effort would be to find some more efficient algorithms for fractal image compression by using the soft computing methods.
1160
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
Fig. 7. (a) Original image, Lena of size 256 × 256. (b) Initial image for the decoder of the fractal compression. (c) Full searching method, MSE computations = 475,799,552, PSNR = 28.91, time = 2589.53 sec. (d) Traditional GA method, MSE computations = 2,688,000, PSNR = 27.44, time = 25.52 sec. (e) Proposed GA method, MSE computations = 1,467,052, PSNR = 27.17, time = 19.48 sec.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]
M.F. Barnsley, S. Demko, Iterated function systems and the global construction of fractals, Proc. Roy. Soc. A 399 (1985) 243–275. M.F. Barnsley, Fractal Everywhere, Academic, New York, 1988. M.F. Barnsley, J.H. Elton, D.P. Hardin, Recurrent iterated function system, in: Constructive Approximation, 1989, pp. 3–31. A.E. Jacquin, Image coding based on a fractal theory of iterated contractive image transformations, IEEE Trans. Image Process. 1 (1) (1992) 18–30. C.T. Wang, et al., Detecting and restoring the tampered images based on iteration-free fractal compression, J. Syst. Software 67 (2) (2003) 131–140. D.R. Chen, et al., Classification of breast ultrasound images using fractal feature, Clin. Imaging 29 (4) (2005) 235–245. T.A. Linnell, F. Derari, Mapping vector accumulator: fractal domain feature for character recognition, Electron. Lett. 40 (22) (2004) 1406–1407. C.H. Li, S.S. Wang, Digital watermarking using fractal image coding, IEICE Trans. Fund. E83-A (6) (2000) 1286–1288. Z. Yao, Fixed point in fractal image compression as watermarking, Int. Conf. Image Process. 2 (2003) 475–478. Y. Fisher, Fractal Image Compression, Theory and Application, Springer-Verlag, New York, 1994. F. Davoine, M. Antonini, J.M. Chassery, M. Barlaud, Fractal image compression based on Delaunay triangulation and vector quantization, IEEE Trans. Image Process. 5 (2) (1996) 338–346. D. Saupe, S. Jacob, Variance-based quadtrees in fractal image compression, Electron. Lett. 33 (1) (1997) 46–48. T.K. Truong, C.M. Kung, J.H. Jeng, M.L. Hsieh, Fast fractal image compression using spatial correlation, Chaos Solitons Fractals 22 (5) (2004) 1071–1076. S. Furao, O. Hasegawa, A fast no search fractal image coding method, Signal Process. Image Commun. 19 (5) (2004) 393–404. J.H. Holland, Adaptation in Natural and Artificial Systems, The University of Michigan Press, Ann Arbor, 1975. D.E. Goldberg, Genetic Algorithms in Search, Optimization, and Machine Learning, Addison–Wesley, New York, 1989.
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
1161
[17] Y. Du, K.T. Chan, Feasibility of applying genetic algorithms in space-time block coded multiuser detection systems, Digital Signal Process. 15 (2005) 161–170. [18] S.H. Doong, C.C. Lai, C.H. Wu, Genetic subgradient method for solving location allocation problems, Appl. Soft Comput. 7 (2007) 373–386. [19] M. Wang, S. Yang, S. Wu, A GA-based UWB pulse waveform design method, Digital Signal Process. 18 (1) (2008) 65–74. [20] S.T. Pan, C.C. Lai, Identification of chaotic systems by neural network with hybrid learning algorithm, Chaos Solitons Fractals 37 (1) (2008) 233–244. [21] S.T. Pan, A canonic-signed-digit coded genetic algorithm for designing finite impulse response digital filter, Digital Signal Process. 20 (2) (2010) 314– 327. [22] S.K. Mitra, C.A. Murthy, M.K. Kundu, Technique for fractal image compression using genetic algorithm, IEEE Trans. Image Process. 7 (4) (1998) 586–593. [23] M.S. Wu, J.H. Jeng, J.G. Hsieh, Schema genetic algorithm for fractal image compression, Eng. Appl. Artif. Intell. 20 (4) (2007) 531–538.
Ming-Sheng Wu received the B.S. degree from the Department of Electronic Engineering of Chung Yuan Christian University, Chung Li, Taiwan, in 1988, and the M.S. and Ph.D. degrees from National Sun Yat-Sen University, Kaohsiung, Taiwan, in 1990 and 2006, respectively, all in Electrical Engineering. He is currently an Associate Professor of the Department of Electrical Engineering, Cheng-Shiu University, Kaohsiung, Taiwan. His research includes image processing and image compression. Yih-Lon Lin received the B.S. degree from the Department of Electronic Engineering, I-Shou University, Kaohsiung, Taiwan, in 1997 and the M.S. and Ph.D. degrees from the Department of Electrical Engineering, National Sun Yat-Sen University, Kaohsiung, Taiwan, in 1999 and 2006, respectively, all in electrical engineering. Currently, he is an Assistant Professor at the Department of Information Engineering, I-Shou University. His research interests include machine learning, multimedia applications, and soft computing.
Contents lists available at ScienceDirect
Digital Signal Processing www.elsevier.com/locate/dsp
Genetic algorithm with a hybrid select mechanism for fractal image compression Ming-Sheng Wu a , Yih-Lon Lin b,∗ a b
Department of Electrical Engineering, Cheng Shiu University, Kaohsiung, Taiwan Department of Information Engineering, I-Shou University, Kaohsiung, Taiwan
a r t i c l e
i n f o
a b s t r a c t
Article history: Available online 4 January 2010
In this paper, a genetic algorithm with a hybrid select mechanism is proposed to speed up the fractal encoder. First, all of the image blocks including domain blocks and range blocks are classified into three classes: smooth; horizontal/vertical edge; and diagonal/subdiagonal edge, according to their discrete cosine transformation (DCT) coefficients. Then, during the GA evolution, the population of every generation is separated into two clans: a superior clan and an inferior clan, according to whether the chromosome type is the same as that of the range block to be encoded or not. The hybrid select mechanism proposed by us is used to select appropriate parents from the two clans in order to reduce the number of MSE computations and maintain the retrieved image quality. Experimental results show that, since the number of MSE computations in the proposed GA method is about half of the traditional GA method, the encoding time for the proposed GA method is less than that of the traditional GA method. For retrieved image quality, the proposed GA method is almost the same as the traditional GA method or only has a little decay. Moreover, in comparison with the full search method, the encoding speed of the proposed GA method is some 130 times faster than that of the full search method, whereas the retrieved image quality is still relatively acceptable. © 2009 Elsevier Inc. All rights reserved.
Keywords: Fractal image compression Discrete cosine transformation Classification method Genetic algorithm Hybrid select mechanism
1. Introduction Fractal image compression is potentially a great coding scheme since it features a high compression ratio and good retrieved image quality. The original idea was proposed by Barnsley [1–3] according to the iterated function system (IFS). In 1992, the first fractal image coding algorithm was realized by Jacquin [4]. The underlying premise of the coding scheme is based on the partitioned iteration function system (PIFS) which utilizes the self-similarity characteristic in a natural image to achieve compression. In addition to image compression, the application of fractals is extremely widespread, in such things as image restoration [5], medical image classification [6], character recognition [7], and watermarking [8,9]. It is therefore little wonder this technique attracts the attention of many researchers. The drawback to fractal image compression is that the encoding process is very time-consuming since, for each range block, a large number of computations of similarity measures must be done to find the best matched domain block. Hence the main research topic of the field is focused on how to reduce the encoding time. In the past, the techniques of partition and classification were used to reduce the number of computations of similarity measures to speedup the encoder [10– 12]. Then the idea of spatial correlation was also added to the fractal encoding algorithm to reduce the search space of the optimal solution. Truong et al. [13] limit the search space for the current range block on the neighborhoods of the
*
Corresponding author. Address for correspondence: No. 1, Sec. 1, Syuecheng Rd., Dashu Township, Kaohsiung County 840, Taiwan. Fax: +886 7 6578944. E-mail addresses: [email protected], [email protected] (Y.-L. Lin).
1051-2004/$ – see front matter doi:10.1016/j.dsp.2009.12.009
©
2009 Elsevier Inc. All rights reserved.
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
1151
matched domain blocks of the neighboring range blocks by utilizing the spatial correlations between neighboring blocks in both the domain pool and the range pool. In comparison with the full search algorithm, this algorithm achieves a 2.6 times speedup ratio yet preserves almost the same retrieved image quality. Furao and Hasegawa [14] present a no search fractal encoding algorithm. In this algorithm, they select a domain block which has the same center with the range block to be encoded as the matched one. This is possible since the two blocks usually have the same texture in a natural image. Utilizing such a characteristic and merging quad-tree technique, the authors have achieved the fastest encoding velocity to date, but unfortunately both the retrieved image quality and the compression ratio suffer as a result. Recent efforts to push the genetic algorithm (GA) to a numerical optimum solution have become highly popular among researchers. GA was developed by Holland (1975) [15] over the course of the 1960s and 1970s and popularized by Goldberg [16]. The genetic algorithm is a global search technique that mimics natural selection and natural genetics. GA is capable of solving many large, complex problems. Experience shows such problems are hard to solve, since their search spaces generally have very rough landscapes riddled with many local optima. Many researchers have used the advantages of GA to achieve some significant results [17–21]. Since the landscape of natural images is also very rough, GA is well suited to the search for the best match in fractal image compression. For example, Mitra [22] proposed a search strategy using GA in which the chromosome consists of the PIFS of all range blocks. In this method, the fitness of a chromosome is measured by the distance between the attractor and the given image. Although a faster speed can be achieved, the retrieved image quality is not good enough because the fitness is defined in terms of the whole image. Wu et al. [23] used a schema genetic algorithm (SGA) to speed up the fractal encoder. In SGA, the genetic operators are adapted according to the schema theorem in the evolutionary process performed on the range blocks. The encoding speed of such a method is about 100 times faster than that of the full search method and can preserve the retrieved image quality. In this paper, a genetic algorithm with a hybrid select mechanism is proposed to speed up fractal encoding. All of the image blocks are classified into three classes: smooth type; diagonal/sub-diagonal edge type; and horizontal/vertical edge type. Two DCT coefficients – the lowest vertical coefficient and the lowest horizontal coefficient – are used to determine which type the block belongs to. Then, for each range block to be encoded, GA is used to find the best matched solution. During the GA evolutionary process, all the chromosomes are classified as one of two clans: a superior clan or an inferior clan. In the superior clan, the type of chromosome is the same as that of the range block. In the inferior clan, it is not. If a chromosome belongs to the superior clan, it is a good candidate for the best match. Hence its fitness value will be calculated. The chromosomes belong to the superior clan are ranked according to their fitness values and placed at the front of the chromosomes belong to the inferior clan. In contrast, if a chromosome belongs to the inferior clan, it is not a good candidate for the best match. Hence its fitness value need not be calculated. The chromosomes belong to inferior clan will be assigned the same rank. Finally, two selection mechanisms – uniform random selection and ranking selection – are merged in the GA mating process to select appropriate parents in order to find a suitable domain block and speedup the encoder. By using the evolutionary technique, we find that the encoding speed can be improved, since the number of fitness computations is reduced. Moreover, since good chromosomes are ranked at the front, good genes can be carried into the next generation during the GA evolutionary process. In this way, the quality of the retrieved image can also be maintained. The rest of the paper is organized as follows. We introduce the conventional fractal image coding scheme in Section 2. Section 3 describes the proposed genetic algorithm in this paper. Some experimental simulations are performed in Section 4 to verify the improvement of our proposed algorithm. Finally, a conclusion is made in Section 5. 2. Fractal image encoding Fractal image compression is based on the local self-similarity property and PIFS. The related implementation process is stated as follows. Suppose the original gray level image f is of size m × m. To simplify the discussion, let m be 256. The domain pool is defined as the set of all possible blocks of size 16 × 16 of the image f , which makes up (256 − 16 + 1) × (256 − 16 + 1) = 58,081 blocks. The range pool is defined to be the set of all non-overlapping blocks of size 8 × 8 of the image f , which makes up (256/8) × (256/8) = 1024 blocks. For each block v from the range pool, the fractal transformation is constructed by searching all of the elements in the domain pool to find the block with the greatest similarity. Let u denote a subsampled domain block which is of the same size as v. The similarity of u and v is measured using mean square error (MSE), which is defined by
MSE =
7 7 1
64
2
u (i , j ) − v (i , j )
(1)
j =0 i =0
The fractal transformation allows Dihedral transformations of the domain blocks, i.e., the 8 orientations of the blocks generated by rotating the blocks counterclockwise at angles 0◦ , 90◦ , 180◦ and 270◦ and flipping with respect to the line y = x, respectively, as shown in Table 1. Thus for a given block from the range pool, there are 58,081 × 8 = 464,648 MSE computations to obtain the block with the greatest similarity from the domain pool. Thus, in total, one needs 1024 × 464,648 = 475,799,552 MSE computations to encode the whole image using this full search compression method.
1152
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
Table 1 The 8 transformations in the Dihedral group. Rotate 0◦
1 0
0 1
0 1
1 0
Rotate 90◦
0 1 −1 0
Rotate 180◦
−1 0
0 −1
Flip with the line x = y from above 1 0 0 −1
0
−1
−1 0
Rotate 270◦
0 −1 1 0
−1 0 0
1
The fractal transformation also allows the adjustment of the contrast p and the brightness q on the sub-sample block u. Thus the similarity is to minimize the quantity d = p · uk + q − v , where uk , 0 k 7 are the 8 orientations of u. By calculus, p and q can be computed directly by
p=
[ N u , v − u , 1 v , 1 ] , [ N u , u − u , 1 2 ]
q=
1 N
v , 1 − p u , 1
= [1 1 . . . 1] T . where N is the number of pixels of the range block and 1 Finally, the position x and y of the domain block, the contrast p, the brightness q, and the orientation k constitute the fractal code of the given range block v. For a 256 × 256 image, both x and y require 8 and 8 bits, respectively, to represent the position of the domain block. For contrast p, brightness q, and the orientation k, 5, 7 and 3 bits are required, respectively. Hence one needs 8 + 8 + 5 + 7 + 3 = 31 bits in total to encode a range block. Finally, as v runs over all 1024 range blocks in the range pool, the encoding process is completed. To decode, one first makes up the 1024 affine transformations from the compression codes. Next, one chooses any initial image and performs the 1024 affine transformations on the image to obtain a new image. The transformation proceeds recursively according to the Contractive Mapping Fixed-Point Theorem and Collage Theorem until the sequence of images converges to the encoded image. The stopping criterion of the recursion is designed in accordance with the user’s application and the final image is the retrieved image of fractal coding. 3. Genetic algorithm for fractal image compression In this section, a genetic algorithm with hybrid select mechanism proposed by us will be introduced. First, by using two DCT coefficients: the lowest vertical coefficient F (1, 0) and the lowest horizontal coefficient F (0, 1), all of the image blocks including domain blocks and range blocks are classified into three classes: smooth type, diagonal/sub-diagonal edge type, and horizontal/vertical edge type according to their texture. Then when the GA is used to encode each range block, a hybrid select mechanism composed of two different kinds of selection mechanisms are used to select the appropriate parents to speed up the encoder. Below, the classification method for the image blocks is introduced first. Second, the proposed GA algorithm designed to reduce encoding time will also be stated in detail. 3.1. The classification algorithm for image blocks The main goal of the classification algorithm is to classify all of the image blocks including domain blocks and range blocks into three classes. By using DCT transformation, an image block can be transferred from the spatial domain to the frequency domain, in which the DCT coefficients located in the upper-left represent the low frequency information of the image block and reflect the rough contour of the image block. In contrast, the DCT coefficients located in the lower-right represent the high frequency information of the image block and reflect the fine texture of the image block. Hence an image block can be determined as belonging to which type, only by using its lower frequency DCT coefficients. Based on the viewpoint stated above, we design a classifier which uses the lowest vertical coefficient F (1, 0) and the lowest horizontal coefficient F (0, 1) to execute the action of classification. Assume f is a given image block of size 8 × 8, F (1, 0) and F (0, 1) can be computed from
√ F (1, 0) =
2
8
√ F (0, 1) =
7
7
2
8
f (i , j ) cos θi
(2)
f (i , j ) cos θ j
(3)
j =0 i =0 7
7
j =0 i =0
where θi = (2i + 1)π /16, i = 0, 1, . . . , 7. The magnitude of F (1, 0) reflects the intensity variation between the left half and right half of image block f and the magnitude of F (0, 1) reflects the intensity variation between the upper half and lower half of image block f . For an image block f , if both | F (1, 0)| and | F (0, 1)| are small, it means that the block f has no clear texture. Hence the landscape on
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
1153
Fig. 1. The classification of H type blocks and non-H type blocks.
block f tends to smoothness. Block f will be classified as smooth type and denoted as S type. Second, if | F (1, 0)| is large and | F (0, 1)| is small, f will have a vertical edge. But if | F (1, 0)| is small and | F (0, 1)| is large, f will have a horizontal edge. According to Dihedral characteristics, for the two cases, f is classified as horizontal/vertical edge type and denoted as H type. Finally, if both | F (1, 0)| and | F (0, 1)| are large, it means that f will have a diagonal or sub-diagonal edge. Such f will be classified as diagonal/sub-diagonal type and denoted as D type. To determine the different types of image blocks, these blocks first execute the DCT transformation to find their F (1, 0) and F (0, 1) values. Then the classification algorithm is used to determine which type these image blocks belong to, according to their respective F (1, 0) and F (0, 1) values. First, H type blocks are separated out. For each image block, we take the absolute values of its F (1, 0) and F (0, 1) to constitute a DCT absolute coefficient pair (| F (1, 0)|, | F (0, 1)|). Then the DCT absolute coefficient pair is undergoes the normalization process. After normalization, the normalized DCT absolute coefficient pair (| F (1, 0)| N , | F (0, 1)| N ) is obtained and is labeled on a Cartesian coordinate system in which | F (1, 0)| N is the x-coordinate and | F (0, 1)| N is the y-coordinate. The result is shown in Fig. 1. Here | F (1, 0)| N and | F (0, 1)| N denote normalized | F (1, 0)| and| F (0, 1)|, respectively. We can find from Fig. 1 that the positions of all the normalized image blocks are located on the unit circle of the first quadrant. According to the texture characteristics of the image block, if the normalized DCT absolute coefficient pair of an image block f is near the point (1, 0), it will have a clear vertical edge since | F (1, 0)| N is large and | F (0, 1)| N is small. Hence the image block f will be classified as H type. For the same reason, if the normalized DCT absolute coefficient pair of an image block f is near the point (0, 1), it will have a clear horizontal edge since | F (1, 0)| N is small and | F (0, 1)| N is large. Such an image block f will also be classified as H type. Two threshold ◦ and T ◦ are used to separate H type blocks from non-H type blocks. Therefore, an image block is classified angles: T H1 H2 ◦ or between T ◦ and 90◦ . Once the H as H type if its position in Fig. 1 is located on the unit circle between 0◦ and T H1 H2 type blocks are determined, the S type blocks are also isolated from the remaining non-H type blocks. For the remaining image blocks, according to their DCT absolute coefficient pairs (| F (1, 0)|, | F (0, 1)|) which have not been normalized, we find an image block whose distance from the origin of the coordinate is greater than the others and label its coefficient pair as (| F (1, 0)|max , | F (0, 1)|max ). Then the DCT absolute coefficient pair (| F (1, 0)|, | F (0, 1)|) of the remaining image blocks is divided by ((| F (1, 0)|max )2 + (| F (0, 1)|max )2 )1/2 . Through the scaling process, the scaled DCT absolute coefficient pair (| F (1, 0)| R , | F (0, 1)| R ) is obtained and is re-labeled on a Cartesian coordinate system in which | F (1, 0)| R is the x-coordinate and | F (0, 1)| R is the y-coordinate. The result is shown in Fig. 2. In Fig. 2, (| F (1, 0)|max, R , | F (0, 1)|max, R ) is the image block whose distance from the origin of the coordinates is greatest. We find that, through the scaling process, the positions of the ◦ and T ◦ . If the position of an image block f is near rest of the image blocks are located in the unit circle between T H1 H2 the origin of the coordinate, it will have a smooth landscape since both | F (1, 0)| R and | F (0, 1)| R are small. Such an image block f will be classified as S type. A threshold circle whose radius is T S is used to separate the S type and non-S type blocks. Hence an image block f is classified as S type if its position is located in the threshold circle. Finally, the non-S type blocks will be classified as D type. For these non-S type blocks, they have a clear diagonal or sub-diagonal edge since their | F (1, 0)| R and | F (0, 1)| R values are large. The classification steps for the image blocks are discussed in detail below: 1. For all the image blocks, calculate the F (1, 0) and F (0, 1) from (2) and (3), respectively, and constitute the DCT absolute coefficient pair (| F (1, 0)|, | F (0, 1)|). 2. Normalize the DCT absolute coefficient pair (| F (1, 0)|, | F (0, 1)|) for all the image blocks to obtain the normalized DCT absolute coefficient pair(| F (1, 0)| N , | F (0, 1)| N ). Label the position of each image block according to its (| F (1, 0)| N , | F (0, 1)| N ) on the Cartesian coordinate system. The result is shown in Fig. 1.
1154
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
Fig. 2. The classification of S type blocks and D type blocks.
Fig. 3. Chromosome formation.
3. Determine the H type blocks from Fig. 1. An image block is classified as H type if its position is located on the unit ◦ or between T ◦ and 90◦ . circle between 0◦ and T H1 H2 4. For the remainder in Step 3, find the pair (| F (1, 0)|max , | F (0, 1)|max ) according to their DCT absolute coefficient pair (| F (1, 0)|, | F (0, 1)|) obtained in Step 1. Divide the (| F (1, 0)|, | F (0, 1)|) of the remaining blocks by ((| F (1, 0)|max )2 + (| F (0, 1)|max )2 )1/2 to obtain the scaled DCT absolute coefficient pair (| F (1, 0)| R , | F (0, 1)| R ). Label their position on the Cartesian coordinate system. The result is shown in Fig. 2. 5. Determine the S type blocks from Fig. 2. An image block is classified as S type if its position is located in a circle whose radius is T S . 6. Classify the rest of the blocks in Step 5 as D type. Once the classification process is complete, for each range block, a genetic algorithm with a hybrid select mechanism proposed by us will be used on the classified image blocks in order to search for the best matched domain block. 3.2. A genetic algorithm with a hybrid select mechanism for FIC Here, a genetic algorithm with hybrid select mechanism is introduced in detail. For each range block to be encoded, the GA algorithm executes GA evolution on the classified domain blocks obtained in Section 3.1 to find the best matched solution. Its main goal is to reduce the amount of computation of similarity measure between range block and sub-sample domain block in order to speed up the encoding speed while maintaining the retrieved image quality. The related evolutionary strategies with respect to the GA are detailed as follows. (1) Chromosome: Since the fractal encoding scheme utilizes the PIFS to encode every range block, one takes the absolute position (x, y ) of a domain block and the Dihedral transformation k to constitute a chromosome. A chromosome is 19 bits in length as shown in Fig. 3, in which x, y and k are the genes of 8, 8 and 3 bits, respectively. (2) Fitness function: The distance between range block and sub-sampled domain block is measured by MSE. The fitness value is defined as the reciprocal of MSE. (3) Initial population: Chromosomes are initialized randomly. (4) Selection: Here, a hybrid selection mechanism which is composed of uniform random selection and ranking selection is adopted to find appropriate parents and reduce the amount of MSE computation. All of the chromosomes are divided into two clans:
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
1155
Fig. 4. Uniform crossover.
a superior clan and an inferior clan. If the type of chromosome is the same as that of the range block to be encoded, it is a good candidate for the best match. It is proof that the gene in the chromosome is good. Hence it should have a higher probability to be selected as a parent in order to generate better offspring. Such a chromosome will be classified as belonging to the superior clan. In contrast, if the type of chromosome is not the same as that of the range block to be encoded, it is not a good candidate for the best match. Its gene is not good and hence the probability it is selected as a parent should be lower. Such a chromosome will be classified as belonging to the inferior clan. For those chromosomes belonging to the superior clan, their genes are good. We hope that these genes can be preserved to next generation to generate better offspring. Hence we let their fitness be calculated and be ranked from high to low according to their fitness values. But for those chromosomes belonging to the inferior clan, their genes are not good. These genes will be eliminated during the GA evolution. Therefore, for reducing the number of MSE computations, their fitness will not be calculated. We assign them the same rank and rank them directly behind the superior clan. For these chromosomes, their rank is equal to the number of chromosomes in the superior clan plus 1. Once all the chromosomes have been ranked, the hybrid select mechanism will be executed in order to select appropriate parents from the ranked chromosomes. First, the ranking selection is carried out to select parents from all the ranked chromosomes. If a chromosome belonging to the superior clan is selected, it will be immediately taken as the parent. But if the select mechanism determines the parent as coming from the inferior clan, the uniform random selection will be invoked in order to select a parent from the inferior clan. For example, suppose the size of the population is 7 and the number of chromosomes belonging to the superior clan is 3. First, the ranking select mechanism is used to find a parent from the population. For the 3 chromosomes in the superior clan, the selected probabilities according to their rank are 7/28, 6/28, and 5/28, in that order. The probability that the parent comes from the inferior clan is 4/28 + 3/28 + 2/28 + 1/28 = 10/28. By use of the random generator, if the random number lies between 0 and 7/28, the chromosome whose rank is 1 will be selected. For the same reason, if the random number lies between 7/28 and 13/28 or between 13/28 and 18/28, the chromosome whose rank is 2 or 3 will be selected, respectively. But if the random number lies between 18/28 and 1, the ranking select mechanism determines that the parent comes from the inferior clan. Hence the uniformly random selection will be invoked to select a parent from the 4 chromosomes in the inferior clan. On average, the selected probability for each chromosome in the inferior clan is (10/28)/4 = 5/56. The fact the parent comes from the inferior clan maintains the diversity of the population. (5) Crossover: The uniform crossover operation, as shown in Fig. 4, is adopted with probability P c to determine whether the operation will be performed or not. The temporary offspring are generated by executing the crossover operation on the two parents. (6) Mutation: The mutation operation with probability P m is applied to the temporary offspring in order to generate the offspring for the next generation. Here the goal of the mutation is to maintain the diversity of the population and to avoid pre-maturity. (7) Stopping criterion: The goal of stopping criterion is to stop the GA evolutionary process. Here, a pre-set number of iterations is adopted to stop the evolution. Connecting the above operational strategies, a genetic algorithm with a hybrid select mechanism is proposed. The detailed evolutionary steps regarding the GA are given as follows: Initially, set the population size, the crossover mask, crossover rate P c , mutation mask, and mutation rate P m . 1. Generate the initial population of chromosomes randomly. 2. Divide the chromosomes into two clans: superior clan and inferior clan according to their types. Calculate the fitness values of the chromosomes belonging to the superior clan. 3. Rank all the chromosomes. First, for the chromosomes belonging to the superior clan, rank them according to their fitness values. Second, let the chromosomes belonging to the inferior clan have the same rank, and their rank is the number of chromosomes in the superior clan plus 1.
1156
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
Table 2 ◦ and T ◦ . The rank condition in part for the determination of T H1 H2 Rank
min( D 1,0 , D 0,1 )
Phase angle
337 338 339 340 341 342 343 344 345 346
0.1921 0.1922 0.1925 0.1934 0.1968 0.1970 0.1980 0.1995 0.2000 0.2008
78.98◦ 11.03◦ 11.05◦ 78.90◦ 11.29◦ 11.30◦ 11.36◦ 11.45◦ 11.48◦ 11.53◦
4. If a pre-set number of iterations is reached, then stop and record the fractal code. Otherwise, go to the next step. 5. Execute the hybrid select mechanism to select parents. First execute the ranking selection. If the selected chromosome belongs to the superior clan, take it immediately as the parent. In contrast, if the selected chromosome belongs to the inferior clan, invoke the uniformly random selection to select a parent from the inferior clan. 6. Perform the uniform crossover to generate a temporary offspring. 7. Perform the mutation operation on the temporary offspring to generate the offspring for the next generation. 8. Go to Step 2. 4. Experimental results The images Lena, Baboon, F16, and Pepper are tested to demonstrate the speed-up rate and retrieved quality of the proposed GA method. For a given image of size 256 × 256, let the size of range block be 8 × 8. The software simulations are performed using BCB v. 6.0 on a Celeron Core2 Duo 2.0 GHz, 2 G RAM, windows XP PC. The quality of the retrieved image g in comparison to the original image f is measured by peak signal to noise ratio (PSNR) given by
PSNR = 10 × log
2552
1 m×n
n−1 m−1 j =0
i =0
( f (i , j ) − g (i , j ))2
where m × n is the image size. 4.1. The determination of threshold values for the classification method All the image blocks will first be classified before the GA evolutionary process is executed. Hence three threshold values: ◦ , T ◦ and T must be obtained in order to execute the classification algorithm. It is obvious that the number of domain T H1 S H2 blocks is much larger than that of range blocks. If we use the domain pool to determine the threshold values, it will produce a lot of computation overhead. Thus, to reduce the computation load, we use the range pool to determine the threshold values instead of the domain pool and partition the range pool into three classes of almost equal numbers of members. First, ◦ and T ◦ in Fig. 1 are obtained to separate H type and non-H type blocks from all the image two threshold angles: T H1 H2 blocks. For each normalized range block, since the more that are near point (1, 0) or point (0, 1), the clearer the vertical or horizontal texture is, we calculate the value min( D 1,0 , D 0,1 ) for all the range blocks and rank them according to value from small to large. Here D 1,0 = ((| F (1, 0)| N − 1)2 + (| F (0, 1)| N − 0)2 )1/2 represents the distance between the normalized range block and point (1, 0) and D 0,1 = ((| F (1, 0)| N − 0)2 + (| F (0, 1)| N − 1)2 )1/2 represents the distance between the normalized range block and point (0, 1). The min( D 1,0 , D 0,1 ) represents the minimum value of both D 1,0 and D 0,1 . Table 2 shows the rank condition in part and the test image is Lena. In the last row of Table 2, the rank of the range block is 346. Its min( D 1,0 , D 0,1 ) is 0.2008 and the phase angle of the line connecting the range block and origin is 11.53◦ . For an image of size 256 × 256, with the size of range block 8 × 8, the number of range blocks is (256/8) × (256/8) = 1024. We let the first ◦ and T ◦ in Fig. 1, T ◦ is assigned to the third of range blocks, about 342, be classified as H type. For the assignment of T H1 H2 H1 middle point between the last rank of the H type blocks and the 1st rank of non-H type blocks. Furthermore, we let both ◦ ◦ . Hence from ◦ the L 1 to x-coordinate and L 2 to y-coordinate have the same span. In other words, we keep T H2 = 90 − T H1 ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ Table 2, we obtain the result T H1 is (11.30 + 11.36 )/2 = 11.33 and T H2 is 90 − 11.33 = 78.67 . If we let the range ◦ is ((90◦ − 78.90◦ ) + 11.29◦ )/2 = 11.20◦ and T ◦ is 90◦ − 11.20◦ = 78.80◦ , blocks of the first 340 be classified as H type, T H1 H2 since the 340th range block is near the y-coordinate. Second, the threshold radius T S is determined from the remaining 1024 − 342 = 682 range blocks. For the same reason, for the remaining range blocks, since the closer to the origin (0, 0), the smoother the characteristics, we calculate the value (| F (1, 0)|2R + | F (0, 1)|2R )1/2 for all the remaining range blocks and rank them according to their value from small to large. Table 3 shows the rank condition in part. In the last row of Table 3, the rank of the range block is 346, from the remaining 682, and the distance from the origin is 0.0822. We take the first half of the range blocks as S type and take the second half of the range blocks as D type. Hence both the number of S type blocks and the number of D type blocks are 682/2 = 341. For the same reason, T S is assigned to the middle point between the last
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
1157
Table 3 The rank condition in part for the determination of T S . Rank
Distance with origin
337 338 339 340 341 342 343 344 345 346
0.0781 0.0793 0.0795 0.0797 0.0803 0.0805 0.0808 0.0813 0.0821 0.0822
Table 4 Three threshold values and the classification result of the domain blocks by using the proposed classification method.
Lena Pepper F16 Baboon
◦ T H1
◦ T H2
TS
The number of H type blocks
The number of S type blocks
The number of D type blocks
11.33◦ 12.49◦ 12.23◦ 13.01◦
78.67◦ 77.51◦ 77.77◦ 76.99◦
0.0804 0.0986 0.0459 0.1504
19,223 19,276 18,782 20,350
11,852 12,818 14,924 14,759
27,006 25,987 24,375 22,972
Fig. 5. The number of chromosomes in the superior clan and inferior clan versus the number of generations for Lena – without elitism.
rank of the S type blocks and the 1st rank of the D type blocks. Therefore T S is (0.0803 + 0.0805)/2 = 0.0804. Table 4 lists the three threshold values and the classification result of the domain blocks by using the proposed classification method for ◦ , T ◦ , and T are 11.33◦ , 78.67◦ , and 0.0804, respectively, and the number of Lena, Pepper, F16, and Baboon. For Lena, T H1 S H2 domain blocks of H type, S type, and D type are 19,223, 11,852, and 27,006, respectively. 4.2. Comparison of proposed GA method and traditional GA method – without elitism Some simulations are done to demonstrate the speed-up ratio and retrieved quality of the proposed GA method in comparison to the traditional GA method. For GA operations, the proposed GA method uses a hybrid select mechanism proposed by us and the traditional GA method uses a ranking select mechanism. The uniform crossover is adapted for both GA methods. Elitism is ignored. At the same time, the GA parameters of the two methods are identical, in which P c is 0.6 and P m is 0.01. Fig. 5 shows the relationship between the number of chromosomes in the superior clan and the number of chromosomes in the inferior clan versus the number of generations for the proposed GA method. The test image is Lena and the size of the population is 250. The two curves in Fig. 5 represent the number of chromosomes in the superior clan and the inferior clan in every generation, respectively, in which the values taken are an average of 1024 range blocks in an image. At the beginning, the number of chromosomes in the superior clan is less than the number of chromosomes in the inferior clan. The number of initial chromosomes for the superior clan and inferior clan are 83.8 and 166.2, respectively, on average. During the GA evolution, the number of chromosomes from the superior clan increases gradually and the number of chromosomes from the inferior clan decreases gradually. At the 8th generation, the number of chromosomes in the superior clan exceeds
1158
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
Table 5 The comparison of proposed GA method and tradition GA method together with full search method – without elitism. Method
Time (sec)
MSE computations
PSNR (dB)
Lena
Full search Proposed GA Traditional GA
2589.53 35.48 47.64
475,799,552 2,831,391 5,120,000
28.91 26.25 26.23
Pepper
Full search Proposed GA Traditional GA
2587.91 32.68 47.93
475,799,552 2,361,127 5,120,000
29.84 26.65 27.10
F16
Full search Proposed GA Traditional GA
2582.92 34.25 47.70
475,799,552 2,653,875 5,120,000
25.24 23.01 23.32
Baboon
Full search Proposed GA Traditional GA
2591.36 33.45 48.32
475,799,552 2,501,905 5,120,000
20.15 19.34 19.47
the number of chromosomes in the inferior clan. Finally, at the 20th generation, the number of chromosomes from the superior clan and inferior clan have changed to 182.9 and 67.1, respectively, on average. We find from Fig. 5 that since the fitness for chromosomes in the inferior clan need not be calculated, we can save about two-thirds of the MSE computations at the GA evolutionary initial stage. Even at a later stage of the GA evolutionary process, we can still save about one-quarter of the MSE computations. On average, in comparison with the traditional GA method, about 45% of MSE computations can be saved. Table 5 lists some comparative results of our proposed GA method, the traditional GA method and the full search method for Lena, Pepper, F16, and Baboon. The size of the population is 250 and the pre-set number of iterations is 20. For a proper comparison, the two GA methods are performed over 10 runs and the results are averaged. For Lena, comparing our proposed GA method with the traditional GA method, the PSNRs of the two methods are almost the same, but for the proposed GA method, the number of MSE computations is only 55% that of the traditional GA method. In terms of the speed up rate, the encoding time of the proposed GA method can save about 25% off that of the traditional GA method, in which the encoding time of our method already includes the time spent on the DCT calculation and classification of all the image blocks. Moreover, for Pepper, F16, and Baboon, the encoding time of the proposed GA method saves about 32%, 28%, and 31%, respectively, from that of the traditional GA method, and the penalty is only 0.13–0.45 dB decay on retrieved image quality. Finally, in comparison with the full search method, the encoding speed of the proposed GA method is about 75 times faster than that of the full search method. 4.3. Comparison of proposed GA method and traditional GA method – with elitism In this section, elitism is considered. The relevant GA operators and parameters are the same as in Section 4.2. Fig. 6 shows both the number of chromosomes in the superior clan and the number of chromosomes in the inferior clan versus the number of generations for the proposed GA method. The test image is also Lena and the size of the population is 250. At the beginning, the number of initial chromosomes for the superior clan and inferior clan are 83.3 and 166.7, respectively, on average. But within a short time, the number of chromosomes in the superior clan rapidly exceeds those in the inferior clan, thus constituting the majority of chromosomes in the population. Finally, at the 20th generation, the number of chromosomes in the superior clan and inferior clan have changed to 213.5 and 36.5, respectively, on average. In theory, such a situation is disadvantageous to our method, since the fitness of chromosomes in the superior clan needs to be calculated. But as a result of elitism, at each generation, half the chromosomes are coming directly from the former generation and almost all belong to the superior clan. For these chromosomes, their fitness is known and need not be re-calculated. Hence for our method, the number of MSE computations can be substantially decreased. Furthermore, for the remaining 125 chromosomes generated by way of the GA evolutionary process, about half belong to the inferior clan and their fitness need not be calculated for our method. But for the traditional GA method, the fitness of all 125 chromosomes need to be calculated. Overall, in comparison with the traditional GA method, about 45% of MSE computations can be saved. The result is almost the same as one that ignores elitism. Table 6 also lists the comparative results of the proposed GA method, traditional GA method and full search method for Lena, Pepper, F16, and Baboon. The size of the population is 250 and the pre-set number of iterations is 20. For proper comparison, the two GA methods are performed over 10 runs and the results are averaged. For the number of MSE computations, the proposed GA method saves about 46–54% from the traditional GA method. In terms of the speed up rate, the encoding time of the proposed GA method is some 24–30% less than that of the traditional GA method, while the penalty is only 0.2–0.43 dB decay on retrieved image quality. Remember too that the encoding time of our method already includes the time spent on the DCT calculation and classification of all the image blocks. Finally, in comparison with the full search method, the encoding speed of our proposed GA method is some 130 times faster than that of the full search method, while the retrieved image quality remains acceptable.
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
1159
Fig. 6. The number of chromosomes in the superior clan and inferior clan versus the number of generations for Lena – with elitism. Table 6 The comparison of proposed GA method and tradition GA method together with full search method – with elitism. Method
Time (sec)
MSE computations
PSNR (dB)
Lena
Full search Proposed GA Traditional GA
2589.53 19.47 25.50
475,799,552 1,463,472 2,688,000
28.91 27.03 27.46
Pepper
Full search Proposed GA Traditional GA
2587.91 18.05 25.70
475,799,552 1,240,140 2,688,000
29.84 27.90 28.29
F16
Full search Proposed GA Traditional GA
2582.92 18.59 25.64
475,799,552 1,329,282 2,688,000
25.24 23.86 24.11
Baboon
Full search Proposed GA Traditional GA
2591.36 18.55 25.83
475,799,552 1,319,639 2,688,000
20.15 19.52 19.72
Results for the retrieved images by our proposed GA, traditional GA, and full search methods are shown in Fig. 7. The retrieved images of the two GA methods are obtained by decoding two compression codes selected respectively from the 10 runs of the two GA methods in Table 6. Fig. 7(a) is the original Lena image and Fig. 7(b) is the initial image used for fractal decode. Figs. 7(c) and 7(d) are the retrieved images using the full search method and traditional GA method, respectively. Fig. 7(e) is the retrieved image using the proposed GA method. 5. Conclusion In this paper, a genetic algorithm with hybrid select mechanism is proposed to speed up the encoding speed of a fractal image. First, all the image blocks are classified into three classes: S type, H type, and D type according to their DCT coefficients. Then, during the GA evolution, the population at every generation is separated into two clans: a superior clan and an inferior clan according to whether the type of chromosome is the same as that of the range block to be encoded or not. The hybrid select mechanism proposed by us is used to select appropriate parents from the two clans in order to reduce the amount of MSE computations while maintaining the retrieved image quality. First, an experiment using GA without elitism is performed. For the number of MSE computations, our proposed GA method requires 46–55% less than that of the traditional GA method. For the speed up rate, the proposed GA method can save about 25–32% from the encoding time of the traditional GA method. Second, for the experiment using GA with elitism, the number of MSE computations of the proposed GA method is some 46–54% fewer than that of the traditional GA method. For the speed up rate, the proposed GA method can save about 24–30% from the encoding time of the traditional GA method. For retrieved image quality, the experiment also shows that, whether the GA utilizes elitism or not, the results of our proposed GA method is almost identical to the traditional GA method, or has only a little decay. Finally, a comparison between our proposed GA method and a full search method is carried out. For GA without elitism, the encoding speed of our method is about 75 times faster than that of the full search method. Even more significantly, for GA with elitism, the encoding speed is some 130 times that of the full search method, while the retrieved image quality remains acceptable. There are, in addition, other soft computing methods such as evolutionary strategy, genetic programming, and particle swarm optimization etc. can be use to speedup the fractal encoder. Our next effort would be to find some more efficient algorithms for fractal image compression by using the soft computing methods.
1160
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
Fig. 7. (a) Original image, Lena of size 256 × 256. (b) Initial image for the decoder of the fractal compression. (c) Full searching method, MSE computations = 475,799,552, PSNR = 28.91, time = 2589.53 sec. (d) Traditional GA method, MSE computations = 2,688,000, PSNR = 27.44, time = 25.52 sec. (e) Proposed GA method, MSE computations = 1,467,052, PSNR = 27.17, time = 19.48 sec.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]
M.F. Barnsley, S. Demko, Iterated function systems and the global construction of fractals, Proc. Roy. Soc. A 399 (1985) 243–275. M.F. Barnsley, Fractal Everywhere, Academic, New York, 1988. M.F. Barnsley, J.H. Elton, D.P. Hardin, Recurrent iterated function system, in: Constructive Approximation, 1989, pp. 3–31. A.E. Jacquin, Image coding based on a fractal theory of iterated contractive image transformations, IEEE Trans. Image Process. 1 (1) (1992) 18–30. C.T. Wang, et al., Detecting and restoring the tampered images based on iteration-free fractal compression, J. Syst. Software 67 (2) (2003) 131–140. D.R. Chen, et al., Classification of breast ultrasound images using fractal feature, Clin. Imaging 29 (4) (2005) 235–245. T.A. Linnell, F. Derari, Mapping vector accumulator: fractal domain feature for character recognition, Electron. Lett. 40 (22) (2004) 1406–1407. C.H. Li, S.S. Wang, Digital watermarking using fractal image coding, IEICE Trans. Fund. E83-A (6) (2000) 1286–1288. Z. Yao, Fixed point in fractal image compression as watermarking, Int. Conf. Image Process. 2 (2003) 475–478. Y. Fisher, Fractal Image Compression, Theory and Application, Springer-Verlag, New York, 1994. F. Davoine, M. Antonini, J.M. Chassery, M. Barlaud, Fractal image compression based on Delaunay triangulation and vector quantization, IEEE Trans. Image Process. 5 (2) (1996) 338–346. D. Saupe, S. Jacob, Variance-based quadtrees in fractal image compression, Electron. Lett. 33 (1) (1997) 46–48. T.K. Truong, C.M. Kung, J.H. Jeng, M.L. Hsieh, Fast fractal image compression using spatial correlation, Chaos Solitons Fractals 22 (5) (2004) 1071–1076. S. Furao, O. Hasegawa, A fast no search fractal image coding method, Signal Process. Image Commun. 19 (5) (2004) 393–404. J.H. Holland, Adaptation in Natural and Artificial Systems, The University of Michigan Press, Ann Arbor, 1975. D.E. Goldberg, Genetic Algorithms in Search, Optimization, and Machine Learning, Addison–Wesley, New York, 1989.
M.-S. Wu, Y.-L. Lin / Digital Signal Processing 20 (2010) 1150–1161
1161
[17] Y. Du, K.T. Chan, Feasibility of applying genetic algorithms in space-time block coded multiuser detection systems, Digital Signal Process. 15 (2005) 161–170. [18] S.H. Doong, C.C. Lai, C.H. Wu, Genetic subgradient method for solving location allocation problems, Appl. Soft Comput. 7 (2007) 373–386. [19] M. Wang, S. Yang, S. Wu, A GA-based UWB pulse waveform design method, Digital Signal Process. 18 (1) (2008) 65–74. [20] S.T. Pan, C.C. Lai, Identification of chaotic systems by neural network with hybrid learning algorithm, Chaos Solitons Fractals 37 (1) (2008) 233–244. [21] S.T. Pan, A canonic-signed-digit coded genetic algorithm for designing finite impulse response digital filter, Digital Signal Process. 20 (2) (2010) 314– 327. [22] S.K. Mitra, C.A. Murthy, M.K. Kundu, Technique for fractal image compression using genetic algorithm, IEEE Trans. Image Process. 7 (4) (1998) 586–593. [23] M.S. Wu, J.H. Jeng, J.G. Hsieh, Schema genetic algorithm for fractal image compression, Eng. Appl. Artif. Intell. 20 (4) (2007) 531–538.
Ming-Sheng Wu received the B.S. degree from the Department of Electronic Engineering of Chung Yuan Christian University, Chung Li, Taiwan, in 1988, and the M.S. and Ph.D. degrees from National Sun Yat-Sen University, Kaohsiung, Taiwan, in 1990 and 2006, respectively, all in Electrical Engineering. He is currently an Associate Professor of the Department of Electrical Engineering, Cheng-Shiu University, Kaohsiung, Taiwan. His research includes image processing and image compression. Yih-Lon Lin received the B.S. degree from the Department of Electronic Engineering, I-Shou University, Kaohsiung, Taiwan, in 1997 and the M.S. and Ph.D. degrees from the Department of Electrical Engineering, National Sun Yat-Sen University, Kaohsiung, Taiwan, in 1999 and 2006, respectively, all in electrical engineering. Currently, he is an Assistant Professor at the Department of Information Engineering, I-Shou University. His research interests include machine learning, multimedia applications, and soft computing.