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

Authors: Sihang Liu; Benoit Tremblais; Phillippe Carre; Nanrun Zhou; Jianhua Wu
Journal: Security and Communication Networks
 (2021)

---

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


---

## Body

## 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 Section2, 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 hinandgin,i=0,1 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.Figure 1
Analysis filter banks for the dual-tree complex wavelet transform [19].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, as(1)GMMx=k112πσ1exp−x22σ12+k212πσ2exp−x22σ22,(2)MRDx=k3xσ32exp−x22σ32+k4xσ42exp−x22σ42,(3)IGDx=λi2πx3exp−λix−μi22μi2x,where σ1–σ4 are the scale parameters, k1–k4 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.


### 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].LetV be the space of all dyadic wavelet transforms of the functions in L2ℝ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 DWT of 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].Figure 2
Alternating orthogonal projections beginning from the zero-element and converging toM=V∩T [12].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.


## 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 hinandgin,i=0,1 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.Figure 1
Analysis filter banks for the dual-tree complex wavelet transform [19].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, as(1)GMMx=k112πσ1exp−x22σ12+k212πσ2exp−x22σ22,(2)MRDx=k3xσ32exp−x22σ32+k4xσ42exp−x22σ42,(3)IGDx=λi2πx3exp−λix−μi22μi2x,where σ1–σ4 are the scale parameters, k1–k4 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.


## 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].LetV be the space of all dyadic wavelet transforms of the functions in L2ℝ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 DWT of 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].Figure 2
Alternating orthogonal projections beginning from the zero-element and converging toM=V∩T [12].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

AJ-stage DT-CWT applied to decompose the original image fx,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. The low-frequency sub-bands are used in the next stage of the DT-CWT analysis decomposition. In tree A, h0x,h1x is applied along the rows and h0y,h1y is applied along the columns, while in tree B, g0x,g1x is applied along the rows and g0y,g1y is applied along the columns, with similar deployment as in trees C and D.Figure 3
DT-CWT decomposition into four wavelet trees in parallel [26].The high-frequency wavelet coefficientsHjfx,y and the low-frequency wavelet coefficients LJfx,y are obtained after the DT-CWT.(4)Hjfx,y=HjRefx,y+&ImaginaryI;HjImfx,y,j=1,2,…,J,LJfx,y=LJRefx,y+&ImaginaryI;LJImfx,y,where i=−1, HjRefx,y and LJRefx,y 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. HjImfx,y and LJImfx,y 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 LJfx,y is only calculated if j=J.


### 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 scale2j.The wavelet-coefficient magnitudeHjfx,y at scale 2jj=1,2,…,J can be calculated from HjRefx,y and HjImfx,y:(5)Hjfx,y=HjRefx,y2+HjImfx,y2.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 ∼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.


### 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 Figure4, 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 σ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(6)Fσi=∫0σiIGDxdx=ξ,where 0<ξ<1 is the integral of the probability density integrated from 0 to σi. Through experiments, an optimal value of ξ was chosen from 0.85 to 0.95 in this paper. The threshold value σi is used to control the number of SPs in the reconstructed image.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. The histograms are taken from the image “Peppers,” with a scale index of 3 and an orientation of +75°.
(a)(b)(c)(d)The multiscale SP matrixSjx,y at scale 2j,j=1,2,…,J is defined as the local maxima in Hjfx,y.(7)Sjx,y=arglmaxx,yHjfx,y|Hjfx,y≥σi,where arglmax finds a local maximum and σi is a threshold.The phaseθ at the multiscale SP can be calculated as [25](8)θSjx,y=tan−1HjImfx,yHjRefx,y.


### 3.4. Proposed Image Reconstruction Scheme

