Image Reconstruction from Multiscale Singular Points Based on the Dual-Tree Complex Wavelet Transform

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 aﬃne space. Theoretical analysis and experimental results show that the proposed IR algorithm is feasible, eﬃcient, and oﬀers 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.


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 scaleinvariant feature transform offers one of the most popular methods for feature extraction [4].
e 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 zerocrossings 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.e 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. is 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.e proposed reconstruction algorithm is more accurate and efficient compared to the Mallat-Hwang method [1] and the CWT method [12].e 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.

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 halfsample phase delay condition with respect to the filters in tree B, as shown in Figure 1, where filters h i (n) and g i (n), (i � 0, 1) have an approximate Hilbert transform relationship (90 °out of phase with each other).
e 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].e 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.
e 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].e combination of magnitude and phase in the complex wavelet domain can represent image features better than a single magnitude in the discrete wavelet domain.

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.e applied model can select the fewest local maxima to obtain an approximate reconstructed image.ese 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.
e 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].e GMM, MRD, and IGD are expressed, respectively, as 2 Security and Communication Networks where σ 1 -σ 4 are the scale parameters, k 1 -k 4 are the weighting factors, λ i is the sharpness parameter, and μ i is the mean value of each high-frequency component i, i � 1, 2, 3, . . ., i ∈ Z.

DWT-Based Related Reconstruction Algorithm.
Image reconstruction from multiscale SPs has been studied in the wavelet coefficient maxima in the DWT domain.e 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 L 2 (R 2 ), 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 M � V ∩ T, 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 Ws of the DWTof signal f is reconstructed by alternating orthogonal projections between the affine space T and the space V; Ws is close to Wf (the original DWT of signal f) [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.

Global Scheme.
A J-stage DT-CWT applied to decompose the original image f(x, y) 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.e low-frequency sub-bands are used in the next stage of the DT-CWTanalysis decomposition.In tree A, h 0 (x), h 1 (x)   is applied along the rows and h 0 (y), h 1 (y)   is applied along the columns, while in tree B, g 0 (x), g 1 (x)   is applied along the rows and g 0 (y), g 1 (y)   is applied along the columns, with similar deployment as in trees C and D.
e high-frequency wavelet coefficients H j [f(x, y)] and the low-frequency wavelet coefficients L J [f(x, y)] are obtained after the DT-CWT. where ] are the highfrequency 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.
] 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 Figure 1: Analysis filter banks for the dual-tree complex wavelet transform [19].

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 2 j .e wavelet-coefficient magnitude |H j [f(x, y)]| at scale 2 j (j � 1, 2, . . ., J) can be calculated from All images in the USC-SIPI image database have been used [31].e wavelet-coefficient magnitude threshold increases as the decomposition stage increases; for instance, the thresholds are ∼10 − 3 , ∼10 − 2 , and ∼10 − 1 for decomposition stages 1, 2, and 3, respectively.After selecting the threshold, the magnitude histograms are ready to fit a particular distribution model.

