# A New Texture Synthesis Algorithm Based on Wavelet Packet Tree

Authors: Hsi Chin Hsin; Tze-Yun Sung; Yaw-Shih Shieh; Carlo Cattani
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

This paper presents an efficient texture synthesis based on wavelet packet tree (TSWPT). It has the advantage of using a multiresolution representation with a greater diversity of bases functions for the nonlinear time series applications such as fractal images. The input image is decomposed into wavelet packet coefficients, which are rearranged and organized to form hierarchical trees called wavelet packet trees. A 2-step matching, that is, coarse matching based on low-frequency wavelet packet coefficients followed by fine matching based on middle-high-frequency wavelet packet coefficients, is proposed for texture synthesis. Experimental results show that the TSWPT algorithm is preferable, especially in terms of computation time.


---

## Body

## 1. Introduction

Texture modeling can be effectively applied to a wide variety of natural surfaces such as plants, furs, skins, minerals, terrains, and fractal materials [1, 2] and is an important issue in cyber-physical systems [3–7]. Numerous techniques have been proposed for texture processing; one may refer to [8] for a complete survey.Given an example texture, the goal of texture synthesis is to produce a visually similar image of any size. One may easily tile small textures to synthesize a larger image; however, there are some blocking effects near the tile edges [9]. Although some smoothing methods were proposed to reduce the blocking effects at the cost of computation time, it seems to be limited in improvements, especially for structured textures [10]. Efros and Leung proposed a neighborhood matching method [11], in which each pixel of the synthesis image was obtained by searching the most similar one in the source image based on the user-defined neighboring pixels. Wei and Levoy took account of the order in which pixels were synthesized and proposed an order-independent search-based texture synthesis algorithm [12]. Instead of synthesizing one pixel at a time, a patch of pixels can be taken as a whole and synthesized by matching the overlap regions between neighboring patches [13–15]. To improve the synthesis speed, the multi-resolution approach has been widely used for texture analysis and synthesis [16, 17]; Fang and Lien developed a rapid image synthesis system by adopting a multi-resolution approach, which consisted of an analysis process and a synthesis process. It consumed most of the computation time in the analysis process; yet the speed of synthesis was very fast [16]. In [17], De Bonet proposed a scheme of generating the synthesis image by sampling the filtered outputs of a texture in the framework of Laplacian pyramid [18]. Burt adopted the Gaussian pyramid [19] to represent both the input texture and the synthesis image at multiple resolutions and synthesized the texture images from lower to higher resolutions [20].Wavelet transform provides an efficient multi-resolution representation [21], in which the higher frequency components of an image are projected onto the shorter basis functions with higher spatial resolutions and the lower frequency components are projected onto the larger basis functions with higher spectral resolutions. Such a compact representation matches the characteristics of human visual system [22]. Wavelet theory has been successfully applied to many applications such as parameter estimation in fractal signals and images [23–29]. Yu et al. randomly sampled blocks of wavelet coefficients from the input texture to substitute that of the synthesis image [30]. Cui et al. adopted a 2-level wavelet transform and generated the synthesis image by minimizing the sum of squared distances between neighboring blocks of wavelet coefficients [31].For images with textures, lots of wavelet coefficients are likely to be significant in the middle-high frequency subbands, which surely demand further decompositions for a more compact representation [32]. Note that wavelet transform only decomposes the low-frequency component of an image at each resolution. However, both the low-frequency and high-frequency components can be decomposed using wavelet packet transform, which provides more bases functions than wavelet transform [33]. In [34], we proposed an efficient scheme to organize the wavelet packet coefficients of an image into hierarchical trees called wavelet packet (WP) trees for image compression. In this paper, an efficient WP tree-based algorithm is proposed for the texture synthesis applications.The remainder of this paper proceeds as follows. In Section2, wavelet transform and wavelet packet transform are briefly reviewed. Section 3 describes the proposed scheme to synthesize a texture image based on the WP trees of an example texture. Experimental results are presented in Section 4. Conclusion can be found in Section 5.


## 2. Review of Wavelet Transform and Wavelet Packet Transform

Wavelet-transform- (WT-) based multiresolution analysis/synthesis has drawn a lot of attention to the signal/image/video applications. The extension of WT known as wavelet packet transform (WPT) provides a much larger family of bases functions with a more compact representation. In this section, a brief review of WT and WPT is given. For a more complete survey, we refer readers to [21].