LetV 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 Hjfx,y and phase θSjx,y of all multiscale SPs Sjx,y at scale 2j,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.The proposed IR algorithm is shown in Figure5, 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.Figure 5
Outline of the proposed algorithm.Step 1.
The updated multiscale wavelet-coefficient matrixSj′x,y is obtained from the proposed interpolation algorithm based on Sjx1,y and the magnitude-modifying algorithm.
Let us define a constantc0 based on the scale index j as(9)c0=c1−j,where an optimal value of c from 0.75 to 0.8 can be chosen experimentally.
For example, the real part of the interpolation valueVinterpotx0,y between two maxima coordinates x1,y and x2,y on the x-axis can be expressed as(10)ReVinterpotx0,y=c0x0−x11−c02x2−x01−c04ReSjx1,y−Vtx1,y+c0x2−x01−c02x0−x11−c04ReSjx2,y−Vtx2,y,where t denotes the times of alternate projections.
Moreover, the real part of the interpolation valueVinterpotx,y0 between two maxima coordinates x,y1 and x,y2 on the y-axis can be expressed as(11)ReVinterpotx,y0=c0y0−y11−c02y2−y01−c04ReSjx,y1−Vtx,y1+c0y2−y01−c02y0−y11−c04ReSjx,y2−Vtx,y2.
The interpolation algorithm for the imaginary part can be obtained similarly based on Equations (10) and (11).
The updated multiscale wavelet-coefficient matrixSj′x,y can be obtained after interpolation and modification of the magnitudes with the phase θSjx,y as follows.
When the interpolation magnitudes are modified on thex-axis,(12)Sj′x,y=Sj′x0−1,y,ifθSj′x0,yπ−0.5<0.126andSj′x0,y>Sj′x0−1,y,Sjj′x0+1,y,ifθSj′x0,yπ−0.5<0.126andSj′x0,y>Sj′x0+1,y.
Subsequently, when the interpolation magnitudes are modified on they-axis,(13)Sj′x,y0=Sj′x,y0−1,ifθSj′x,y0π−0.5>0.376andSj′x,y0>Sj′x,y0−1,Sj′x,y0+1,ifθSj′x,y0π−0.5>0.376andSj′x,y0>Sjj′x,y0+1.
After completing the interpolation and modification process, the updated multiscale wavelet-coefficient matrixSj′x,y is obtained. This step represents the projection process from V to T.Step 2.
By combining the low-frequency dataLJfx,y and the high-frequency data Sj′x,y, the reconstructed image f′x,y is obtained with the aid of the inverse DT-CWT.Step 3.
After placing the reconstructed imagef′x,y into the structure of DT-CWT for decomposition, the new high-frequency data Hj′x,y 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) off′x,y no longer increases. The DT-CWT of f′x,y converges to the element of the set M.
The PSNR is defined as(14)PSNR=10log10M×N×2552∑x=1M∑y=1Nfx,y−f′x,y2,where M × N is the size of the original image.


## 3.1. Global Scheme

AJ-stage DT-CWT applied to decompose the original image fx,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. The low-frequency sub-bands are used in the next stage of the DT-CWT analysis decomposition. In tree A, h0x,h1x is applied along the rows and h0y,h1y is applied along the columns, while in tree B, g0x,g1x is applied along the rows and g0y,g1y is applied along the columns, with similar deployment as in trees C and D.Figure 3
DT-CWT decomposition into four wavelet trees in parallel [26].The high-frequency wavelet coefficientsHjfx,y and the low-frequency wavelet coefficients LJfx,y are obtained after the DT-CWT.(4)Hjfx,y=HjRefx,y+&ImaginaryI;HjImfx,y,j=1,2,…,J,LJfx,y=LJRefx,y+&ImaginaryI;LJImfx,y,where i=−1, HjRefx,y and LJRefx,y 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. HjImfx,y and LJImfx,y 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 LJfx,y is only calculated if j=J.


## 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 scale2j.The wavelet-coefficient magnitudeHjfx,y at scale 2jj=1,2,…,J can be calculated from HjRefx,y and HjImfx,y:(5)Hjfx,y=HjRefx,y2+HjImfx,y2.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 ∼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.


