# Image Denoising Based on Dilated Singularity Prior

Authors: Shaoxiang Hu; Zhiwu Liao; Wufan Chen
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

In order to preserve singularities in denoising, we propose a new scheme
by adding dilated singularity prior to noisy images. The singularities are detected by
canny operator firstly and then dilated using mathematical morphology for finding
pixels “near” singularities instead of “on” singularities. The denoising results for pixels
near singularities are obtained by nonlocal means in spatial domain to preserve singularities
while the denoising results for pixels in smooth regions are obtained by EM
algorithm constrained by a mask formed by downsampled spatial image with dilated
singularity prior to suiting the sizes of the subbands of wavelets. The final denoised results
are got by combining the above two results. Experimental results show that the scheme
can preserve singularity well with relatively high PSNR and good visual quality.


---

## Body

## 1. Introduction

How to preserve edges and textures in image denoising is a very difficult and important topic. Many efforts focus on this topic for a long time. The most straight way to image denoising while preserving singularities is adaptive methods [1, 2], which suggest that denoising should be carried on local regions. However, preserving singularities requires the size of local area to be small while denoising requires the size of local area to be large, which makes it become a trade-off problem.Recently, Buades et al. argue that the small size of local area is not a requirement for preserving singularities in image denoising. That is, if correct similar pixels can be found through their similar local areas, singularity preserving denoising can be achieved by simple mean smoothing [3]. Nonlocal means has very good performance in maintaining the borders and textures than most of existing methods. Based on it, some improved algorithms are proposed [4–10]. However, its computation cost is very high.Although some denoising methods on spatial domain, such as nonlocal means and its improvements, nonlinear anisotropic diffusion, and fractional calculus can provide satisfied denoising results, most of them suffer from high computation burden [9, 11–20]. Thus, some researchers suggest that wavelets can be used to reduce the computation cost and suppress noises.Some works on wavelet domain propose that threshold values should be selected according to statistical rules; the type of these methods is called shrinkage denoising. In [21], the threshold is determined by SURE shrinkage, which has better performance than threshold selected by Baysian rule. However, all coefficients in a subband are considered as independent identical distributed (iid) variables, and this is considered an unreal assumption.In order to improve the performance of shrinkages, some schemes focus on how to describe the residual dependencies among wavelet coefficients [22–26]. In [22], the authors present a tree structure to describe the dependency among the wavelet coefficients. Recently, [23] proposes that trivariates should be used for three subbands of wavelet coefficients.Both methods can provide improved denoising results than independent model for wavelet coefficients, but they do not consider the difference between pixels near singularities and in smooth regions. That is, they think that the coefficients near singularities are the same as they are in smooth regions. Intuitively, it is wrong.We suggest that the dependency among wavelet coefficients should be considered both on their local energy and their object labels [26]. However, object labels cannot be assigned correctly in serious noise especially in high-frequency subbands.Most recently, following block similarity proposed in [3], a denoised method based on collaborative filtering and wavelets is proposed. It obtains great success in its high value of peak signal noise ratio (PSNR). However, estimating real levels both of singularity points and smooth points with large regions leads it has very high computation complexity.In summary, there are two popular methods to preserve singularities in denoising: one is local area similarity proposed in [3] and then improved by many efforts; the other is adaptive denoising in small local regions. The former has good performance in preserving singularities but suffers from its high computation complexity while the latter has low computation complexity and good performance in smooth regions but blurs the singularities.It reminds us that if we can distinguish points near singularities from the points in smooth regions correctly and design different denoising schemes for these two types of points separately, singularities can be preserved in denoising. Following this idea, dilated singularity prior is posed to find points “near” singularities instead of “on” singularities. After distinguishing pixels near singularities from pixels in smooth regions, the real level for pixels in smooth regions can be estimated by a smooth version of noisy image. By combining these two different denoising methods together, we can obtain good denoising results.The detailed framework will be introduced as follows. The second section in this paper will introduce backgrounds, the third section introduces how to build mask, and the fourth section describes the denoising framework. Section5 is the experimental results and discussion. Section 6 gives conclusions, and, finally, there is the acknowledgment part.


## 2. Backgrounds

In this section, we will introduce terms, notations, and concepts for image denoising, wavelets, and GMM.


### 2.1. Image Denoising