Multiscale Singular Points Selection.
It is meaningful to choose the appropriate model to fit the magnitudes of the DT-CWT coefficients.e 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.e 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 σ i can be set to filter the DT-CWT wavelet coefficients in each high-frequency component.Based on Equation (3), in particular, σ i can be defined as where 0 < ξ < 1 is the integral of the probability density integrated from 0 to σ i .rough experiments, an optimal value of ξ was chosen from 0.85 to 0.95 in this paper.e threshold value σ i is used to control the number of SPs in the reconstructed image.e multiscale SP matrix S j (x, y) at scale 2 j , (j � 1, 2, . . ., J) is defined as the local maxima in where arglmax finds a local maximum and σ i is a threshold.e phase θ at the multiscale SP can be calculated as [25] θ 4 Security and Communication Networks multiscale SPs S j (x, y) at scale 2 j , (j � 1, 2, . . ., J) but must also be located at the corresponding positions in the scale space.By letting the set M � V ∩ T 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.e proposed IR algorithm is shown in Figure 5, and the details of the reconstruction process are described below.
e 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.
e updated multiscale wavelet-coefficient matrix S j ′ (x, y) is obtained from the proposed interpolation algorithm based on S j (x 1 , y) and the magnitude-modifying algorithm.
Let us define a constant c 0 based on the scale index j as where an optimal value of c from 0.75 to 0.8 can be chosen experimentally.
For example, the real part of the interpolation value V t interpo (x 0 , y) between two maxima coordinates (x 1 , y) and (x 2 , y) on the x-axis can be expressed as Security and Communication Networks where t denotes the times of alternate projections.Moreover, the real part of the interpolation value V t interpo (x, y 0 ) between two maxima coordinates (x, y 1 ) and (x, y 2 ) on the y-axis can be expressed as e interpolation algorithm for the imaginary part can be obtained similarly based on Equations (10) and (11).
e updated multiscale wavelet-coefficient matrix S j ′ (x, y) can be obtained after interpolation and modification of the magnitudes with the phase θ S j (x,y) as follows.
When the interpolation magnitudes are modified on the x-axis, After completing the interpolation and modification process, the updated multiscale wavelet-coefficient matrix S j ′ (x, y) is obtained.
is step represents the projection process from V to T.
Step 2. By combining the low-frequency data L J [f(x, y)] and the high-frequency data S j ′ (x, y), the reconstructed image f ′ (x, y) is obtained with the aid of the inverse DT-CWT.
Step 3.After placing the reconstructed image f ′ (x, y) into the structure of DT-CWT for decomposition, the new highfrequency data H j ′ (x, y) can be obtained as the projection from T to V.
e previous steps are repeated until the peak signalto-noise ratio (PSNR) of f ′ (x, y) no longer increases.e DT-CWT of f ′ (x, y) converges to the element of the set M.
e PSNR is defined as where M × N is the size of the original image.

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) e mean squared error (MSE) is defined as where f(x, y) denotes the original image and f ′ (x, y) denotes the reconstructed image.
(2) e structural similarity (SSIM) is defined as where μ f is the mean value of f, μ f′ is the mean value of f ′ , σ 2 f is the variance of f, and σ 2 f is the variance of f ′ corresponding to the special case of C 1 � C 2 � 0 [32].e greater the value of SSIM, the better the quality of the reconstructed image.e optimized scale for the DT-CWT analysis and synthesis process is selected from scale 2 J .e 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.e 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.

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 2 2 in the ±15 °, ±45 °, and ±75 °directions, along with the reconstructed image.
e low-frequency images used for reconstruction are shown in the second row of Figure 12, which were reconstructed from the low-frequency band L 3 [f(x, y)]; 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.

Security and Communication Networks
However, these reconstructed images retain the local features of the original image.
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.
e PSNR values for all five reconstructed images are plotted in Figure 7.
e 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. is 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."House" "Boat" "Elaine" "Peppers" "Barbara"

Data Requirements and Efficiency Analysis of the Reconstruction Scheme.
e execution time is an important consideration in IR. e 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.

Denoising Effect with the Proposed Algorithm.
e 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: Security and Communication Networks 9 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.e 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  "Barbara" "House" Figure 13: SSIM values with different methods: the five points on the abscissa denote the reconstructed images of "House," "Boat," "Elaine," "Peppers," and "Barbara."proposed algorithm retains the sharpness of the singularities from polluted images, but smooths the reconstructed images in the spatial domain.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.erefore, 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.

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.eoretical 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.

Figure 4 :
Figure 4: Comparison of three probability density function (PDF) candidates: (a) original image; (b-d) curve fitting for the inverse Gaussian distribution (IGD), Gaussian mixture model (GMM), and mixture Rayleigh distribution (MRD), respectively.e histograms are taken from the image "Peppers," with a scale index of 3 and an orientation of +75 °.

Figure 8 :
Figure 8: Relationship between the peak signal-to-noise ratio (PSNR) of the reconstructed image and the number of iterations for each image for the scale index J � 3.

Figure 9 :Figure 10 :Figure 11 :
Figure 9: Relationship between PSNR of the reconstructed image and the scale index j for five iterations per image.

Table 2 :
Analysis of the IR scheme with scale index J � 3 and five iterations.
Figure14: Time required for and proportion of each procedure of "Boat" reconstruction.e time required for the five iterations of IR is 5.901 s for the maximum scale index J � 3.

Table 3
gives the PSNR and SSIM values for different noisy original images, where the PSNR values can reflect the noise intensity, and Table4gives the PSNR values for different reconstructed images with increasing noise parameters.Moreover, Table5gives the SSIM values of different reconstructed images as the noise intensity increases.

Table 3 :
PSNR (dB) and SSIM of original images with various intensities of white Gaussian noise.

Table 4 :
PSNR (dB) of reconstructed images with various intensities of white Gaussian noise added to the original image (values in brackets are PSNR gains relative to PSNRs of the original image).

Table 5 :
SSIM of reconstructed images with various intensities of white Gaussian noise added to the original image (values in brackets are SSIM gains relative to SSIMs of the original image).