### 2.1. Wavelet Transform

WT has a variety of desirable properties, for example, joint space-spatial frequency localization, self-similarity across subbands of the same orientation, and energy clustering within each subband. For a discrete signal at resolutionℓ, Sℓ(n), the wavelet transform is defined as(2.1)Sl+1(n)=∑kh(2n-k)⋅Sl(k),Dl+1(n)=∑kg(2n-k)⋅Sl(k),
where(2.2)h(n)=〈ϕ,ϕ-1,-n〉,g(n)=〈ψ,ϕ-1,-n〉,ϕ-1,-n(x)=2-1/2ϕ(2-1x-n).ψ is a (mother) wavelet, ϕ is the corresponding scaling function, Sℓ+1(n) is the approximation signal at the next coarser resolution ℓ+1, Dℓ+1(n) is the detail information between resolutions ℓ and ℓ+1, h(n) and g(n) are low-pass filter and high-pass filter, respectively, and 〈·,·〉 is an inner product operator. Sℓ(n) can be exactly reconstructed from Sℓ+1(n) and Dℓ+1(n) by using the following inverse wavelet transform (IWT):(2.3)Sl(n)=∑kh̃(n-2k)⋅Sl+1(k)+∑kg̃(n-2k)⋅Dl+1(k),
where h̃(n)=h(-n) and g̃(n)=g(-n).For the image applications, 2D WT can be obtained by using a tensor product of two 1D WTs, horizontally followed by vertically or vice versa. Specifically, letLL0(m,n) be an image at the finest resolution 0, where m and n are indices for the vertical and horizontal orientations, respectively. The 2D WT of LL0(m,n) is as follows:(2.4)LL1(m,n)=∑i∑jh(i)⋅h(j)⋅LL0(2m-i,2n-j),LH1(m,n)=∑i∑jh(i)⋅g(j)⋅LL0(2m-i,2n-j),HL1(m,n)=∑i∑jg(i)⋅h(j)⋅LL0(2m-i,2n-j),HH1(m,n)=∑i∑jg(i)⋅g(j)⋅LL0(2m-i,2n-j),
where LL1(m,n) is the approximation image at the next coarser resolution 1, LH1(m,n), HL1(m,n), and HH1(m,n) are the detail images in the vertical, horizontal, and diagonal orientations, respectively. Figure 1 shows a 2-level 2D WT, where subbands LHℓ, HLℓ, HHℓ; ℓ=1,2, and LL2 are delimited by solid lines. Similarly, 2D IWT can be obtained by using the tensor product of two 1D IWTs.Figure 1
2-level wavelet transform (solid boundaries), wavelet packet transform (dashed boundaries), and a wavelet packet tree with wavelet packet coefficients connected by arrows.


### 2.2. Wavelet Packet Transform

2D WT is only to decompose the lowest-frequency subband of an image in an iterative manner. More specifically, only the scaling coefficients are decomposed from higher to lower resolutions. However, for the texture applications, wavelet coefficients in the middle- and high-frequency subbands are likely to be significant, which needs to be taken into account to improve the multiresolution representation.As one can see, both the low-frequency scaling coefficients and high-frequency wavelet coefficients of a signal, at any resolution, can be decomposed, which leads to wavelet packet transform (WPT), and a much larger family of bases functions can be produced [33]. Moreover, 2D WPT can be obtained by using a tensor product of two 1D WPTs. Figure 1 shows a 2-level 2D WPT, where all the wavelet subbands LH1, HL1, and HH1, at resolution 1, are further decomposed into wavelet packet subbands delimited by dashed lines.


## 2.1. Wavelet Transform