An imagex is regarded as a realization of a random field X. In these terms, the image denoising problem can be rephrased as follows: given a noisy image, estimate for real gray level of each pixel.Let the original image be{xi,j},i,j=1,…,N, where N is some integer power of 2 and (i,j) is the pixel coordinate of the image. The image has been corrupted by additive noise and one observation(2.1)yi,j=xi,j+ni,j,
where ni,j is an iid variable as normal N(0,σn2) and independent to xi,j. The goal is to remove the noise from yi,j and obtain an estimate x̂i,j of xi,j which minimizes the mean squared error (MSE):(2.2)MSE(x̂)=1N2∑i,j=1N[x̂i,j-xi,j]2.
Thus, the minimum MSE estimate would be the conditional mean estimate of a Gaussian signal with Gausssian noise(2.3)E[Xi,j∣Yi,j]=σXi,j2σXi,j2+σn2yi,j,
where Xi,j and Yi,j are the random variables of xi,j and yi,j respectively.However, only one observation (a noisy image) provides for us to estimate real gray levels for the image. Thus we have to share statistical information among image pixels. One common method used in adaptive denoising assuming pixels in a local region is iid variables. However, since local regions have singularities, iid variable assumption leads to blurred edges and textures.Intuitively, pixels with more similar real gray levels have higher probability for identical distribution. Thus more plausible way for selecting iid pixels converts to find pixels with similar real gray levels in a local region.However, finding pixels with similar real gray levels in noise is very difficult. In this paper, we use singularities detected in spatial domain combined with a mask built on LL subband of wavelets to help us find similar pixels in a local region. The whole framework will be given in Section4.


### 2.2. Wavelets and Gaussian Mixture Model (GMM) [21]

LetW=wx, W denote the matrix of wavelet coefficients of x, where w is the 2D dyadic orthogonal wavelet transform operator. It is convenient to group wavelet coefficients into groups or subbands of different scales and orientations, where, for example, the label HL1 refers to those coefficients at the first scale of decomposition which are the output of the high-pass filter in the horizontal direction and the low-pass filter in the vertical direction. The subbands HHk,HLk,LHk, k=1,2,…,J are called the details, where k is the scale, with J being the largest (or coarsest) scale in the decomposition, and a subband at scale k has size N/2k×N/2k. The subband LLJ is the low-resolution residual. Wavelets lead to a natural structure of the wavelet coefficients into three subbands representing the horizontal, vertical, and diagonal edges.The coefficients of DWTs of many real-world images tend to be sparse, with just a few nonzero coefficients containing most of the image energy. This property, combined with our view of the image as a realization from a family or distribution of signal, leads to a simple model for an individual wavelet coefficient. Hence, the marginal density of each wavelet coefficient is typically described by a peaky and heavy-tailed non-Gaussian density.Such densities are well approximated by a two-component Gaussian mixture model [21], and its overall pdf is given by(2.4)fWi,j(ωi,j)=∑m=01PSi,jYWi,j∣Si,j(ωi,j∣Si,j=m).In general, a 2-state Gaussian mixture model for a random variableWi,j consists of(1)
discrete random state variableSi,j taking the value Si,j∈{0,1} according to the probability mass function (pmf);(2)
means and variances of Gaussian distributions.In most applications of mixture models, the valueωi,j, can be observed, but the value of the state variable is not; we say that the value of Si,j is hidden. Although each wavelet coefficient Wi,j is conditionally Gaussian given its state variable Si,j, the wavelet coefficient has an overall non-Gaussian density due to the randomness of Si,j.


### 2.3. EM Algorithm

EM algorithm is a kind of Max-likelihood Estimation (MLE) algorithm whose target is to find a set of hidden states to maximize the probability of observations. If the initial values for parameters are known, EM iterates between estimating the probability for the state (Expectation) and updating the model given the state probabilities (Maximization).The normal form of the E step and M step of EM algorithm is as follows.Expectation (E):(2.5)PSi,j∣Wi,j(m∣ωi,j)=PSi,j(m)PWi,j|Si,j(ωi,j∣m)∑m=01PSi,j(m)PWi,j∣Si,j(ωi,j∣m).Maximization (M):(2.6)PSi,j(m)=1N∑t,kPSt,k|Wt,k(m∣ωt,k),σi,j,m2=1NPSi,j(m)∑t,kωt,k2PSt,k|Wt,k(m∣ωt,k),
where N is the samples of an image.


## 2.1. Image Denoising

An imagex is regarded as a realization of a random field X. In these terms, the image denoising problem can be rephrased as follows: given a noisy image, estimate for real gray level of each pixel.Let the original image be{xi,j},i,j=1,…,N, where N is some integer power of 2 and (i,j) is the pixel coordinate of the image. The image has been corrupted by additive noise and one observation(2.1)yi,j=xi,j+ni,j,
where ni,j is an iid variable as normal N(0,σn2) and independent to xi,j. The goal is to remove the noise from yi,j and obtain an estimate x̂i,j of xi,j which minimizes the mean squared error (MSE):(2.2)MSE(x̂)=1N2∑i,j=1N[x̂i,j-xi,j]2.
Thus, the minimum MSE estimate would be the conditional mean estimate of a Gaussian signal with Gausssian noise(2.3)E[Xi,j∣Yi,j]=σXi,j2σXi,j2+σn2yi,j,
where Xi,j and Yi,j are the random variables of xi,j and yi,j respectively.However, only one observation (a noisy image) provides for us to estimate real gray levels for the image. Thus we have to share statistical information among image pixels. One common method used in adaptive denoising assuming pixels in a local region is iid variables. However, since local regions have singularities, iid variable assumption leads to blurred edges and textures.Intuitively, pixels with more similar real gray levels have higher probability for identical distribution. Thus more plausible way for selecting iid pixels converts to find pixels with similar real gray levels in a local region.However, finding pixels with similar real gray levels in noise is very difficult. In this paper, we use singularities detected in spatial domain combined with a mask built on LL subband of wavelets to help us find similar pixels in a local region. The whole framework will be given in Section4.