## 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 Figure4, 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 σ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(6)Fσi=∫0σiIGDxdx=ξ,where 0<ξ<1 is the integral of the probability density integrated from 0 to σi. Through experiments, an optimal value of ξ was chosen from 0.85 to 0.95 in this paper. The threshold value σi is used to control the number of SPs in the reconstructed image.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. The histograms are taken from the image “Peppers,” with a scale index of 3 and an orientation of +75°.
(a)(b)(c)(d)The multiscale SP matrixSjx,y at scale 2j,j=1,2,…,J is defined as the local maxima in Hjfx,y.(7)Sjx,y=arglmaxx,yHjfx,y|Hjfx,y≥σi,where arglmax finds a local maximum and σi is a threshold.The phaseθ at the multiscale SP can be calculated as [25](8)θSjx,y=tan−1HjImfx,yHjRefx,y.


## 3.4. Proposed Image Reconstruction Scheme

LetV 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 Hjfx,y and phase θSjx,y of all multiscale SPs Sjx,y at scale 2j,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.The proposed IR algorithm is shown in Figure5, 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.Figure 5
Outline of the proposed algorithm.Step 1.
The updated multiscale wavelet-coefficient matrixSj′x,y is obtained from the proposed interpolation algorithm based on Sjx1,y and the magnitude-modifying algorithm.
Let us define a constantc0 based on the scale index j as(9)c0=c1−j,where an optimal value of c from 0.75 to 0.8 can be chosen experimentally.
For example, the real part of the interpolation valueVinterpotx0,y between two maxima coordinates x1,y and x2,y on the x-axis can be expressed as(10)ReVinterpotx0,y=c0x0−x11−c02x2−x01−c04ReSjx1,y−Vtx1,y+c0x2−x01−c02x0−x11−c04ReSjx2,y−Vtx2,y,where t denotes the times of alternate projections.
Moreover, the real part of the interpolation valueVinterpotx,y0 between two maxima coordinates x,y1 and x,y2 on the y-axis can be expressed as(11)ReVinterpotx,y0=c0y0−y11−c02y2−y01−c04ReSjx,y1−Vtx,y1+c0y2−y01−c02y0−y11−c04ReSjx,y2−Vtx,y2.
The interpolation algorithm for the imaginary part can be obtained similarly based on Equations (10) and (11).
The updated multiscale wavelet-coefficient matrixSj′x,y can be obtained after interpolation and modification of the magnitudes with the phase θSjx,y as follows.
When the interpolation magnitudes are modified on thex-axis,(12)Sj′x,y=Sj′x0−1,y,ifθSj′x0,yπ−0.5<0.126andSj′x0,y>Sj′x0−1,y,Sjj′x0+1,y,ifθSj′x0,yπ−0.5<0.126andSj′x0,y>Sj′x0+1,y.
Subsequently, when the interpolation magnitudes are modified on they-axis,(13)Sj′x,y0=Sj′x,y0−1,ifθSj′x,y0π−0.5>0.376andSj′x,y0>Sj′x,y0−1,Sj′x,y0+1,ifθSj′x,y0π−0.5>0.376andSj′x,y0>Sjj′x,y0+1.
After completing the interpolation and modification process, the updated multiscale wavelet-coefficient matrixSj′x,y is obtained. This step represents the projection process from V to T.Step 2.
By combining the low-frequency dataLJfx,y and the high-frequency data Sj′x,y, the reconstructed image f′x,y is obtained with the aid of the inverse DT-CWT.Step 3.
After placing the reconstructed imagef′x,y into the structure of DT-CWT for decomposition, the new high-frequency data Hj′x,y 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) off′x,y no longer increases. The DT-CWT of f′x,y converges to the element of the set M.
The PSNR is defined as(14)PSNR=10log10M×N×2552∑x=1M∑y=1Nfx,y−f′x,y2,where 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(15)MSE=1MN∑x=1M∑y=1Nfx,y−f′x,y2,wherefx,y denotes the original image and f′x,y denotes the reconstructed image.(2)
The structural similarity (SSIM) is defined as(16)SSIMf,f′=2μfμf′+C12σff′+C2μf2+μf′2+C1σf2+σf′2+C2,whereμf is the mean value of f, μf′ is the mean value of f′, σf2 is the variance of f, and σf2 is the variance of f′ corresponding to the special case of C1=C2=0 [32]. The greater the value of SSIM, the better the quality of the reconstructed image.Figure 6
(a) Original “House” image; (b–g) SPs in the directions of +15°, +45°, and +75°, -−15°, −45°, and −75° at scale22; (h) reconstructed image.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 7
PSNR values with different methods: the five points on the abscissa denote the reconstructed images of “House,” “Boat,” “Elaine,” “Peppers,” and “Barbara.”The optimized scale for the DT-CWT analysis and synthesis process is selected from scale2J. 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.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 indexJ = 3.Figure 9
Relationship between PSNR of the reconstructed image and the scale indexj for five iterations per image.Figure 10
Relationship between the peak structural similarity (SSIM) of the reconstructed image and the number of iterations for each image for the scale indexJ = 3.Figure 11
Relationship between SSIM of the reconstructed image and the scale indexj for five iterations per image.