WT has a variety of desirable properties, for example, joint space-spatial frequency localization, self-similarity across subbands of the same orientation, and energy clustering within each subband. For a discrete signal at resolutionℓ, Sℓ(n), the wavelet transform is defined as(2.1)Sl+1(n)=∑kh(2n-k)⋅Sl(k),Dl+1(n)=∑kg(2n-k)⋅Sl(k),
where(2.2)h(n)=〈ϕ,ϕ-1,-n〉,g(n)=〈ψ,ϕ-1,-n〉,ϕ-1,-n(x)=2-1/2ϕ(2-1x-n).ψ is a (mother) wavelet, ϕ is the corresponding scaling function, Sℓ+1(n) is the approximation signal at the next coarser resolution ℓ+1, Dℓ+1(n) is the detail information between resolutions ℓ and ℓ+1, h(n) and g(n) are low-pass filter and high-pass filter, respectively, and 〈·,·〉 is an inner product operator. Sℓ(n) can be exactly reconstructed from Sℓ+1(n) and Dℓ+1(n) by using the following inverse wavelet transform (IWT):(2.3)Sl(n)=∑kh̃(n-2k)⋅Sl+1(k)+∑kg̃(n-2k)⋅Dl+1(k),
where h̃(n)=h(-n) and g̃(n)=g(-n).For the image applications, 2D WT can be obtained by using a tensor product of two 1D WTs, horizontally followed by vertically or vice versa. Specifically, letLL0(m,n) be an image at the finest resolution 0, where m and n are indices for the vertical and horizontal orientations, respectively. The 2D WT of LL0(m,n) is as follows:(2.4)LL1(m,n)=∑i∑jh(i)⋅h(j)⋅LL0(2m-i,2n-j),LH1(m,n)=∑i∑jh(i)⋅g(j)⋅LL0(2m-i,2n-j),HL1(m,n)=∑i∑jg(i)⋅h(j)⋅LL0(2m-i,2n-j),HH1(m,n)=∑i∑jg(i)⋅g(j)⋅LL0(2m-i,2n-j),
where LL1(m,n) is the approximation image at the next coarser resolution 1, LH1(m,n), HL1(m,n), and HH1(m,n) are the detail images in the vertical, horizontal, and diagonal orientations, respectively. Figure 1 shows a 2-level 2D WT, where subbands LHℓ, HLℓ, HHℓ; ℓ=1,2, and LL2 are delimited by solid lines. Similarly, 2D IWT can be obtained by using the tensor product of two 1D IWTs.Figure 1
2-level wavelet transform (solid boundaries), wavelet packet transform (dashed boundaries), and a wavelet packet tree with wavelet packet coefficients connected by arrows.


## 2.2. Wavelet Packet Transform

2D WT is only to decompose the lowest-frequency subband of an image in an iterative manner. More specifically, only the scaling coefficients are decomposed from higher to lower resolutions. However, for the texture applications, wavelet coefficients in the middle- and high-frequency subbands are likely to be significant, which needs to be taken into account to improve the multiresolution representation.As one can see, both the low-frequency scaling coefficients and high-frequency wavelet coefficients of a signal, at any resolution, can be decomposed, which leads to wavelet packet transform (WPT), and a much larger family of bases functions can be produced [33]. Moreover, 2D WPT can be obtained by using a tensor product of two 1D WPTs. Figure 1 shows a 2-level 2D WPT, where all the wavelet subbands LH1, HL1, and HH1, at resolution 1, are further decomposed into wavelet packet subbands delimited by dashed lines.


## 3. Wavelet-Packet-Tree-Based Texture Synthesis

Wavelet packet transform provides more bases functions, which leads to a more compact representation in comparison with wavelet transform. For the image-coding applications, we had proposed an efficient scheme to organize the wavelet packet (WP) coefficients of an image into hierarchical trees called WP trees [32]. In this section, we explore the key features of WP trees and propose a WP-tree-based algorithm for texture synthesis.


### 3.1. Wavelet Packet Tree

