A GEOMETRY GUIDED IMAGE DENOISING SCHEME

During image denoising, it is often difficult to balance between the removal of noise and the preservation of contrast and fine features, especially when the noise is excessive. We propose to efficiently balance the two using segmentation and more general geometry extraction transforms. Explained in Nonlocal-means (NL-means) framework, we introduce mutual position function to ensure averaging is only taken over pixels in the same segmentation phase, and provide a convolution kernel and weight function selection scheme to further improve the performance. To address unreliable segmentation due to more excessive noise, we use feature extraction transform that is more general than segmentation and less sensitive to noise. Unlike most denoising approaches that only work for one type of noise and/or involve heuristic parameter tuning, the proposed method comes with an automatic parameter selection scheme, and can be easily adapted for various types of noise, ranging from Gaussian, Poisson, Rician to Ultrasound noise etc. Comparison with the original NL-means as well as ROF, BM3D, K-SVD on various simulated data, MRI and SEM images indicates potentials of our method.

1. Introduction. Image denoising is a crucial process to remove noise from images for direct visualization or further image analysis. The biggest challenge in image denoising is the balance between the removal of noise and the preservation of salient features, like edges and fine structures. Parameter selection is also an issue, especially when the noise is inhomogeneous. There have been some works aiming to remove noise while preserving edges. These include PDE/variational methods [1,2,3,4], Bayesian approach [5], adaptive smoothing [6], bilateral filtering [7,8,9], hybrid methods [10], NL-means [11] and its improvements and variants [12,13,14,15,16,17,18,19,22,23]. However, most existing methods assume that the noise is of a specific type and/or homogeneous (with fixed parameter(s)). Since real images usually contain various types of noise, the application scope of these methods is narrow. They fail to preserve features well when the noise is excessive, thereby leading to inefficient recovery for heavily contaminated images. In this paper, we propose a fully automatic scheme and its variants to remove inhomogeneous noise of various types while preserving image geometry. This scheme outperforms most of the existing methods, especially when there is excessive noise.
The proposed method is explained using NL-means framework. Instead of taking the local structures into account, NL-means method explores the similarity and repetitive information in the whole image, i.e., non-locally. Suppose the noise-free image u(x) is defined in the bounded domain Ω ⊂ R 2 . Let the observed noisy image be v = u+n with n an additive white Gaussian noise. Then the NL-means algorithm states that the estimated intensity of pixel x is a Gaussian weighted average of the intensities of all the pixels in the image: ∫ Ω w(x, y)v(y)dy (1) where w(x, y) = exp and C(x) = ∫ Ω w(x, y)dy.
Here G α (t) is a Gaussian convolution kernel with standard deviation α, and h acts as a filtering parameter. If we denote D(z) = e −z 2 (a kernel function used to define weights w), then the weight w(x, y) is computed as D In practice, the integration is substituted by summation, and to save the computation cost, two local windows are employed in the implementation. One is the similarity window N (0), a square neighborhood of the origin, that substitutes R 2 in the definition of weight w. The other one is search window S(x) that substitutes the whole image domain Ω in the weighted average equation (1). Since enlarging the search domain beyond a certain size will not guarantee a dramatic improvement of performance but at the cost of computation, the search window is fixed as 15×15 in practice [11]. To remove excessive noise, a high h value is usually used in NL-means, but it ends up with loss of fine features especially edges, leading to blurry boundaries or loss of contrast. To solve these problems, instead of averaging over the whole search window as in the original NL-means, we only average over the same phase. More specifically, as illustrated on the right of Figure 1, we call the blue object in the middle one phase (corresponding to level set function ϕ < 0 [24]) and the white background another phase (level set function ϕ > 0). Suppose A is a pixel in the background with excessive noise (corresponding to phase ϕ > 0) and is close to the border ϕ = 0. Our proposed method removes noise at A by only using points like B located in the same phase as A, but not C or D which are in different phases and whose intensities are very different than that of A. Contributions of C, D to denoising at A lead to loss of contrast there. In comparison, in the original NLmeans, even though the neighborhoods of A and C are not that close, but a larger h value brings ∥v(A)−v(C)∥ 2 2,a /h 2 closer to 0, causing weight w(A, C) to be large. As a result, pixel C contributes significantly to the smoothing of A, causing decrease of contrast at A.
The paper has the following four main contributions: • Present segmentation boosted NL-means (SNL-means) that utilized segmentation to remove excessive additive Gaussian noise while preserving contrast. Automatic selection of parameter is also addressed for piecewise constant and smooth images respectively. To reduce the computational burden, block-wise implementation of SNL-means is discussed. • Develop transform based NL-means (TNL-means) that preserves more image geometry when segmentation is not reliable due to more excessive noise. Figure 1. Illustration of the motivation of the proposed method, where a level set function ϕ is used to represent phases.
• Provide mathematical validation to block-wise SNL-means from Bayesian perspective and pixel-wise SNL-means from nonparametric regression perspectives. Selection of convolution kernel and weight function is also addressed in detail. • Generalize SNL-means to Poisson, Rician and Ultrasound noise removal. The rest of the paper is organized as follows. Section 2 develops pixel-wise and block-wise SNL-means filter which takes full advantage of segmentation with automatic filtering parameter configuration. TNL-means is also presented to substitute SNL-means when the noise is too excessive to trust segmentation out of the noisy observation or its pre-smoothed variants. Section 3 focuses on mathematical justification of SNL-means in terms of Bayesian inference and nonparametric regression. To address other types of noise including Poisson noise, Rician noise and multiplicative noise, several variants of our SNL-means are stated in Section 4. The numerical implementation and more results are shown in Section 5, followed by conclusion in Section 6.
2. Proposed Scheme. In this section, as done similarly in the original NL-means, we assume that the observed noisy image v is related to the underlying clean image u by: v = u + n with n independent Gaussian noise. Sections 2.1 and 2.2 explain how to incorporate segmentation and geometry extraction transform respectively to remove excessive noise while preserving geometries. Generalization to removal of other noise types will be presented in Section 4.