### 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. Figure6 shows the SPs of the “House” image at scale 22 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 L3fx,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. However, these reconstructed images retain the local features of the original image.Figure 12
First row: original images: (a) “Boat,” (b) “Elaine,” (c) “Barbara,” and (d) “Peppers.” Second row: low-frequency images subsampled by a factor of 23 along their rows and columns. Third row: reconstructed images.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Table1 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].Table 1
Structural similarity (SSIM), mean squared error (MSE), and PSNR (dB) values of the reconstructed images.
SchemeImageNo. Of SPsRatio of SPs over original image (%)SSIMMSEPSNRProposed“Boat”197549.100.8383104.9727.92“Elaine”149447.260.859354.2130.79“Peppers”189218.780.864288.5228.66Ref. [12]“Boat”195379.020.8123134.4526.72“Elaine”151227.330.8201117.8627.24“Peppers”190328.820.8165122.3427.12Ref. [33]“Boat”201319.240.8013198.7825.13“Elaine”150427.300.8196163.1326.05“Peppers”191248.850.8125186.9325.27The bold values show better results than the other methods.Figure 13
SSIM values with different methods: the five points on the abscissa denote the reconstructed images of “House,” “Boat,” “Elaine,” “Peppers,” and “Barbara.”As can be seen in Table1, 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.Table 2
Analysis of the IR scheme with scale indexJ=3 and five iterations.
SchemeImageTime (s)PSNR (dB)SSIMProposed algorithm“Boat”5.90127.920.8383“Elaine”5.37930.790.8593“Peppers”5.63428.660.8642Ref. [12]“Boat”7.55626.720.8123“Elaine”7.23227.240.8201“Peppers”7.44127.120.8165Ref. [33]“Boat”6.02425.130.8013“Elaine”5.85726.050.8196“Peppers”5.93225.270.8125Figure 14
Time required for and proportion of each procedure of “Boat” reconstruction. The time required for the five iterations of IR is 5.901 s for the maximum scale indexJ = 3.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:(17)C′=C+kG,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.Figure 15
Comparison of denoising performances with various intensities of white Gaussian noise. First row: original image with added white Gaussian noise,k = (a) 0.01, (b) 0.1, (c) 0.2, and (d) 0.5. Second row: reconstructed image with k = (e) 0.01, (f) 0.1, (g) 0.2, and (h) 0.5.
(a)(b)(c)(d)(e)(f)(g)(h)Table 3
PSNR (dB) and SSIM of original images with various intensities of white Gaussian noise.
k = 0.01k = 0.05k = 0.1PSNRImage “Boat”20.0719.0717.03Image “Elaine”20.0519.1917.24Image “Peppers”20.1719.1417.06SSIMImage “Boat”0.81170.80310.7821Image “Elaine”0.81630.80780.7872Image “Peppers”0.81310.80510.7832Table 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).
Image“Boat”“Elaine”“Peppers”k0.010.050.10.010.050.10.010.050.1Proposed algorithm24.85 (4.78)22.45 (3.38)18.81 (1.78)25.77 (5.72)23.06 (3.87)19.22 (1.98)25.03 (4.86)22.54 (3.50)18.87 (1.81)Mean - value filter [34]21.53 (1.48)18.61 (-0.46)17.82 (0.79)23.53 (3.48)22.54 (3.35)18.92 (1.68)23.25 (3.08)21.81 (2.67)18.34 (1.28)Median-value filter [35]22.67 (2.60)20.23 (1.16)18.57 (1.54)24.25 (4.20)22.91 (3.72)19.05 (1.81)24.70 (4.53)21.79 (2.65)18.63 (1.57)Wiener filter [36]22.58 (2.51)20.02 (0.95)18.14 (1.11)24.65 (4.60)22.73 (3.54)19.15 (1.91)24.61 (4.44)21.85 (2.71)18.42 (1.36)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).
Image“Boat”“Elaine”“Peppers”k0.010.050.10.010.050.10.010.050.1Proposed algorithm0.8321 (0.0204)0.8196 (0.0165)0.8003 (0.0181)0.8354 (0.0191)0.8211 (0.0133)0.8101 (0.0229)0.8331 (0.0200)0.8223 (0.0172)0.7991 (0.0159)Mean - value filter [24]0.8142 (0.0025)0.7989 (-0.0042)0.7867 (0.0046)0.8241 (0.0078)0.8193 (0.0115)0.8023 (0.0151)0.8267 (0.0136)0.8171 (0.0120)0.7956 (0.0124)Median-value filter [25]0.8210 (0.0093)0.8132 (0.0054)0.7978 (0.0157)0.8286 (0.0123)0.8203 (0.0125)0.8031 (0.0159)0.8316 (0.0185)0.8210 (0.0159)0.7981 (0.0149)Wiener filter [26]0.8178 (0.0071)0.8112 (0.0034)0.7892 (0.0071)0.8312 (0.0149)0.8196 (0.0118)0.8054 (0.0176)0.8312 (0.0181)0.8201 (0.0150)0.7953 (0.0121)As can be seen in Tables3–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.