## 2.2. Wavelets and Gaussian Mixture Model (GMM) [21]

LetW=wx, W denote the matrix of wavelet coefficients of x, where w is the 2D dyadic orthogonal wavelet transform operator. It is convenient to group wavelet coefficients into groups or subbands of different scales and orientations, where, for example, the label HL1 refers to those coefficients at the first scale of decomposition which are the output of the high-pass filter in the horizontal direction and the low-pass filter in the vertical direction. The subbands HHk,HLk,LHk, k=1,2,…,J are called the details, where k is the scale, with J being the largest (or coarsest) scale in the decomposition, and a subband at scale k has size N/2k×N/2k. The subband LLJ is the low-resolution residual. Wavelets lead to a natural structure of the wavelet coefficients into three subbands representing the horizontal, vertical, and diagonal edges.The coefficients of DWTs of many real-world images tend to be sparse, with just a few nonzero coefficients containing most of the image energy. This property, combined with our view of the image as a realization from a family or distribution of signal, leads to a simple model for an individual wavelet coefficient. Hence, the marginal density of each wavelet coefficient is typically described by a peaky and heavy-tailed non-Gaussian density.Such densities are well approximated by a two-component Gaussian mixture model [21], and its overall pdf is given by(2.4)fWi,j(ωi,j)=∑m=01PSi,jYWi,j∣Si,j(ωi,j∣Si,j=m).In general, a 2-state Gaussian mixture model for a random variableWi,j consists of(1)
discrete random state variableSi,j taking the value Si,j∈{0,1} according to the probability mass function (pmf);(2)
means and variances of Gaussian distributions.In most applications of mixture models, the valueωi,j, can be observed, but the value of the state variable is not; we say that the value of Si,j is hidden. Although each wavelet coefficient Wi,j is conditionally Gaussian given its state variable Si,j, the wavelet coefficient has an overall non-Gaussian density due to the randomness of Si,j.


## 2.3. EM Algorithm

EM algorithm is a kind of Max-likelihood Estimation (MLE) algorithm whose target is to find a set of hidden states to maximize the probability of observations. If the initial values for parameters are known, EM iterates between estimating the probability for the state (Expectation) and updating the model given the state probabilities (Maximization).The normal form of the E step and M step of EM algorithm is as follows.Expectation (E):(2.5)PSi,j∣Wi,j(m∣ωi,j)=PSi,j(m)PWi,j|Si,j(ωi,j∣m)∑m=01PSi,j(m)PWi,j∣Si,j(ωi,j∣m).Maximization (M):(2.6)PSi,j(m)=1N∑t,kPSt,k|Wt,k(m∣ωt,k),σi,j,m2=1NPSi,j(m)∑t,kωt,k2PSt,k|Wt,k(m∣ωt,k),
where N is the samples of an image.


## 3. Designing Mask

Just as discussed in Section2.1, image denoising using statistical theory has to find similar pixels as iid variables. However, pixels near singularities have very different statistical properties to pixels in smooth regions. Thus, if two types of pixels, pixels near singularities and pixels in smooth regions, can be parted correctly, iid variables can be selected from the pixels with identical type labels. Based on this idea, we present a method to find similar pixels by posing dilated singularity prior. That is, singularity prior and its dilated version are used to find pixels near singularities while other pixels are considered as in smooth regions.However, polluted gray levels will also lead to noisy wavelet coefficients in the details of wavelets, which hampers us to detect singularities and measure the similarity among wavelet coefficients correctly. Therefore, posing singularity prior and finding similar pixels in the details of wavelets cannot be carried on directly using the details of wavelets.One solution for these difficulties is to design a mask and then apply it to the details of wavelets for image denoising. In this section, we will discuss the method to design the mask in detail.


### 3.1. Dilated Singularity Prior