The WP coefficients of a sequence of wavelet coefficientsDℓ(n), at resolution ℓ, are computed by(3.1)D̃l,1(n)=∑kh(2n-k)⋅Dl(k),D̃l,2(n)=∑kg(2n-k)⋅Dl(k),
which can be efficiently rearranged and concatenated by(3.2)D̃l(n)={D̃l,1(n2);evenn,D̃l,2(n-12);oddn.
After the rearrangement and concatenation above, the dyadic relationship of wavelet coefficients across subbands is still valid and can be used to construct the dyadic WP trees of a signal. Similarly, the (2D) WP trees of an image can be obtained by rearranging the WP coefficients horizontally followed by vertically or vice versa. Figure 1 shows a 2-level WP tree with arrows connecting the related WP coefficients.The key idea behind the construction of WP trees is based on the spatial relationships of WP coefficients. It has the same structure, that is, quad-tree structure, as the conventional wavelet trees. Furthermore, the number of high-energy wavelet coefficients can be significantly reduced through the use of wavelet packet transform. Take the rice image shown in Figure5 as an example. The cumulative energy distribution (CED) of wavelet coefficients or WP coefficients is given by(3.3)CED(n)=∑i=1n|C(i)|∑i=1N|C(i)|×100%,
where |C(i)|,i=1,2,…,N, is the sorted magnitudes of wavelet coefficients or WP coefficients in descending order, n in is the number of coefficients, and N is the total number of coefficients. Figure 2 shows the CED curves of wavelet coefficients and WP coefficients, where the horizontal and vertical axes are the number of coefficients and percentage of energy, respectively. It is noted that the energy clustering of WP coefficients is more compact than wavelet coefficients. As a result, the WP-based representation is preferable to the wavelet-based representation for texture images.Figure 2
Cumulative energy distributions (CEDs) of the wavelet coefficients (solid line) and wavelet packet coefficients (dashed line) of the rice image shown in Figure5.


### 3.2. Proposed Algorithm

As noted, the low-frequency WP coefficients retain the global information of an image, and the high-frequency WP coefficients contain the local detail. It is desirable to coarsely synthesize an image based on the low-frequency WP coefficients and then tune the intermediate synthesis result based on the high-frequency WP coefficients. Motivated by the fact above, we propose an efficient WP-tree-based texture synthesis algorithm using a two-step process: a coarse searching followed by a fine tuning. Figure2 depicts a flowchart of the proposed algorithm. It is presented in steps as follows.Step 1 (initialization).
Decompose the source image by wavelet packet transform, rearrange the high-frequency WP coefficients, and construct the WP trees. Randomly select a WP tree from the source image, which is replicated in the upper left corner of the synthesis image.Step 2 (coarse matching).
For every WP tree to be synthesized, starting from the upper left corner to the lower right corner of the synthesis image, search the candidate WP trees from the source image by using a coarse matching with a tolerance as follows:(3.4)CWPj={WPsouce,i∣err⁡LFN(i,j)≤tolj,∀i},tolj=min⁡i⁡err⁡LFN(i,j)⋅(1+Tr),err⁡LFN(i,j)=∑p∈Ni,q∈Njdist⁡LFN(p,q),dist⁡LFN(p,q)=∑m,n∈LFN(WPsource,p(m,n)-WPsynthesis,q(m,n))2,
where CWPj is the set of candidate WP trees obtained from the source image for matching the synthesis WP tree WPsynthesis,j, with the tolerance tolj, Tr is a given threshold, distLFN(p,q) is the distance between WPsource,p and WPsynthesis,q based on the low-frequency nodes (LFNs) of WP trees, and Ni and Nj are neighbors of WPsource,i and WPsynthesis,j, respectively. Causal neighborhoods were used in our experiments.Step 3 (fine matching).
After the coarse matching in Step2, the following fine matching is used to find the best WP tree based on the high-frequency nodes (HFNs) of the candidate WP trees.                                    (3.5)WPsynthesis,j=WPsource,i;i=arg(min⁡i⁡err⁡HFN(i,j)),err⁡HFN(i,j)=∑p∈Ni,q∈Njdist⁡HFN(p,q),dist⁡HFN(p,q)=∑m,n∈HFN(WPsource,p(m,n)-WPsynthesis,q(m,n))2,
where WPsource,i∈CWPjin (3.5).Step 4.
Repeat Step2 followed by Step 3 until all the WP trees to be synthesized are obtained from the WP trees of the source image.Step 5.
Finally, the synthesis image is obtained by taking the inverse wavelet packet transform of the synthesis WP trees.To reduce the synthesis time, one can easily modify Steps2 and 3 by using patches of WP trees instead of single WP trees. Moreover, it is noted that lots of high-frequency WP coefficients are not significant, and only a small portion of middle-frequency WP coefficients are sufficient for the fine matching in Step 3. Thus, the proposed texture synthesis based on wavelet packet tree (TSWPT) algorithm is simple and computationally efficient. Flowchart of the TSWPT algorithm is shown in Figure 3.Figure 3
Flowchart of the TSWPT algorithm.


## 3.1. Wavelet Packet Tree

The WP coefficients of a sequence of wavelet coefficientsDℓ(n), at resolution ℓ, are computed by(3.1)D̃l,1(n)=∑kh(2n-k)⋅Dl(k),D̃l,2(n)=∑kg(2n-k)⋅Dl(k),
which can be efficiently rearranged and concatenated by(3.2)D̃l(n)={D̃l,1(n2);evenn,D̃l,2(n-12);oddn.
After the rearrangement and concatenation above, the dyadic relationship of wavelet coefficients across subbands is still valid and can be used to construct the dyadic WP trees of a signal. Similarly, the (2D) WP trees of an image can be obtained by rearranging the WP coefficients horizontally followed by vertically or vice versa. Figure 1 shows a 2-level WP tree with arrows connecting the related WP coefficients.The key idea behind the construction of WP trees is based on the spatial relationships of WP coefficients. It has the same structure, that is, quad-tree structure, as the conventional wavelet trees. Furthermore, the number of high-energy wavelet coefficients can be significantly reduced through the use of wavelet packet transform. Take the rice image shown in Figure5 as an example. The cumulative energy distribution (CED) of wavelet coefficients or WP coefficients is given by(3.3)CED(n)=∑i=1n|C(i)|∑i=1N|C(i)|×100%,
where |C(i)|,i=1,2,…,N, is the sorted magnitudes of wavelet coefficients or WP coefficients in descending order, n in is the number of coefficients, and N is the total number of coefficients. Figure 2 shows the CED curves of wavelet coefficients and WP coefficients, where the horizontal and vertical axes are the number of coefficients and percentage of energy, respectively. It is noted that the energy clustering of WP coefficients is more compact than wavelet coefficients. As a result, the WP-based representation is preferable to the wavelet-based representation for texture images.Figure 2
Cumulative energy distributions (CEDs) of the wavelet coefficients (solid line) and wavelet packet coefficients (dashed line) of the rice image shown in Figure5.


## 3.2. Proposed Algorithm

As noted, the low-frequency WP coefficients retain the global information of an image, and the high-frequency WP coefficients contain the local detail. It is desirable to coarsely synthesize an image based on the low-frequency WP coefficients and then tune the intermediate synthesis result based on the high-frequency WP coefficients. Motivated by the fact above, we propose an efficient WP-tree-based texture synthesis algorithm using a two-step process: a coarse searching followed by a fine tuning. Figure2 depicts a flowchart of the proposed algorithm. It is presented in steps as follows.Step 1 (initialization).
Decompose the source image by wavelet packet transform, rearrange the high-frequency WP coefficients, and construct the WP trees. Randomly select a WP tree from the source image, which is replicated in the upper left corner of the synthesis image.Step 2 (coarse matching).
For every WP tree to be synthesized, starting from the upper left corner to the lower right corner of the synthesis image, search the candidate WP trees from the source image by using a coarse matching with a tolerance as follows:(3.4)CWPj={WPsouce,i∣err⁡LFN(i,j)≤tolj,∀i},tolj=min⁡i⁡err⁡LFN(i,j)⋅(1+Tr),err⁡LFN(i,j)=∑p∈Ni,q∈Njdist⁡LFN(p,q),dist⁡LFN(p,q)=∑m,n∈LFN(WPsource,p(m,n)-WPsynthesis,q(m,n))2,
where CWPj is the set of candidate WP trees obtained from the source image for matching the synthesis WP tree WPsynthesis,j, with the tolerance tolj, Tr is a given threshold, distLFN(p,q) is the distance between WPsource,p and WPsynthesis,q based on the low-frequency nodes (LFNs) of WP trees, and Ni and Nj are neighbors of WPsource,i and WPsynthesis,j, respectively. Causal neighborhoods were used in our experiments.Step 3 (fine matching).
After the coarse matching in Step2, the following fine matching is used to find the best WP tree based on the high-frequency nodes (HFNs) of the candidate WP trees.                                    (3.5)WPsynthesis,j=WPsource,i;i=arg(min⁡i⁡err⁡HFN(i,j)),err⁡HFN(i,j)=∑p∈Ni,q∈Njdist⁡HFN(p,q),dist⁡HFN(p,q)=∑m,n∈HFN(WPsource,p(m,n)-WPsynthesis,q(m,n))2,
where WPsource,i∈CWPjin (3.5).Step 4.
Repeat Step2 followed by Step 3 until all the WP trees to be synthesized are obtained from the WP trees of the source image.Step 5.
Finally, the synthesis image is obtained by taking the inverse wavelet packet transform of the synthesis WP trees.To reduce the synthesis time, one can easily modify Steps2 and 3 by using patches of WP trees instead of single WP trees. Moreover, it is noted that lots of high-frequency WP coefficients are not significant, and only a small portion of middle-frequency WP coefficients are sufficient for the fine matching in Step 3. Thus, the proposed texture synthesis based on wavelet packet tree (TSWPT) algorithm is simple and computationally efficient. Flowchart of the TSWPT algorithm is shown in Figure 3.Figure 3
Flowchart of the TSWPT algorithm.


## 4. Experimental Results

In the first experiment, the size of source images is128×128, which are shown in the 1st column of Figure 4; the size of synthesis images is 256×256. The TSWPT algorithm is applied to patches of WP trees in order to reduce computation time. The size of patches is 11×11, and the width of overlapped regions between neighboring patches is 3. The root nodes of WP trees, that is, WP coefficients in the lowest frequency subband, are used in the coarse matching step to retain the global appearance of the input source image. Since the essential components of textures are mainly distributed in the middle-frequency subbands, the highest-frequency WP coefficients, that is, the leaf nodes of WP trees, can be ignored in the fine matching step. The biorthogonal 9/7 wavelet is used. The threshold value, Tr, is set to 0.1.Figure 4
Visual comparison. Source images (1st column), synthesis images using Efros’ algorithm [13] (2nd column), Cui’s algorithm [31] (3rd column), and the TSWPT algorithm (4th column).Figure 5
Synthesis images at double the size of the source image (top) using Efros’ algorithm [13] (middle) and the TSWPT algorithm (bottom).Two widely used algorithms, Efros’ algorithm [13] and Cui’s algorithm [31], were used for comparisons with the TSWPT algorithm. The synthesis results are shown in the 2nd, 3rd, and 4th columns of Figure 4, respectively. Visual inspection shows that TSWPT is comparable to Cui’s algorithm and is preferable to Efros’ algorithm. All the algorithms were implemented by Matlab without optimization in source codes. Table 1 shows the computation times running on PC with CPU of 1.7 GHz. It is noted that TSWPT is superior to both Efros’ and Cui’s algorithms.Table 1
Comparison of computation times.
Efros and freeman [13]Cui et al. [31]TSWPTImage 13.784 s0.703 s0.416 sImage 23.893 s0.854 s0.425 sImage 33.939 s0.810 s0.404 sImage 43.914 s0.737 s0.372 sImage 53.891 s0.731 s0.392 sIn the second experiment, the sizes of source image and synthesis image are192×128 and 384×256, respectively. All the settings are the same as the first experiment. The source image and synthesis results using Efros’ algorithm and TSWPT are shown in Figure 5. It is shown that there are some blocking effects in the synthesis image using Efros’ algorithm. However, blocking effects are likely to be eliminated by using the TSWPT algorithm due largely to the filtering operations of inverse wavelet packet transform in Step 5. Figure 6 shows the computation times at various enlargement rates of image size. As one can see, the TSWPT algorithm outperforms both Efros’ and Cui’s algorithms.Figure 6
Synthesis times of the texture image (Figure5) using Efros’ algorithm [13] (solid line), Cui’s algorithm [31] (dashed line), and the TSWPT algorithm (dotted line).


## 5. Conclusion

The multi-resolution approach is suitable for texture synthesis in terms of computation time. Wavelet packet transform provides more bases functions than wavelet transform and therefore produces a more compact representation. We adopt wavelet packet transform to analyze the input texture and organize wavelet packet coefficients to form hierarchical trees called wavelet packet trees. The low-frequency nodes of wavelet packet trees contain the global characteristics of an image; the high-frequency nodes contain the local details. Thus, we propose texture synthesis based on wavelet packet tree (TSWPT). It has the advantage of saving computation time dramatically, and moreover, no training process is needed. Given a128×128 texture, experimental results show that the computation time for synthesizing an 256×256 image is only a fraction of a second.



---

*Source: 305384-2012-05-17.xml*