Segmentation Boosted NL-means filter (SNL-means).
2.1.1. Mutual Position Function. We introduce the mutual position function m(x, y) to describe whether pixels x and y lie in the same phase as follows: Phase information can be obtained from image segmentation. Since segmentation is not the focus of this paper, we just apply Mumford-Shah segmentation method [1]. We use Chan-Vese method [25], which is a level set implementation of Mumford-Shah and works well for piecewise constant images. If there are only two phases in an image, one level set function ϕ is enough to describe them: {x|ϕ(x) > 0} and {x|ϕ(x) < 0} while {x|ϕ(x) = 0} is the object contour. When there are four phases, two level set functions ϕ 1 , ϕ 2 are needed [37]. For instance, {x|ϕ 1 (x) > 0, ϕ 2 (x) > 0} represents one phase. By the four color theorem from graph theory, two level set functions suffice to distinguish all objects in an image. Let H(t) be the Heaviside function defined as then we could define m(x, y) as below: where ϕ(x) is the level set function obtained by the 2 phase Chan-Vese segmentation model and ϕ 1 , ϕ 2 are obtained by the multi-phase Chan-Vese segmentation model.

SNL-means.
With the mutual position function m(x, y), we define a new weight as where the normalizing constant Here v(N (x)) represents the observed image patch in similarity window N (x) centered at pixel x and ∥·∥ 2,α represents the Euclidean distance weighted by Gaussian convolution kernel G α . The segmentation boosted nonlocal means filter is then defined as: S K (x) is a spatial neighborhood of x that contains at least K pixels in the same phase as x. For a fixed pixel x near the boundary, m(x, y) tends to be 0 for some pixel y, leading to insufficient number of contributing pixels if the size of the search window S(x) is fixed (nonadaptive). To ensure similar number of contributing y's in search window of every x, we look for an adaptive window S K (x) as follows. At first we define the candidate based on the segmentation and square windows: and k corresponds to the radius of square windows. To ensure comparing at least K pixels, we set our adaptive search window S K (x) as S l (x) where and |S k (x)| the cardinality of S k (x). If the phase containing x has fewer pixels than K, we define this whole phase as S K (x).