Pixels near singularities are different to the pixels in smooth regions. Thus if we can find the pixels near singularities correctly and handle them separately, the singularities can be preserved in denoising. On the other side, smoothing only among pixels in smooth regions, which is implemented by abandoning pixels near singularities, can help us reduce the outliers in estimate. Moreover, parting pixels to different groups can help us design different denoising methods separately according to their different statistical properties. Thus it can get satisfying denoising results for both groups.However, in noise, singularities in the details of wavelets cannot be detected correctly. Moreover, the low-resolution residual of wavelets is an oversmoothing version for the noisy image, which blurs many important weak singularities, such as some textures and degraded edges. Therefore, singularities cannot be detected from the subbands of the wavelets directly. In this paper, the singularities are located in spatial domain using canny operator and then downsampled to suit the size of different subbands of wavelets.It should be indicated that although some pixels near singularities have similar real gray levels to nearby pixels in smooth regions, they have different statistical natures to the pixels in nearby smooth region. Thus, we need find pixelsnear singularities instead of on singularities. Based on this consideration, the detected singularities by canny are dilated using mathematical morph operator.In summary, procedures for posing singularity prior are(1)
detecting singularities in the noisy image using canny operator;(2)
dilating singularities using mathematical morphology operator;(3)
downsampling singularity image obtained in procedure 2 to form singularity prior for different subbands of wavelets.One example for these three steps of noisy512×512 barbara whose noise variance is 0.01 is shown in Figures 1(a)–1(c).The steps for locating singularity prior and making the mask. (a)–(c) three steps for locating singularity prior in Section3.1; (d) the mask in Section 3.2.
(a)
Step 1: the singularities detected by canny in noisy barbara(b)
Step 2: the dilated singularities(c)
Step 3: the downsampled dilated singularity prior(d)
The mask composed byLL1 and singularity prior


### 3.2. Designing Mask

After parting pixels into two groups by posing dilated singularity prior, we have to measure the similarity for two groups correctly in noise. However, even having posed singularity prior, measuring similarity among pixels in smooth regions correctly in the details of wavelets is very difficult because of noise influence in high frequency.Fortunately, the low-resolution residual of wavelets can provide a downsampled smoother version for noisy images. Thus, if the singularity prior is posed to the low-resolution residual of wavelets, the similarity among pixels in smooth regions can be measured correctly by the coefficients of the low-resolution residual. Based on this discussion, the mask for the details of wavelets is designed as posing singularity prior to the low resolution of wavelets.Figures1(a)–1(c) show the steps for building the singularity prior for pixels in smooth regions, and the final mask is shown in Figure 1(d).


## 3.1. Dilated Singularity Prior

Pixels near singularities are different to the pixels in smooth regions. Thus if we can find the pixels near singularities correctly and handle them separately, the singularities can be preserved in denoising. On the other side, smoothing only among pixels in smooth regions, which is implemented by abandoning pixels near singularities, can help us reduce the outliers in estimate. Moreover, parting pixels to different groups can help us design different denoising methods separately according to their different statistical properties. Thus it can get satisfying denoising results for both groups.However, in noise, singularities in the details of wavelets cannot be detected correctly. Moreover, the low-resolution residual of wavelets is an oversmoothing version for the noisy image, which blurs many important weak singularities, such as some textures and degraded edges. Therefore, singularities cannot be detected from the subbands of the wavelets directly. In this paper, the singularities are located in spatial domain using canny operator and then downsampled to suit the size of different subbands of wavelets.It should be indicated that although some pixels near singularities have similar real gray levels to nearby pixels in smooth regions, they have different statistical natures to the pixels in nearby smooth region. Thus, we need find pixelsnear singularities instead of on singularities. Based on this consideration, the detected singularities by canny are dilated using mathematical morph operator.In summary, procedures for posing singularity prior are(1)
detecting singularities in the noisy image using canny operator;(2)
dilating singularities using mathematical morphology operator;(3)
downsampling singularity image obtained in procedure 2 to form singularity prior for different subbands of wavelets.One example for these three steps of noisy512×512 barbara whose noise variance is 0.01 is shown in Figures 1(a)–1(c).The steps for locating singularity prior and making the mask. (a)–(c) three steps for locating singularity prior in Section3.1; (d) the mask in Section 3.2.
(a)
Step 1: the singularities detected by canny in noisy barbara(b)
Step 2: the dilated singularities(c)
Step 3: the downsampled dilated singularity prior(d)
The mask composed byLL1 and singularity prior


## 3.2. Designing Mask

After parting pixels into two groups by posing dilated singularity prior, we have to measure the similarity for two groups correctly in noise. However, even having posed singularity prior, measuring similarity among pixels in smooth regions correctly in the details of wavelets is very difficult because of noise influence in high frequency.Fortunately, the low-resolution residual of wavelets can provide a downsampled smoother version for noisy images. Thus, if the singularity prior is posed to the low-resolution residual of wavelets, the similarity among pixels in smooth regions can be measured correctly by the coefficients of the low-resolution residual. Based on this discussion, the mask for the details of wavelets is designed as posing singularity prior to the low resolution of wavelets.Figures1(a)–1(c) show the steps for building the singularity prior for pixels in smooth regions, and the final mask is shown in Figure 1(d).


## 4. The Framework

