Abstract
The representation of an image with several multiscale singular points has been the main concern in image processing. Based on the dual-tree complex wavelet transform (DT-CWT), a new image reconstruction (IR) algorithm from multiscale singular points is proposed. First, the image was transformed by DT-CWT, which provided multiresolution wavelet analysis. Then, accurate multiscale singular points for IR were detected in the DT-CWT domain due to the shift invariance and directional selectivity properties of DT-CWT. Finally, the images were reconstructed from the phases and magnitudes of the multiscale singular points by alternating orthogonal projections between the CT-DWT space and its affine space. Theoretical analysis and experimental results show that the proposed IR algorithm is feasible, efficient, and offers a certain degree of denoising. Furthermore, the proposed IR algorithm outperforms other classical IR algorithms in terms of performance metrics such as peak signal-to-noise ratio, mean squared error, and structural similarity.
1. Introduction
Since the human visual system processes images in a multiscale manner, representative features in the wavelet domain are a suitable choice for image reconstruction (IR) [1]. Among these representative features, singularities and irregular structures often carry the most important information in signals and images [2]. In images, intensity discontinuities provide the location of object contours, which is particularly relevant for recognition purposes [3]. Multiscale singular points (SPs), which represent unique structures in images, have many applications in image processing. For instance, the well-known scale-invariant feature transform offers one of the most popular methods for feature extraction [4]. The Harris-Laplace scale space has been explored for the detection of SPs based on the normalized Laplacian operator in the Gaussian scale space [5].
In a similar research on SPs in scale space, the problem of minimal representation has been thoroughly explored due to the compression made by Desolneux and Leclaire [6]. Another example can be found in the heat diffusion-based IR algorithm investigated by Elder and Zucker [7]. Lillhom and Nielsen [8] described an entropy-based variational reconstruction method. Inspired by Lillholm et al. [8, 9], Kanters et al. showed that images can be reconstructed based on the information content of SPs, with only a few hundred points in Gaussian scale space [10].
Hummel et al. formulated a new IR method from zero-crossings in the scale space [11]. Mallat and Zhong found an algorithm to reconstruct an approximate image based on local maxima of wavelet coefficients [12]. Moreover, Dalal and Triggs [13] and Allen and Kon [14] showed how to obtain approximate reconstructions under some additional hypotheses. However, fast wavelet transforms, such as the discrete wavelet transform (DWT), have poor directional selectivity, thus they are less suitable for estimating anisotropic properties. To obtain better directionality, an alternative approach is to implement approximate versions of continuous wavelets, such as Gabor wavelets. Unfortunately, doing so significantly increases the computational cost. In addition, the Mallat–Zhong algorithm has neither decomposition invariance nor the best directional selectivity.
In view of the aforementioned drawbacks of DWT in IR from SPs, this paper proposes a model for SP selection in the domain of dual-tree complex wavelet transform (DT-CWT). DT-CWT is widely used in image processing, such as image quality assessment [15], image denoising [16], image watermarking [17], and IR [18]. While retaining multiresolution characteristics and local time-frequency analysis capabilities, DT-CWT also has shift invariance and good directional selectivity [19]. Moreover, the high-frequency part of the two-dimensional (2D) DT-CWT reflecting detailed features can be described with ±15°, ±45°, and ±75° attributes in six directions [20]. As a result, DT-CWT can reflect image variations in different directions at different scales, thereby enabling better characterization of directional image features. Here, a new IR algorithm in the monogenic wavelet domain is also proposed [21], which accomplishes the extraction of multiscale SPs and an IR simultaneously.
To the best of our knowledge, few studies have been conducted on reconstruction algorithms in the DT-CWT domain [22]. Given this situation and the above discussion, a new image reconstruction algorithm is proposed here that involves less data and offers a certain degree of denoising. The main contributions of our proposed algorithm are as follows:(1)Multiscale SPs are obtained from DT-CWT wavelet coefficients, whereby accurate SPs can be extracted by combining shift invariance with directional selectivity.(2)Based on the corresponding magnitudes of the wavelet coefficients, the proposed SP selection model extracts only some of the multiscale SPs. This model reduces the amount of data required for IR and offers a certain degree of denoising.(3)Based on the magnitudes and phases of multiscale SPs, an interpolation-projection algorithm is proposed. The proposed reconstruction algorithm is more accurate and efficient compared to the Mallat–Hwang method [1] and the CWT method [12].
The rest of the paper is organized as follows. In Section 2, the principle of DT-CWT and the related IR algorithm are reviewed. Section 3 describes the proposed IR scheme in detail. Simulation details and experimental results are given in Section 4. Finally, a brief conclusion is drawn in Section 5.
2. Related Work
2.1. Dual-Tree Complex Wavelet Transform
Proposed by Kingsbury [23, 24], DT-CWT is an enhancement of DWT. Its filters employed in two trees are designed in such a way that the aliasing in one branch of the first tree is approximately cancelled by the corresponding branch in the second tree [25]. Kingsbury focused on designing a two-channel filter bank in which the filters in tree A satisfy the half-sample phase delay condition with respect to the filters in tree B, as shown in Figure 1, where filters have an approximate Hilbert transform relationship (90° out of phase with each other). The Hilbert transform leads to the CWT having the same advantages as the Fourier transform, such as (i) no oscillations, (ii) shift invariance, and (iii) being nonaliasing.
Directional selectivity was achieved in the 2D case by combining the outputs of the filter bank such that the equivalent complex filters have support in only one quadrant of the frequency plane [26]. The wavelet functions from the two trees play the role of the real and imaginary parts of a complex analytic wavelet. In this way, the 2D DT-CWT has six directions at approximately ±15°, ±45°, and ±75° [19].
With six 2D analytic wavelet sub-bands, the DT-CWT has six pairs of magnitude and phase information as features for one stage. The magnitude information indicates the discontinuities and mutability, which are used as crucial features, and for IR, the magnitude often carries important information. Phase is another important feature that provides structural information, such as edges, textures, and ridges. Furthermore, in typical natural images, the phase contains more structural information than the magnitude [27], and a rigid translation of image structures leads to a consistent phase shift [28]. The combination of magnitude and phase in the complex wavelet domain can represent image features better than a single magnitude in the discrete wavelet domain.
2.2. Statistical Modeling of the Magnitude of DT-CWT Wavelet Coefficients
Statistical modeling of wavelet coefficients is necessary for selecting SPs and is designed to eliminate those multiscale SPs with low magnitude. The applied model can select the fewest local maxima to obtain an approximate reconstructed image. These local maxima are selected using the modulus of the coefficients in the DT-CWT domain with respect to the neighborhood.
A Gaussian mixture model (GMM) is applied, which is a mixture of two Gaussian distributions [18]. One Gaussian distribution models the coefficients around zero, while the other accounts for the higher magnitude coefficients that constitute the singularities. Mixture Rayleigh distribution (MRD) was chosen to model the magnitude [29]; this involves a mixture of two Rayleigh distributions to obtain a better approximation than with a single Rayleigh distribution. However, both of the above-mentioned models have an excessive number of parameters. The inverse Gaussian distribution (IGD) [30] has fewer errors and parameters than the MRD and GMM, so an adaptive IGD model can be applied to fit the probability density function (PDF) of the DT-CWT coefficient magnitude based on the Kullback-Leibler divergence (KLD) [15]. The GMM, MRD, and IGD are expressed, respectively, aswhere – are the scale parameters, – are the weighting factors, is the sharpness parameter, and is the mean value of each high-frequency component , .
2.3. DWT-Based Related Reconstruction Algorithm
Image reconstruction from multiscale SPs has been studied in the wavelet coefficient maxima in the DWT domain. The maxima of the wavelet coefficient at multiscale contain important information about the image [12].
Let V be the space of all dyadic wavelet transforms of the functions in , where V begins from the zero-element. Let the set T be all of the specific functions. Such functions not only take the values of all multiscale maxima but must also be located at the corresponding positions in the scale space. Let a set , which obtains the elements of M that minimize the Sobolev norm. Alternating projections can be applied between V and T, which are orthogonal with respect to the Sobolev norm. As shown in Figure 2, an approximation of the DWT of signal f is reconstructed by alternating orthogonal projections between the affine space T and the space V; is close to (the original DWT of signal ) [1].
It is difficult to reconstruct the original image from the multiscale SPs of its DWT by a simple inverse transformation, but these SPs can be used to reconstruct approximate images via a deliberately designed iterative construction algorithm. In this case, the iteration is used to eliminate possible instabilities of the solution at a small additional computational cost.
3. Proposed Approach
3.1. Global Scheme
A -stage DT-CWT applied to decompose the original image requires four separable wavelet transforms in parallel [19]. As shown in Figure 3, each stage of the DT-CWT generates 4 × 4 sub-bands comprising 3 × 4 high-frequency ones and 1 × 4 low-frequency ones. The low-frequency sub-bands are used in the next stage of the DT-CWT analysis decomposition. In tree A, is applied along the rows and is applied along the columns, while in tree B, is applied along the rows and is applied along the columns, with similar deployment as in trees C and D.
The high-frequency wavelet coefficients and the low-frequency wavelet coefficients are obtained after the DT-CWT.where , and are the high-frequency and low-frequency coefficients of the real trees, respectively, and these coefficients are the addition or subtraction of the outputs of trees A and B in Figure 3. and are the high-frequency and low-frequency coefficients of the imaginary trees, respectively, and these coefficients are the addition or subtraction of the outputs of trees C and D. Note that is only calculated if .
3.2. Preprocessing
A thresholding preprocess is designed to remove near-zero coefficients while retaining as much image information as possible. Near-zero coefficients can lead to errors and make it difficult to model the histogram of wavelet-coefficient magnitudes at each scale .
The wavelet-coefficient magnitude at scale can be calculated from and :
All images in the USC-SIPI image database have been used [31]. The wavelet-coefficient magnitude threshold increases as the decomposition stage increases; for instance, the thresholds are ∼, ∼, and ∼ for decomposition stages 1, 2, and 3, respectively. After selecting the threshold, the magnitude histograms are ready to fit a particular distribution model.
3.3. Multiscale Singular Points Selection
It is meaningful to choose the appropriate model to fit the magnitudes of the DT-CWT coefficients. The GMM, MRD, and IGD models were compared with the KLD values of the original distribution and the fitted models. All images in the USC-SIPI image database were tested. As shown in Figure 4, comparing these three candidate models involves different PDF fitting curves. The average KLD values were ∼0.05 for the GMM model, ∼0.04 for the MRD model, and ∼0.02 for the IGD model; accordingly, we chose the IGD model with the smallest KLD value. Using the fitted IGD model, a threshold can be set to filter the DT-CWT wavelet coefficients in each high-frequency component. Based on Equation (3), in particular, can be defined aswhere is the integral of the probability density integrated from 0 to . Through experiments, an optimal value of was chosen from 0.85 to 0.95 in this paper. The threshold value is used to control the number of SPs in the reconstructed image.
(a)
(b)
(c)
(d)
The multiscale SP matrix at scale is defined as the local maxima in .where finds a local maximum and is a threshold.
The phase at the multiscale SP can be calculated as [25]
3.4. Proposed Image Reconstruction Scheme
Let V be the set of a J-stage DT-CWT of an image and let the set T be all of the specific functions. Such functions must not only satisfy the magnitude and phase of all multiscale SPs at scale but must also be located at the corresponding positions in the scale space. By letting the set and the alternate projection begin from the zero-element of V, the algorithm converges strongly to the element of M whose Sobolev norm is minimal.
The proposed IR algorithm is shown in Figure 5, and the details of the reconstruction process are described below. The differences between the Mallat–Hwang algorithm and the proposed reconstruction algorithm are similar, except that the latter converges to an approximate reconstructed image within a few iterations.
Step 1. The updated multiscale wavelet-coefficient matrix is obtained from the proposed interpolation algorithm based on and the magnitude-modifying algorithm.
Let us define a constant based on the scale index aswhere an optimal value of from 0.75 to 0.8 can be chosen experimentally.
For example, the real part of the interpolation value between two maxima coordinates and on the -axis can be expressed aswhere t denotes the times of alternate projections.
Moreover, the real part of the interpolation value between two maxima coordinates and on the -axis can be expressed asThe interpolation algorithm for the imaginary part can be obtained similarly based on Equations (10) and (11).
The updated multiscale wavelet-coefficient matrix can be obtained after interpolation and modification of the magnitudes with the phase as follows.
When the interpolation magnitudes are modified on the -axis,Subsequently, when the interpolation magnitudes are modified on the -axis,After completing the interpolation and modification process, the updated multiscale wavelet-coefficient matrix is obtained. This step represents the projection process from V to T.
Step 2. By combining the low-frequency data and the high-frequency data , the reconstructed image is obtained with the aid of the inverse DT-CWT.
Step 3. After placing the reconstructed image into the structure of DT-CWT for decomposition, the new high-frequency data can be obtained as the projection from T to V.
Step 4. The previous steps are repeated until the peak signal-to-noise ratio (PSNR) of no longer increases. The DT-CWT of converges to the element of the set M.
The PSNR is defined aswhere M × N is the size of the original image.
4. Simulation Results and Analysis
4.1. Parameters Setup
To analyze the proposed IR algorithm, experiments were conducted on a computer (1.80 GHz, i5-8250U CPU, 16.00 GB RAM) using MATLAB 2016b. Selected from the USC-SIPI image database [26], the test images were the grayscale “House” image of size 256 × 256 [as shown in Figure 6a] and the grayscale “Boat,” “Elaine,” “Barbara” and “Peppers” of size 512 × 512 (as shown in the first row of Figure 7). To measure the quality of the reconstructed images, as well as the PSNR defined in Equation (15), two other image quality assessment criteria were considered, namely.(1)The mean squared error (MSE) is defined as where denotes the original image and denotes the reconstructed image.(2)The structural similarity (SSIM) is defined as where is the mean value of , is the mean value of , is the variance of , and is the variance of corresponding to the special case of [32]. The greater the value of SSIM, the better the quality of the reconstructed image.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
The optimized scale for the DT-CWT analysis and synthesis process is selected from scale . The relationship between the number of iterations and the PSNR of the reconstructed image is shown in Figure 8, and the relationship between the scale index j and the PSNR of the reconstructed image is shown in Figure 9. Comparing the PSNR values at various points in Figures 8 and 9 gives the optimal scale index J and the number of iterations. The PSNR of the reconstructed image decreases at J ≧ 3 due to the loss of data after multiple DT-CWT decompositions. As shown in Figures 10 and 11, the SSIM values show consistent results with the PSNR values relative to the number of iterations. After five or six iterations, the PSNR and SSIM values do not increase significantly, so the iterative process is stopped to make a trade-off between model performance and computational cost. In conclusion, it is possible to choose the optimal scale index J = 3 and the optimal number of iterations of 5.
4.2. Analysis of Image Reconstruction Schemes
How to code images with a minimum number of bits for transmission or storage is always a challenge and requires a large number of simulations to evaluate the proposed reconstruction algorithm. Figure 6 shows the SPs of the “House” image at scale in the ±15°, ±45°, and ±75° directions, along with the reconstructed image. The low-frequency images used for reconstruction are shown in the second row of Figure 12, which were reconstructed from the low-frequency band ; the reconstructed images from the low-frequency band and the multiscale SPs are shown in the third row. Comparing the images in Figures 6(b)–6(g) shows that the SPs in each high-frequency component retain some unique features of the image in different directions. Due to the loss of some high-frequency coefficients, the reconstructed images appear smoother than the original images. However, these reconstructed images retain the local features of the original image.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
Table 1 gives the results of the image quality assessment of the three test images (“Boat,” “Elaine,” and “Peppers”) when using the methods of Mallat and Zhong [12] and Hua and Orchard [33] and the proposed method. The PSNR values for all five reconstructed images are plotted in Figure 7. The SSIM values for all five reconstructed images are plotted in Figure 13. Under the same experimental conditions, the proposed algorithm gives the best reconstruction quality than the other two algorithms [12, 33].
As can be seen in Table 1, the proposed method achieves good results in terms of performance metrics such as SSIM and PSNR with fewer SPs compared to other state-of-the-art methods. This implies that the proposed method properly reconstructs the detailed information of the source image by extracting accurate SPs with the properties of shift invariance and good directional selectivity.
4.3. Data Requirements and Efficiency Analysis of the Reconstruction Scheme
The execution time is an important consideration in IR. The running times of the proposed reconstruction algorithm and similar ones [12, 33] are given in Table 2. Under the same conditions, the percentage of pixels used (including SPs and data points in lower frequency bands) in the original image shows that the proposed reconstruction algorithm is more efficient, while achieving higher PSNR values. In the proposed scheme, the time-consuming processes are (i) coefficients fitting, (ii) SPs selection, and (iii) interpolation-projection algorithm. Figure 14 shows the reconstruction time required for each part, which indicates that the proposed IR algorithm is feasible and efficient.
Compared to other similar reconstruction algorithms, the proposed IR algorithm is capable of obtaining a satisfactory approximate solution based on more representative SPs in fewer iterations.
4.4. Denoising Effect with the Proposed Algorithm
The shift invariance and the directionality of DT-CWT make it advantageous in many areas of image processing, such as denoising [16]. With the automatic threshold to control the number of retained coefficients in the wavelet domain, the proposed reconstruction algorithm offers a certain degree of denoising. Gaussian noise is added to the original image to test the denoising effect as follows:where C’ and C denote the original image with and without noise, respectively, G indicates white Gaussian noise with zero mean and unit standard deviation, and k is the coefficient of noise intensity. The results of the reconstructed image “Elaine” with white Gaussian noise added to the original image are shown in the first row of Figure 15 for k = 0.01, 0.1, 0.2, and 0.5, respectively. Although the quality of the reconstructed images decreases with increasing noise parameters, the reconstructed images are still identifiable, as shown in the second row of Figure 15. Furthermore, the proposed algorithm retains the sharpness of the singularities from polluted images, but smooths the reconstructed images in the spatial domain. Table 3 gives the PSNR and SSIM values for different noisy original images, where the PSNR values can reflect the noise intensity, and Table 4 gives the PSNR values for different reconstructed images with increasing noise parameters. Moreover, Table 5 gives the SSIM values of different reconstructed images as the noise intensity increases.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
As can be seen in Tables 3–5, the performance of the proposed image reconstruction algorithm outperforms the other algorithms. Specifically, the reconstructed image of the proposed scheme achieves the best denoising results in terms of both PNSR and SSIM because the multiscale points of interest in the CT-DWT domain retain the relevant image singularity information, while other classical filtering methods, such as mean filtering [34], median filtering [35], and Wiener filtering [36], have the ability to filter noise but do not have the ability to preserve the singularity in image. Therefore, the proposed image reconstruction algorithm is good at retaining finer image details while removing noise to a certain extent. Indeed, there are many advanced image denoising methods that outperform our proposed scheme. It should be emphasized that our proposed algorithm is mainly intended for image reconstruction from SPs in the DT-CWT domain rather than image denoising.
5. Conclusion
A multiscale singular point-based image reconstruction is proposed for the first time in the DT-CWT domain. DT-CWT offers properties of shift invariance and more choices of directional selectivity, enabling it to better characterize multiscale singularities. An automatic singular point selection model was built based on the inverse Gaussian distribution model. By combining the wavelet-coefficient model fitting and iterative projection operations, IR is completed with the optimally selected multiscale singular points. Theoretical analysis and experimental results show that the proposed image reconstruction scheme not only reduces the amount of information involved but also has a certain degree of denoising effect. Moreover, the proposed IR algorithm outperforms some other algorithms in terms of performance metrics, such as PSNR, MSE, and SSIM. Future work may involve exploring the reconstruction of color images with multiscale singular points in the DT-CWT domain.
Data Availability
The databases of tested images were taken from the following database USC-SIPI, Signal and Image Processing Institute, The University of Southern California, http://sipi.usc.edu/database/
Conflicts of Interest
The authors declare no conflicts of interest.
Acknowledgments
This work was partly supported by the National Natural Science Foundation of China (Grant nos. 62041106 and 61662047).