Filtering Parameter h.
Selection of the filtering parameter h plays a vital role in the proposed SNL-means algorithm. Here we describe how to automatically select the parameter h for piecewise constant images and piecewise smooth images separately.
• Piecewise Constant Images For the sake of simplicity, suppose the underlying clean image u : Ω → R 2 is piecewise constant with 2 phases such that one level set function ϕ is enough to describe them. Let where c 1 , c 2 are constants. Under the assumption of additive noise , we have where n i (x) are independent and identically distributed (i.i.d.) within Ω i , i = 1, 2. Using the uniformly minimum variance unbiased (UMVU) estimator [26], we get the best estimator for noise variance in each phase as Then the filtering parameter h i used in the i-th phase (i = 1, 2) is chosen as h i = √ Var(n i ). In the multi-phase case, we still employ the unbiased estimation of the noise variance in each phase to obtain the phase-wise filtering parameter h's.

• Piecewise Smooth Images
Let the underlying image u(x) be piecewise smooth with independent and spatially variant noise n(x) where n(x) has mean zero and variance σ 2 (x). By the assumption, the mean and variance for the conditional probability are as below: In this case, there are no closed forms for the mean and variance and thereby we employ the nonparametric dual model to estimate the noise level locally [27]: at pixel x, we use the nearest four-pixel neighboring set, denoted by N 4 (x), as its sample space and obtain the Gasser-Sroka-Jennen-Steinmetz pseudo-residual ) which guarantees that E(ϵ(x) 2 ) = σ 2 (x). By taking the average to correct the bias, we define the unbiased local variance estimator as: where β is between 0 and 1 and is set as .5 in our experiments and the results are not sensitive to the choice of β.
The key ingredient in the proposed method SNL-means is the mutual position function defined from segmentation. The improvement over the original NL-means is significant when the segmentation is accurate. In Figure 2, we compare the results of NL-means and SNL-means (fed with ground truth segmentation) on a noisy MRI brain image with intensity range [0,230] and additive zero-mean Gaussian noise with 10% (σ = 23) noise. It is observed that the SNL-means leads to sharper reconstructed image or anatomically clearer separation of different phases (gray matter, white matter and cerebral spinal fluid).