In this section, we will discuss the framework for the proposed method, whose flow chart is shown in Figure2. In Figure 2, the squares without shadow represent the method used in this framework while the shaded squares represent the processing results for the methods, and ⊕ represents combining two squares together. The dashed frame with a title frame “Mask” is the steps for designing masks for different pixel groups presented in Section 3.Figure 2
The flowchart of proposed method.The proposed method is based on a novel idea of preserving singularities by processing pixels near singularities and in smooth regions separately. The steps for denoising to different types of pixels are shown in Figure2 as the two dashed frames with number 1 and 2 in two triangles where number 1 indicates the denoising steps for pixels in smooth regions and number 2 indicates denoising steps for pixels near singularities. Both methods will be discussed detailedly in this section.


### 4.1. Denoising for Pixels Near Singularities

In order to preserve singularities, the pixels near singularities located by dilated singularity prior should be processed alone. Since LL subband of wavelets will blur the singularities while the details of wavelets will magnify the influence of noises, the denoising for pixels near singularities should be carried on spatial domain directly. In this subsection, we will present the denoising method for pixels near singularities.For each pixelyi,j of image, after imposing dilated singularity prior, it has a one and only one label fi,j and fi,j∈{0,1}, where 0 represents “non-singularity” and 1 represents “near singularity.”Essentially, the most important procedure for image denoising is how to select similar points from image pixels. Similar pixels selected according to their block similarity are proposed in [3] recently. In this framework, the block is defined as a square with fixed size centered at the consideration pixel. In our method, the similarity among two pixels yi,j, ys,t with labels fi,j=1 and fs,t=1 can be measured by two 7×7 blocks centered at (i,j) and (s,t):(4.1)S(yi,j,ys,t)=∑k=-3,l=-3k=3,l=3‖yi-k,j-l-ys-k,t-l‖2e-(k2+l2)/2σ2.Thus, for each pixelyi,j with fi,j=1, searches in a 21×21 window to find its similar pixels whose similarity between itself and yi,j is below a predefined threshold T:(4.2)χ(s,t)={1:S(yi,j,ys,t)≤T,0:S(yi,j,ys,t)>T.Thus the estimate of the real gray level ofyi,j is(4.3)ŷi,j=1n∑k,l=-10k,l=10yi-k,j-lχ(i-k,j-l)fi-k,j-l,
where(4.4)n=∑k,l=-10k,l=10χ(i-k,j-l)fi-k,j-l.Sinceχ(i-k,j-l) and fi-k,j-l are two indicator functions, we can combine them as(4.5)χ(i-k,j-l)={1:S(yi,j,yi-k,j-l)≤T,fi-k,j-l=1,0:otherwise.Therefore, (4.3) becomes(4.6)ŷi,j=1n∑k,l=-10k,l=10yi-k,j-lχ(i-k,j-l),
where(4.7)n=∑k,l=-10k,l=10χ(i-k,j-l).According to the above discussion, the denoising steps for pixels near singularities are as follows.(1)
Initialization: give T, compute fi,j for all image pixels according to dilated singularity prior.(2)
For an image pixel(i,j)withfi,j=1
(a)
fork=-10:10, l=-10:10computeS(yi,j,yi-k,j-l) using (4.1), for all fi-k,j-l=1,(b)
ifS≤T, χ(i-k,j-l)=1otherwise,χ(i-k,j-l)=0,(c)
computeŷi,j using (4.6).(3)
Repeat step 2(a)–(c) until all image pixels have been processed.


### 4.2. Denoising for Pixels in Smooth Regions

Wavelets, which are based on the idea that a linear, invertible transform will represent the image by the sparse wavelet coefficients whose structure is “simpler" to process, are a powerful tool to reduce the complexity.However, in noise, there are two difficulties: one is the singularities cannot be located correctly in all subbands; the other is how to find similar pixels in the details of wavelets. In order to solve the above two problems, a mask will be built in the LL subbands and spatial domain together whose detailed algorithm is discussed in Section3. That is, locate prior in spatial domain and then downsample and pose it to the LL subband to form the mask. Therefore, the real values of the details of wavelets can be estimated based on the mask.LetWi,j,kt denote a coefficient of wavelet, where t∈{0,1,2,3} and 0, 1, 2, 3 represent the LL, HL, LH, and HH subband, respectively, and k represents the decomposition scale.Since the mask is built by posing singularity prior to the LL subband of each scale, we can give a labelιi,j,k for a pixel Wi,j,k0 in LL subband and ιi,j,k∈{0,1}, where 0 represents “non-singularity” and 1 represents “near singularity”. Thus ιi,j,ks combined with Wi,j,k0s forms the mask to measure the similarity between each pair of wavelet coefficients.In our method,Wi,j,kt, t=0,1,2,3, is considered as a random variable with GMM distribution(4.8)PWi,j,kt(ωi,j,kt)=∑m=01PSi,j,ktYWi,j,kt∣Si,j,kt(ωi,j,kt∣Si,j,kt=m),
where i,j=1,…,N/2k, t=0,1,2,3.Thus if the parametersΘ={PSi,j,kt=m,σm,i,j,kt,m=0,1} where σm,i,j,kt,m=0,1 are two standard deviations of Gaussian components, have known. The real value of Wi,j,kt can be estimated as(4.9)ω̂i,j,kt=∑m=01P(Si,j,kt=m∣ω,Θ)[(σm,i,j,kt)2-σn2]+(σm,i,j,kt)2,
where(4.10)[(σm,i,j,kt)2-σn2]+={[(σm,i,j,kt)2-σn2]:[(σm,i,j,kt)2-σn2]≥0,0:otherwise.However, only one observation is forW; thus we have to share statistical information to estimate parameter Θ. Generally, the statistical information should be shared among iid variables. Pixels with small spatial distance and similar real values are iid variables which is a plausible assumption. Then we will discuss how to find pixels with small distance and similar real values in the details of wavelets. Since the following steps are processed in a subband and one scale, we will omit index t for t=1,2,3, k for above notations for the sake of simple explanation.The similar pixels ofWi,j can be selected from a (2ν+1)×(2ν+1) window with measure of similarity(4.11)γ(ωi,j,ωs,t)=‖ωi,j0-ωs,t0‖2×ϕ(ωi,j0,ωs,t0),
where ωi,j0 and ωs,t0 are values of (i,j) and (s,t) in the LL subband, respectively; the ϕ(ωi,j0,ωs,t0) is a binary value function(4.12)ϕ(ωi,j0,ωs,t0)={1:ιi,j=ιs,t,0:otherwise.Then by using similarity measure defined in (4.11) we can find similar points in a (2ν+1)×(2ν+1) centered at (i,j)(4.13)Φi,j={(s,t)∣γ(ωi,j,ωs,t)≠0,γ(ωi,j,ωs,t)≤τ,|s-i|≤ν,|t-j|≤ν}.Thus the parameterΘ can be estimated using EM algorithm in the collection Φi,j.Expectation (E):(4.14)PSi,j∣Wi,j(m∣ωi,j)=PSi,j(m)PWi,j∣Si,j(ωi,j∣m)∑m=01PSi,j(m)PWi,j|Si,j(ωi,j∣m).Maximization (M):(4.15)PSi,j(m)=1|Φi,j|∑s,t,(s,t)∈ϕPSs,t|Ws,t(m∣ωs,t),σi,j,m2=1|Φi,j|PSi,j(m)∑s,t,(s,t)∈ϕωs,t2PSs,t∣Ws,t(m∣ωs,t),
where |Φi,j| is the number of elements in Φi,j.In summary, the steps of denoising for pixels in smooth regions are(1)
Initialization: input scale J, dilated singularity prior F whose size is N×N, and ν.(2)
Build marks for the details of different scale using downsample (see Section3).(3)
Find similar pixels for eachWi,j,kt: find similar points using (4.11) and then compute ϕ.(4)
EM algorithm, for each Wi,j,kt, estimate parameter Θ using EM algorithm in (4.14)-(4.15).(5)
Repeat Steps(3)-(4)until all wavelet coefficients have been processed.(6)
Inverse wavelet transform: obtain denoised results using inverse wavelet transform for all processed coefficients.(7)
Get final denoised image: the final denoised image can be got by combining both denoising image together; that is, the real values for pixels near singularities are estimated using the method in Section 4.1, and the values for other pixels are estimated using steps (1)–(6).


## 4.1. Denoising for Pixels Near Singularities

In order to preserve singularities, the pixels near singularities located by dilated singularity prior should be processed alone. Since LL subband of wavelets will blur the singularities while the details of wavelets will magnify the influence of noises, the denoising for pixels near singularities should be carried on spatial domain directly. In this subsection, we will present the denoising method for pixels near singularities.For each pixelyi,j of image, after imposing dilated singularity prior, it has a one and only one label fi,j and fi,j∈{0,1}, where 0 represents “non-singularity” and 1 represents “near singularity.”Essentially, the most important procedure for image denoising is how to select similar points from image pixels. Similar pixels selected according to their block similarity are proposed in [3] recently. In this framework, the block is defined as a square with fixed size centered at the consideration pixel. In our method, the similarity among two pixels yi,j, ys,t with labels fi,j=1 and fs,t=1 can be measured by two 7×7 blocks centered at (i,j) and (s,t):(4.1)S(yi,j,ys,t)=∑k=-3,l=-3k=3,l=3‖yi-k,j-l-ys-k,t-l‖2e-(k2+l2)/2σ2.Thus, for each pixelyi,j with fi,j=1, searches in a 21×21 window to find its similar pixels whose similarity between itself and yi,j is below a predefined threshold T:(4.2)χ(s,t)={1:S(yi,j,ys,t)≤T,0:S(yi,j,ys,t)>T.Thus the estimate of the real gray level ofyi,j is(4.3)ŷi,j=1n∑k,l=-10k,l=10yi-k,j-lχ(i-k,j-l)fi-k,j-l,
where(4.4)n=∑k,l=-10k,l=10χ(i-k,j-l)fi-k,j-l.Sinceχ(i-k,j-l) and fi-k,j-l are two indicator functions, we can combine them as(4.5)χ(i-k,j-l)={1:S(yi,j,yi-k,j-l)≤T,fi-k,j-l=1,0:otherwise.Therefore, (4.3) becomes(4.6)ŷi,j=1n∑k,l=-10k,l=10yi-k,j-lχ(i-k,j-l),
where(4.7)n=∑k,l=-10k,l=10χ(i-k,j-l).According to the above discussion, the denoising steps for pixels near singularities are as follows.(1)
Initialization: give T, compute fi,j for all image pixels according to dilated singularity prior.(2)
For an image pixel(i,j)withfi,j=1
(a)
fork=-10:10, l=-10:10computeS(yi,j,yi-k,j-l) using (4.1), for all fi-k,j-l=1,(b)
ifS≤T, χ(i-k,j-l)=1otherwise,χ(i-k,j-l)=0,(c)
computeŷi,j using (4.6).(3)
Repeat step 2(a)–(c) until all image pixels have been processed.


## 4.2. Denoising for Pixels in Smooth Regions

Wavelets, which are based on the idea that a linear, invertible transform will represent the image by the sparse wavelet coefficients whose structure is “simpler" to process, are a powerful tool to reduce the complexity.However, in noise, there are two difficulties: one is the singularities cannot be located correctly in all subbands; the other is how to find similar pixels in the details of wavelets. In order to solve the above two problems, a mask will be built in the LL subbands and spatial domain together whose detailed algorithm is discussed in Section3. That is, locate prior in spatial domain and then downsample and pose it to the LL subband to form the mask. Therefore, the real values of the details of wavelets can be estimated based on the mask.LetWi,j,kt denote a coefficient of wavelet, where t∈{0,1,2,3} and 0, 1, 2, 3 represent the LL, HL, LH, and HH subband, respectively, and k represents the decomposition scale.Since the mask is built by posing singularity prior to the LL subband of each scale, we can give a labelιi,j,k for a pixel Wi,j,k0 in LL subband and ιi,j,k∈{0,1}, where 0 represents “non-singularity” and 1 represents “near singularity”. Thus ιi,j,ks combined with Wi,j,k0s forms the mask to measure the similarity between each pair of wavelet coefficients.In our method,Wi,j,kt, t=0,1,2,3, is considered as a random variable with GMM distribution(4.8)PWi,j,kt(ωi,j,kt)=∑m=01PSi,j,ktYWi,j,kt∣Si,j,kt(ωi,j,kt∣Si,j,kt=m),
where i,j=1,…,N/2k, t=0,1,2,3.Thus if the parametersΘ={PSi,j,kt=m,σm,i,j,kt,m=0,1} where σm,i,j,kt,m=0,1 are two standard deviations of Gaussian components, have known. The real value of Wi,j,kt can be estimated as(4.9)ω̂i,j,kt=∑m=01P(Si,j,kt=m∣ω,Θ)[(σm,i,j,kt)2-σn2]+(σm,i,j,kt)2,
where(4.10)[(σm,i,j,kt)2-σn2]+={[(σm,i,j,kt)2-σn2]:[(σm,i,j,kt)2-σn2]≥0,0:otherwise.However, only one observation is forW; thus we have to share statistical information to estimate parameter Θ. Generally, the statistical information should be shared among iid variables. Pixels with small spatial distance and similar real values are iid variables which is a plausible assumption. Then we will discuss how to find pixels with small distance and similar real values in the details of wavelets. Since the following steps are processed in a subband and one scale, we will omit index t for t=1,2,3, k for above notations for the sake of simple explanation.The similar pixels ofWi,j can be selected from a (2ν+1)×(2ν+1) window with measure of similarity(4.11)γ(ωi,j,ωs,t)=‖ωi,j0-ωs,t0‖2×ϕ(ωi,j0,ωs,t0),
where ωi,j0 and ωs,t0 are values of (i,j) and (s,t) in the LL subband, respectively; the ϕ(ωi,j0,ωs,t0) is a binary value function(4.12)ϕ(ωi,j0,ωs,t0)={1:ιi,j=ιs,t,0:otherwise.Then by using similarity measure defined in (4.11) we can find similar points in a (2ν+1)×(2ν+1) centered at (i,j)(4.13)Φi,j={(s,t)∣γ(ωi,j,ωs,t)≠0,γ(ωi,j,ωs,t)≤τ,|s-i|≤ν,|t-j|≤ν}.Thus the parameterΘ can be estimated using EM algorithm in the collection Φi,j.Expectation (E):(4.14)PSi,j∣Wi,j(m∣ωi,j)=PSi,j(m)PWi,j∣Si,j(ωi,j∣m)∑m=01PSi,j(m)PWi,j|Si,j(ωi,j∣m).Maximization (M):(4.15)PSi,j(m)=1|Φi,j|∑s,t,(s,t)∈ϕPSs,t|Ws,t(m∣ωs,t),σi,j,m2=1|Φi,j|PSi,j(m)∑s,t,(s,t)∈ϕωs,t2PSs,t∣Ws,t(m∣ωs,t),
where |Φi,j| is the number of elements in Φi,j.In summary, the steps of denoising for pixels in smooth regions are(1)
Initialization: input scale J, dilated singularity prior F whose size is N×N, and ν.(2)
Build marks for the details of different scale using downsample (see Section3).(3)
Find similar pixels for eachWi,j,kt: find similar points using (4.11) and then compute ϕ.(4)
EM algorithm, for each Wi,j,kt, estimate parameter Θ using EM algorithm in (4.14)-(4.15).(5)
Repeat Steps(3)-(4)until all wavelet coefficients have been processed.(6)
Inverse wavelet transform: obtain denoised results using inverse wavelet transform for all processed coefficients.(7)
Get final denoised image: the final denoised image can be got by combining both denoising image together; that is, the real values for pixels near singularities are estimated using the method in Section 4.1, and the values for other pixels are estimated using steps (1)–(6).


## 5. Experimental Results

In order to compare our method with state-of-the-art methods and denoising methods using wavelets, the proposed method is compared with nonlocal means [3], BM3D [10], BLS-GSM [21], and HMT [22]. Nonlocal means [3] and BM3D [10] are two well-known state-of-the-art methods. That is, nonlocal means and BM3D find similar points using the block similarity between a pair of pixels, which is proposed firstly in [3]. BM3D designs a complex hierarchical scheme to get very high PSNR, which is usually used as benchmark in many denoising scheme.Moreover, HMT which is presented in [22] is a famous wavelet denoising method capturing dependency among the details of wavelets using a tree structure while BLS-GSM is a wavelet denoising method using adaptive shrinkage and GMM.Since the objective for presenting the denoised results of our method only focuses on whether adding singularity prior can improve the denoising performance, we use some simple tools for testing the potential of our new framework. The wavelet used in this paper is one scale “Haar”; similarity thresholdτ, ν in Φi,j is from 14 to 28 and 3, respectively. Even by using this simple form, the proposed method has good performance both in PSNR and visual quality (see Figure 3). That is, its average PSNR only is 0.2 db lower than the nonlocal means [3]; the average PSNR for lena is 0.4 db lower than the BM3D, and for Barbara, it is 1.4 db lower than BM3D (see Table 1).Table 1
PSNR Comparison of HMT [22], BLS-GSM [21], nonlocal [3], BM3D [10], and the proposed method with simulated gaussian noisy image.
ImageStandard deviation of noiseHMT[22]BLS-GSM[21]Nonlocal[3]BM3D[10]Proposed  methodLena(512×512)1034.0835.2335.6535.9135.511531.8933.5034.1234.2633.852030.4132.2632.8233.0532.662529.3631.3031.8932.0931.71Barbara(512×512)1031.6233.1334.0634.9333.741529.2830.7631.7833.0531.672027.8129.1030.4131.7230.162526.6127.8029.2930.6429.00The noisy image is shown in (a) whose standard deviation of noise is 25, and the denoised results using different methods are shown in (b)–(f).
(a)(b)(c)(d)(e)(f)Although the values of PSNR for proposed method are lower than both nonlocal means and BM3D, considering both of them being state-of-art denoising methods and their computation complexity, PSNR of proposed method is an acceptable result. Moreover, the main objective for this paper only focuses on designing a new scheme to improve performance in singularity preserving by adding singularity prior in image denoising. We think the experimental results can demonstrate the potential of proposed method. Based on this framework, some different methods can be used to improve the performance further. For example, the performance can be improved by processing pixels near singularities using BM3D, using different wavelets, designing different denoising methods for pixels in smooth regions, or by designing more suitable singularity location tools, and so forth. However, discussion about how to improve our framework further is beyond the scope of this paper.


## 6. Conclusions

In this paper, we propose a new scheme to preserve singularities in image denoising by adding dilated singularities to noisy image. This scheme even that uses simple tools can get satisfying denoising results both in PSNR and visual quality. Since the scheme is based on the fact that pixels near singularities have different statistical properties to the pixels in the smooth regions, it is designed by parting the pixels into two groups using dilated singularity prior and handled pixels with different type labels separately. Dilated singularity prior used to locate the pixels “near” singularities instead of “on” singularities can help us obtain satisfied denoising results and reduce computation complexity. Moreover, denoising in smooth regions can be carried on a smoother version of the noisy image, called mask, since the pixels near singularities have been excluded. Based on the above well-designed framework, the proposed method can obtain good denoising results.



---

*Source: 767613-2012-05-08.xml*