## 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(15)MSE=1MN∑x=1M∑y=1Nfx,y−f′x,y2,wherefx,y denotes the original image and f′x,y denotes the reconstructed image.(2)
The structural similarity (SSIM) is defined as(16)SSIMf,f′=2μfμf′+C12σff′+C2μf2+μf′2+C1σf2+σf′2+C2,whereμf is the mean value of f, μf′ is the mean value of f′, σf2 is the variance of f, and σf2 is the variance of f′ corresponding to the special case of C1=C2=0 [32]. The greater the value of SSIM, the better the quality of the reconstructed image.Figure 6
(a) Original “House” image; (b–g) SPs in the directions of +15°, +45°, and +75°, -−15°, −45°, and −75° at scale22; (h) reconstructed image.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 7
PSNR values with different methods: the five points on the abscissa denote the reconstructed images of “House,” “Boat,” “Elaine,” “Peppers,” and “Barbara.”The optimized scale for the DT-CWT analysis and synthesis process is selected from scale2J. 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.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 indexJ = 3.Figure 9
Relationship between PSNR of the reconstructed image and the scale indexj for five iterations per image.Figure 10
Relationship between the peak structural similarity (SSIM) of the reconstructed image and the number of iterations for each image for the scale indexJ = 3.Figure 11
Relationship between SSIM of the reconstructed image and the scale indexj for five iterations per image.