Transform based NL-means filter.
In practice, one does not have the exact segmentation of the ground truth but only the noisy observation. When the observed noisy image has excessive or spatially variant noise, it is more difficult to distinguish different phases. As the noise level grows, image segmentation gets more and more challenging: fine features are usually missed and some artificial edges are created. In this case, we have the following three options: • Pre-smooth. Pre-smooth the noisy image to get a cleaner image, on which we apply segmentation. This approach is the easiest and works only when the noise level is low. But the segmentation depends on the selection of the pre-smooth approach and is not accurate when noise is excessive. Most of the existing denoising approaches tend to loose salient features while removing excessive noise. Loss of salient features in the pre-smoothed image causes segmentation errors. • Iterative denoising. Alternately segment and smooth images using SNLmeans in an iterative way: apply segmentation on the noisy observation, then apply segmentation boosted NL-means on the original noisy image to get one cleaner image, which is segmented again for better phase information that in the next round helps remove noise from the original noisy image. Repeat the procedure until the adjacent smoothing results are close. This method is intuitive but expensive due to iterations. • Geometrical structure extraction transform. Substitute segmentation by some more general transform that is less sensitive to noise but detect more intrinsical geometric structures. Then replace mutual position function defined on segmentation by some new function (again part of the weight definition) describing geometric structure similarity. Segmentation only detects edges, while more general geometric structure extraction transforms detect more geometries including edges. In what follows, we explain this method in detail.
The idea is to apply certain transform catching the main features of the image and to reduce the impact of noise during the restoration step. We demonstrate the idea using contourlet transform, but any other general transform could be used. Contourlet is a pyramidal directional transform which represents C 2 natural images efficiently using contour segments [28]. Due to the contour-like edges in many real world images, the contourlet transform is appropriate to be incorporated into the computation of the local weight. The proposed transform based NL-means (TNLmeans) scheme is with weight w(x, y) the product of weighted intensity similarity and geometry similarity defined as below where C(x) is the normalization factor, η 1 , η 2 are balancing constants and h 1 , h 2 are respective filtering parameters. In our later experiments, h 1 and h 2 are set as the standard deviation of the noise while η 1 , η 2 are adjustable values within (0, 1). Here T is designated as the combination of three steps: contourlet transform, thresholding in the transform domain, and inverse contourlet transform back to the image domain. In Figure 3 and Table 1 we demonstrate that TNL-means outperforms thresholding based contourlet approach and NL-means respectively. We add the same additive Gaussian noise with σ = 23 as in Figure 2 to the clean MRI brain image. We can see in Figure 3 that the contourlet restored image after thresholding in the transform domain captures certain geometries of the original image while having much noise remained (middle left). The proposed TNL-means method takes advantage of the geometric information to build the weight and removes the noise more efficiently (right two images). In Table 1, the comparisons with NL-means for three different noise levels 10% (σ = 24.6), 15% (σ = 36.9), 20% (σ = 49.3) are listed. We use three quantitative measures: peak-to-noise ratio (PSNR), root of mean squared error (RMSE) and structural similarity (SSIM) [29]. It can be seen that the proposed method consistently gives the highest PSNR, lowest RMSE and the highest SSIM. The higher the noise level, the more the TNL-means outperforms.
In practice, depending on the structure of the image, contourlet transform based weight may not be the most appropriate candidate, and other filter bank based transforms may be more applicable, e.g. curvelet transform [20] and shearlet transform [21].  In this subsection, we present the block-wise SNLmeans filter. For natural images with large size, pixel-wise SNL-means filtering will cost too much time on comparing the weighted patches. To lower the computation burden, we follow the idea in [22] to develop block-wise SNL-means. The idea is to use SNL-means to obtain estimate of the clean image in a block by taking the weighted average of the most similar blocks and then get the recovered intensity at a specific pixel x by averaging over all the estimators of intensities at x given by all the blocks containing x. Rather than choosing all the blocks in the search window, we requires the centers of the contributing blocks to lie in the same phase as that of the target block. Specifically, we first partition the whole image domain Ω into overlapping blocks Then we combine the information from segmentation to the nonlocal means formula: where and C(B i ) is a normalizing constant. Here M (B i , B j ) is a block-wise mutual position function and for simplicity we only apply the centers of the two blocks of interest: where x i and y j are the center pixels of blocks B i and B j .
Let B(x) = {B i |x ∈ B i ⊂ Ω} be the set of all the blocks including x, the estimated intensity at each pixel is calculated as below: However, there is a trade-off between the denoising performance and the computation complexity. Block-wise implementation will lower the computation time tremendously but at the cost of the quality which will be shown in Section 5 Table  7.
3. Mathematical validations. In this section, we provide mathematical validations, in more general setup, of the proposed block-wise SNL-means from Bayesian perspective (Section 3.1) and the pixel-wise SNL-means version using non-parametric regression (Section 3.2). Based on the validations, we generalize the proposed scheme to deal with images with other noise types, e.g., Poisson, Rician and Ultrasound noise in Section 4.

Bayesian Interpretation of Block-wise SNL-means.
In [14], Kervrann et al. derived the NL-means filter in perspective of Bayesian framework which states the relationship between the weight function in NL-means filter and the probability density function (pdf) of the white Gaussian noise. We extend their derivation to general additive noise. Various local convolution kernels relating the structural prior of images are also considered.
As before, we assume v = u + n where n is any general independent noise and u is piecewisely constant or smooth. Let {Ω i } i be a partition of the image domain Ω such that For every x ∈ Ω, we denote the vectorized u and v taken in a neighborhood of x Letû be an estimator for u, and define the expectation of loss associated with this estimator as where Λ ⊆ R k is the range of u(x).
By choosing L(u(x),û(x)) = |u(x) −û(x)| 2 and letting p be a probability, the optimal estimator is computed as below: Let S K (x) be the neighborhood of x defined as in Section 2.1.2. By the assumption that patches with centers within the same phase follow the same distribution, the neighboring patches with centers in S K (x) can be considered as the realizations of the patch u(x). Assuming u(x) is uniformly distributed in {u(y) | y ∈ S K (x)} which replaces Λ, i.e., there is no preference to choose one neighboring patch from another, we have As K → ∞ and thereby the number of the neighboring patches increases, the law of large numbers asserts: ∑

p(v(x)|u(x))p(u(x)).
Then one gets an approximated u opt as below:
Considering v(y) as the realization of u(y), we get an estimation of u in the similarity window N (x) Next, we assume n is a Gaussian random field and n(x) is a multivariate Gaussian random variable therein. Note that v(x) = v(y) + n(x, y), with n(x, y) = u(x) − u(y) + n(x) − n(y), (y ̸ = x).
If x and y are very close to each other within the same phase, then u(x) ≈ u(y). In addition, to embody certain regularity of the original image between u(x) and u(y), we assume that n(x, y) = R −1 E where E ∼ N (0, σ 2 I) and R ∈ R k×k is symmetric positive definite. Thus we get ) .

A GEOMETRY GUIDED IMAGE DENOISING SCHEME 11
It's obvious that Var(n(x, y)) = (R T R) −1 and ∥R(v(x) − v(y))∥ 2 is a semi-metric measuring the closeness of two image patches. In the original NL-means, R T R is chosen as a square matrix whose diagonal is the vectorized 2D Gaussian convolution kernel G α without taking phase information into account. In this case, equation (6) leads to the same formula as (5) in Section 2.3. The diagonal entries in R T R could come from Gaussian convolution kernels or more generally ones like rectangular disk, circular disk, average (constant) as shown in Figure 4. We have tested some piecewise constant images and piecewise smooth images using the listed kernels and found that for piecewise constant images, flatter kernels enhance denoising performance. However for most images with more intensity variations, more exponential-like kernels are preferred. The denoising assessment for each kernel in terms of PSNR for a 2-phase toy image and a T1 MRI image are shown in Table 2.

Nonparametric Regression Understanding of Pixel-wise SNL-means.
Next we understand the proposed pixel-wise SNL-means filter in terms of local smooth regression. Define adaptive search window S K (x l ) for a fixed pixel x l as in Section 2.1.2. To get an estimate of the clean image intensity at x l , we use sample pairs composed of intensities at each pixel in the search window and those in their similarity windows, i.e., we use sample pairs , and Y i ∈ R k−1 the observed column-wise stacked intensity vector of the neighbors (similarity window of size √ k × √ k) of pixel x i except at x i . The Nadaraya-Watson kernel-weighted average [30] provides an estimate of u at x l asû where the weight function Here R is a diagonal matrix defined as in the previous section and D(·) is a kernel function satisfying the conditions: (i) positivity: D(t) ≥ 0 (ii) decay as distance increases: lim t→∞ D(t) = 0 and lim t→0 D(t) = 1.
(iii) unit volume: the above local linear regression can be actually obtained by solving a weighted least squares problem: which is the discretization of In the original NL-means algorithm, the weight function is defined through selecting D(t) as exponential function e −t 2 and R is based on discretized Gaussian convolution kernel. But different combinations of kernel function D and convolution kernel lead to different weight functions. Selection of convolution kernel has been addressed in the previous section. Next, we focus on selection of kernel function D on removing additive Gaussian noise. Table 3 lists five kernel functions D without normalizing constant and Figure 5 shows their respective graphs. Note D applied on gives weights before normalization. For a fixed h, the principle of selecting D is as follows: for piecewise constant images, choose D that is flatter and more spread out, e.g. Cauchy; while for images with large intensity variation in a small neighborhood, choose D that is more steep in the range of h , e.g., exponential one if the range is inside [1, 1.5]. The reason is that for piecewise constant images, the intensity is close to one constant for all pixels in the same phase. Combining with the fact that larger number of contributing pixels leads to a more robust estimate (for i.i.d. random variables r i 's, the variance of ∑ N i=1 r i /N approaches to 0 as N goes to infinity), we expect most of the weights to be nonzero and close to each other, which can be done by choosing a kernel function D that is flat. From the graphs in Figure 5, it is obvious that Cauchy function is the flattest and is thus the best choice for piecewise constant images. On the other hand for images with various intensities, we only expect pixels with close intensities to contribute more. In terms of weights, we only expect those with smaller ∥R(Y l − Y j )∥ 2 (h is fixed) to get larger values. Therefore, a kernel that is steep in the range of is preferred. We test a piecewise constant toy image (h fixed as σ = 18.2) and a T1 MRI image (h fixed as σ = 25.5) with different weight functions (fixing the rectangular disk as R), and we list the results as in Table 4. For the piecewise constant toy image, Cauchy function shows significantly better performance than the others. While for  the T1 MRI image that has more intensity variations, the range of [1, 1.5], so as analyzed above, the result using exponential function is the best. Furthermore, the similar conclusion will be obtained in the case of Rician noise with the same test images. 4. Generalization to other noise types. Section 2 addresses in detail how SNLmeans removes Gaussian noise automatically. To deal with non-Gaussian noise, we provide two options here. One is to insert a specific noise type into the generic Bayesian formula (6) in Section 3.1 and see whether a closed form is computable. For some distributions (e.g., Poisson, binomial and negative binomial), another option is to transform the given noise to i.i.d. Gaussian noise, then apply the SNL-means scheme thoroughly discussed in Section 2 to remove the additive Gaussian noise, and finally transform the processed image back to the original image space. But special care must be taken to correct the bias resulting from the inverse transform.
the image contaminated by Poisson noise is transformed to an image with approximately homogeneous additive Gaussian noise. In the transformed domain, our proposed method (Section 2.1) with estimation of the unknown additive noise level (Section 2.1.3) can be implemented. To minimize the undesired bias generated by the transformation, the asymptotic unbiased inverse transform is used In a word, our proposed denoising scheme on Poisson noise removal iŝ 4.2. Rician noise. In practice, images may be complex valued and medical technicians visualize the magnitude of the images. If the noise is well modeled by a Gaussian probability density function in the real and imaginary parts of the complex data respectively, then the noise in the associated magnitude image can be modeled as Rician noise. This noise model is used in a wide variety of data, such as MRI, SAR [32]. Let v ∼ R(u, σ) be the observed magnitude image obeying Rician distribution with parameters u ≥ 0 and σ ≥ 0, where J n represents the modified Bessel function of order n, This distribution reduces to a Rayleigh distribution if u = 0. From the expression of J n (x), it's easy to see that A GEOMETRY GUIDED IMAGE DENOISING SCHEME 15 The mean and variance of v therefore are, respectively, where ] .
The second moment is therefore implying a denoising scheme for remove Rician noise. Based on the observation that w(x, y) = 1 and the nonlocal estimator proposed by Wiest-Daesslé et al. [19], we present the proposed scheme on Rician noise removal 4.3. Ultrasound noise. In ultrasound imaging, the observed noisy image v is modeled as v = u + √ u · n [33,34] where u is the underlying clean image, and noise is √ u · n with n ∼ N (0, σ 2 ), i.e., the noise depends on u. Then v|u ∼ N (u, uσ 2 ) i.e., ) .
Following the reasoning in Section 3.1 and plugging the above p(v(x)|u(x)) into the generic Bayesian formula (6), we get where ./ represents the entry-wise division and h can be chosen around σ.