## 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. Figure6 shows the SPs of the “House” image at scale 22 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 L3fx,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. However, these reconstructed images retain the local features of the original image.Figure 12
First row: original images: (a) “Boat,” (b) “Elaine,” (c) “Barbara,” and (d) “Peppers.” Second row: low-frequency images subsampled by a factor of 23 along their rows and columns. Third row: reconstructed images.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Table1 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].Table 1
Structural similarity (SSIM), mean squared error (MSE), and PSNR (dB) values of the reconstructed images.
SchemeImageNo. Of SPsRatio of SPs over original image (%)SSIMMSEPSNRProposed“Boat”197549.100.8383104.9727.92“Elaine”149447.260.859354.2130.79“Peppers”189218.780.864288.5228.66Ref. [12]“Boat”195379.020.8123134.4526.72“Elaine”151227.330.8201117.8627.24“Peppers”190328.820.8165122.3427.12Ref. [33]“Boat”201319.240.8013198.7825.13“Elaine”150427.300.8196163.1326.05“Peppers”191248.850.8125186.9325.27The bold values show better results than the other methods.Figure 13
SSIM values with different methods: the five points on the abscissa denote the reconstructed images of “House,” “Boat,” “Elaine,” “Peppers,” and “Barbara.”As can be seen in Table1, 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.Table 2
Analysis of the IR scheme with scale indexJ=3 and five iterations.
SchemeImageTime (s)PSNR (dB)SSIMProposed algorithm“Boat”5.90127.920.8383“Elaine”5.37930.790.8593“Peppers”5.63428.660.8642Ref. [12]“Boat”7.55626.720.8123“Elaine”7.23227.240.8201“Peppers”7.44127.120.8165Ref. [33]“Boat”6.02425.130.8013“Elaine”5.85726.050.8196“Peppers”5.93225.270.8125Figure 14
Time required for and proportion of each procedure of “Boat” reconstruction. The time required for the five iterations of IR is 5.901 s for the maximum scale indexJ = 3.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:(17)C′=C+kG,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.Figure 15
Comparison of denoising performances with various intensities of white Gaussian noise. First row: original image with added white Gaussian noise,k = (a) 0.01, (b) 0.1, (c) 0.2, and (d) 0.5. Second row: reconstructed image with k = (e) 0.01, (f) 0.1, (g) 0.2, and (h) 0.5.
(a)(b)(c)(d)(e)(f)(g)(h)Table 3
PSNR (dB) and SSIM of original images with various intensities of white Gaussian noise.
k = 0.01k = 0.05k = 0.1PSNRImage “Boat”20.0719.0717.03Image “Elaine”20.0519.1917.24Image “Peppers”20.1719.1417.06SSIMImage “Boat”0.81170.80310.7821Image “Elaine”0.81630.80780.7872Image “Peppers”0.81310.80510.7832Table 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).
Image“Boat”“Elaine”“Peppers”k0.010.050.10.010.050.10.010.050.1Proposed algorithm24.85 (4.78)22.45 (3.38)18.81 (1.78)25.77 (5.72)23.06 (3.87)19.22 (1.98)25.03 (4.86)22.54 (3.50)18.87 (1.81)Mean - value filter [34]21.53 (1.48)18.61 (-0.46)17.82 (0.79)23.53 (3.48)22.54 (3.35)18.92 (1.68)23.25 (3.08)21.81 (2.67)18.34 (1.28)Median-value filter [35]22.67 (2.60)20.23 (1.16)18.57 (1.54)24.25 (4.20)22.91 (3.72)19.05 (1.81)24.70 (4.53)21.79 (2.65)18.63 (1.57)Wiener filter [36]22.58 (2.51)20.02 (0.95)18.14 (1.11)24.65 (4.60)22.73 (3.54)19.15 (1.91)24.61 (4.44)21.85 (2.71)18.42 (1.36)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).
Image“Boat”“Elaine”“Peppers”k0.010.050.10.010.050.10.010.050.1Proposed algorithm0.8321 (0.0204)0.8196 (0.0165)0.8003 (0.0181)0.8354 (0.0191)0.8211 (0.0133)0.8101 (0.0229)0.8331 (0.0200)0.8223 (0.0172)0.7991 (0.0159)Mean - value filter [24]0.8142 (0.0025)0.7989 (-0.0042)0.7867 (0.0046)0.8241 (0.0078)0.8193 (0.0115)0.8023 (0.0151)0.8267 (0.0136)0.8171 (0.0120)0.7956 (0.0124)Median-value filter [25]0.8210 (0.0093)0.8132 (0.0054)0.7978 (0.0157)0.8286 (0.0123)0.8203 (0.0125)0.8031 (0.0159)0.8316 (0.0185)0.8210 (0.0159)0.7981 (0.0149)Wiener filter [26]0.8178 (0.0071)0.8112 (0.0034)0.7892 (0.0071)0.8312 (0.0149)0.8196 (0.0118)0.8054 (0.0176)0.8312 (0.0181)0.8201 (0.0150)0.7953 (0.0121)As can be seen in Tables3–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.



---

*Source: 6752486-2021-12-21.xml*