More Experimental Results.
Besides those results shown in the previous sections, here we further test our algorithm and its variants on a variety of scalar images and also compare these results with those of the standard NL-means algorithm and the other state-of-the-art denoising methods. Equipped with the geometric information obtained from segmentation, appropriate convolution kernel and weight function based on the image to be denoised, the proposed framework shows the prominence on edge preservation and noise removal. All tests were run under Windows 7 Professional and MATLAB v7.12 (R2011a) on a laptop with an Intel Core 2 Duo CPU at 2.8 GHz and 3 GB of memory. Both qualitative and quantitative comparison are provided throughout this section. For quantitative comparison, among many other criteria, two assessment measures are used here. One is peak signal-to-noise ratio (PSNR) which is defined in decibels (dB) as PSNR = 10 log 10 255 2 MSE where MSE is a mean squared error between the noisy image and the clean image. This is the main criterion used later to compare the denoising performance. Another one is signal to noise ratio (SNR) which quantifies the recovery by comparing between the "ground truth" u(x) and the reconstructed imageû(x). We adopt the SNR formula in [38].
It's used only for assessing the performance of the ultrasound denoising. Four sets of tests are done to show the significant improvements of our proposed algorithms. Note that if not specified explicitly, pixel-wise implementation is always performed and the clean images for our tests range in [0, 255].
• The first set is to show the advantage of the pixel-wise SNL-means on piecewise constant and smooth image denoising. We start with a two-phase piecewise constant image containing three shapes contaminated by white Gaussian noise with different known variances. We apply our proposed SNL-means filter (2), the original NL-means algorithm, total variational method (ROF) [3], block-matching and 3D filtering (BM3D) [35] and K-SVD [36] to the same image. Optimal parameters are selected for all. In the proposed method, filtering parameter h is set as the given noise level, and since it is a piecewise constant image, as discussed in Section 3.2, we replace the exponential weight function used in the original NL-means by the Cauchy weight function. Table 5 compares in terms of PSNR as noise variance changes. It can be seen that across all the noise levels, the proposed method consistently and significantly performs better than all the other methods. Figure 6 shows the visual comparisons when the noise level is 50. Next, we test a multi-phase piecewise constant Shepp-Logan phantom image. Multiple phases are obtained by using multi-phase Chan-Vese segmentation model [37]. To fetch more accurate phase information, a pre-smoothed version of the input image is favored for segmentation. Here BM3D is applied for pre-smoothing. We employ our 6-phase SNL-means filter (2) with phase-wisely defined h (Section 2.1.3). From Figure 7 and Table 6, one can see that the proposed model significantly outperforms the others both quantitatively and qualitatively. In Figure 8, we test a brain MRI image corrupted by inhomogeneous additive Gaussian noise using iterative 3-phase piecewise smooth denoising scheme mentioned in the beginning of Section 2.2. Unlike the NL-means using one universal filtering parameter h, our proposed method mentioned in Section 2.1.3 estimates the inhomogeneous noise level automatically and defines the spatially variant h accordingly. The implementation starts from the NL-means filtering with low h and ends after two iterations of segmentation-smoothing. One can see the original NL-means removes some fine features while iterative denoising removes mostly random noise.
• The second test is to illustrate the quality/speed balance for piexel-wise and block-wise SNL-means. We test using the 2-phase toy image in Figure 6 with 10% noise (σ = 27.3). There are two factors in block-wise SNL-means that directly affect its quality and speed. One is the block size described by its radius r and the other one is the distance between two adjacent blocks denoted by d. We vary both r and d to get different block-wise SNL-means results and compare their PSNR and computation time (in seconds) with that of SNL-means. The comparison is listed in Table 7, in which pixel-wise and block-wise SNL-means are abbreviated as SNL and BSNL respectively. The two numbers behind "BSNL" represent r, d respectively. It is not surprising to see that pixel-wise SNL-means takes more time than all the block-wise SNL-means but gives the best quality. And as r or d increases, block-wise SNL-means is faster but pays the price of lower recovery quality.  (7) with NL-means based filter (i.e., u =Ã −1 (N L(A(v)))) in Figure 9. Even though the PSNR improvement is not much, a close visual comparison shows that the proposed method provides better contrast, especially the buildings in the background. Second row: inhomogeneous noise added; difference from the noisy image by NL-means and our method.
Next we add Rician noise on an eagle image as described in Section 4.2 and test our proposed method (8). After segmenting the image into two phases, we applied exponential and Cauchy weight functions for the phases inside and outside of the eagle respectively. As seen from Figure 10 our proposed method results in more noise suppression in the background and clearer outline of the eagle.
For removing speckle noise in ultrasound images, we simulate the noisy image in the way as described in Section 4.3. Results of our proposed method (9) and NL-means based ultrasound denoising filter are shown in Figure  11. Since Bayesian-based NL-means filter shows better denoising performance than some competitive state-of-the-art methods: Lee's filter, Kuan's filter, SPF and SRAD [38], we restrict the comparison to in between NL-means filter and our method. Consistently, the proposed method leads to a reconstruction closer to the ground truth.
• Finally, we apply the SNL-means filter on a real clinical data. Figure 12 shows the result of SNL-means on a scanning electron microscopy (SEM) image of mouse axons showing myelin sheath, mitochondria and microtubules. The sharp denoising result makes it easier to investigate the mitochondria (amount and structure) and microtubules which are affected by multiple sclerosis (MS).
6. Conclusion. The proposed method extracts reliable geometric information and incorporates it to guide the denoising. In the framework of NL-means, we described how to remove excessive inhomogeneous noise while preserving fine features. To speed up computation, the blockwise implementation of our proposed filter is described as well. We also adapted the idea to remove excessive noise of various types. It appears that coupling geometry extraction and nonlocal means leads to better denoising from both quantitative and qualitative aspects. In spite of getting stateof-the-art recoveries, we believe that there is large room for further improvements by exploiting the structures of edges, as well as other image geometry, in a more effective way.