# Spatial Images Feature Extraction Based on Bayesian Nonlocal Means Filter and Improved Contourlet Transform

Authors: Pengcheng Han; Junping Du
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

Spatial images are inevitably mixed with different levels of noise and distortion. The contourlet transform can provide multidimensional sparse representations of images in a discrete domain. Because of its filter structure, the contourlet transform is not translation-invariant. In this paper, we use a nonsubsampled pyramid structure and a nonsubsampled directional filter to achieve multidimensional and translation-invariant image decomposition for spatial images. A nonsubsampled contourlet transform is used as the basis for an improved Bayesian nonlocal means (NLM) filter for different frequencies. The Bayesian model adds a sigma range in imagea priori operations, which can be more effective in protecting image details. The NLM filter retains the image edge content and assigns greater weight to similarities for edge pixels. Experimental results both on standard images and spatial images confirm that the proposed algorithm yields significantly better performance than nonsubsampled wavelet transform, contourlet, and curvelet approaches.


---

## Body

## 1. Introduction

In spatial rendezvous and docking, spatial images are obtained by multisource remote sensors. Spatial images are inevitably mixed with different levels of noise and distortion. The accurate image feature extraction will be helpful for spatial object recognition and can directly influence the success of spatial rendezvous and docking [1, 2]. Image feature extraction of spatial images is based on the definition of image features; to some extent, it can be said that it is based on sensitivity changes to image grayscale values for the human eye. Multidimensional image representation can process images for the sparsest representation, especially for 2D image signals [3, 4]. This approach identifies optimal high-dimensional function representation for an image and yields superior image-processing results for an effective solution. A nonlocal means (NLM) filter uses redundant image information on the basis that structural similarity superimposed on pixel noise is random and noise can be effectively removed using weighted averages [5, 6]. Compared to traditional statistical filtering methods, NLM filtering overcomes the constraint of the local neighborhood and extends pixel similarity to block-based similarity, so it is very suitable to deal with spatial images.In this paper, we use a nonsubsampled pyramid structure and a nonsubsampled directional filter to achieve multidimensional and translation-invariant image decomposition for spatial images. A nonsubsampled contourlet transform is used as the basis for an improved Bayesian nonlocal means (NLM) filter for different frequencies. The Bayesian model adds a sigma range in imagea priori operations, which can be more effective in protecting image details. The NLM filter retains the image edge content and assigns greater weight to similarities for edge pixels. Experimental results both on standard images and spatial images confirm that the proposed algorithm yields significantly better performance than nonsubsampled wavelet transform, contourlet, and curvelet approaches.The rest of this paper is organized as follows. Section2 describes multidimensional image decomposition, with a focus on contourlet and nonsubsampled contourlet transforms (NSCTs). Section 3 outlines application of an NLM filter and proposes an improved NLM algorithm based on a Bayesian model. Section 4 applies the improved NLM filter to NSCT, especially NSDFB, to process image features for further extraction. Section 5 compares feature extraction results for the proposed algorithm and other algorithms. Section 6 concludes the paper.


## 2. Contourlet Transform Decomposition

### 2.1. Multidimensional Image Decomposition

The target of image multidimensional representation is to provide a description of image with less characteristic information. The wavelet transform is a classic image multidimensional representation algorithm that has a good effect on image edge points [7, 8]. However, the wavelet transform can capture only limited direction information in the horizontal, vertical, and diagonal directions, as shown in the left side of Figure 1. It is difficult to express image smoothness contours; a better image representation is shown in the right side of Figure 1.Figure 1
Multidimensional image decomposition.Other well-known multidimensional image decomposition algorithms include bandlets, brushlets, edge multidimensional transform, complex wavelets, and wedgelet. However, these algorithms require image edge detection and then summarize a representative adaptive coefficient. A decomposition algorithm that can transform an image into fixed decomposition coefficients is desirable. These coefficients can then be used in a broader context that does rely on edge detection alone but also includes better directional image decomposition.In 2004, Candès and Donoho proposed a curvelet transform that uses a value approximation algorithm for a continuous 2D spatial domain and adds a smooth signal on the basis of a 1D Fourier transform [9]. The best approximation deviation is O((logM)3M-2) for curvelet and O(M-1) for wavelet transforms. The curvelet transform is first applied to a continuous signal and then combines a multidimensional filter and ridgelet transformation. A second curvelet transform is based on frequency segments and extreme judgment. The curvelet transform is universally applicable to continuous signals, but there will be parallel noise in discrete fields [10]. It is also biased in directional image decomposition. The reason is that the typical rectangular sampling mode leads to a priori geometric deviation in decomposition of discrete image signals, especially in the horizontal and vertical directions. This limitation prompted researchers to develop a new multiscale decomposition algorithm that does not depend on edge detection and can decompose images in cross-scale multidimensions.


### 2.2. Contourlet Transform

The contourlet transform is a multidimensional decomposition algorithm proposed by Do and Vetterli in 2005 [11]. The transform can be directly used for multidimensional decomposition of discrete image signals. It has a dual filter for image decomposition and yields a smoother sparse representation of the original image. The two filters are a Laplacian pyramid (LP) filter [12, 13] and a directional filter bank (DFB) [14]. The LP yields nonconsecutive image points, and then the DFB connects consecutive points into a nonlinear structure. The process is shown in Figure 2.Figure 2
Contourlet transform.A subsample contourlet transform uses a relevance factorM for image subsampling at each decomposition level. A 2D filter is evolved from the 1D filter. For complete image reconstruction, the following relationship holds for the 1D filter:(2.1)M0(z)N0(z)+M1(z)N1(z)=2,Sm={D(2l-1,2),for0≤m≤2l-1-1,D(2,2l-1),for2l-1≤m≤2l-1,
where M0(z) and M1(z) represent low- and high-pass analysis filters, and N0(z) and N1(z) represent low- and high-pass synthesis filters, respectively. Downsampling matrices Sm are shown above. For 2D complete decomposition,(2.2)M0(M′(z))N0(M′(z))+M1(M′(z))N1(M′(z))=2.M′(z) represents discrete properties of the heterogeneous domain, which can help to reduce the filter complexity from O(N2) to O(N). In two dimensions, the first DFB step is construction of the spectrum in the frequency domain using two-channel quincunx filter banks to decompose an image into horizontal and vertical directions [15]. The DFB equivalent parallel family is(2.3){Dk(l)[n-Sm(l)]}0≤m≤2l,m∈L2.Cutting operations on both directions for the decomposition spectrum provide 2D directional and segmental image decomposition. Like the discrete wavelet transform, the discrete downsampling contourlet transform is shift-invariant [16].


### 2.3. Nonsubsampled Contourlet Transform

NSCT is a fast implementation of the contourlet transform that provides a shift-invariant and multidimensional image representation [17]. Compared with subsampled contourlet transforms, NSCT is closer to the nonredundant wavelet transform [18]. NSCT uses a 2D nonsubsampled filter bank and can be expressed as(2.4)M0(z)N0(z)+M1(z)N1(z)=1,
where M(z) represents a 2D filter of the z transform, M0(z) and M1(z) represent 2D low- and high-pass analysis filters, and N0(z) and N1(z) represent 2D low- and high-pass synthesis filters, respectively. There are also other limitations for the filter design.NSCT involves two steps: multidimensional representation and directional decomposition. Multidimensional representation is achieved by nonsubsampled pyramid decomposition. This step is similar to the 1D discrete nonsubsampled wavelet transform (NSWT), which uses the à trous method [19]. Compared to NSWT, NSCT uses a nonsubsampled 2D filter. The frame bound of an NSCT directional decomposition is(2.5)P1≤|M0(eε)|+|M1(eε)|︸t(eε)≤P2,P1=essinf⁡t(eε),P2=esssup⁡t(eε),ε∈[-π,π]2.After decomposition of the first layer, the sampling filter banks provide multiscale decomposition of the underlying properties. The process for two-layer nonsubsampled pyramid decomposition is shown in Figure3.Figure 3
Non-subsampled pyramid decomposition.The frequency domain for layerj supported by a low-pass filter is [-(π/2j),-(π/2j)]; the replacement domain is from [-(π/2j-1),-(π/2j-1)] to [-(π/2j),-(π/2j)], which is supported by a high-pass filter.Each step in NSWT image decomposition involves three directions. The total image redundancy is3J+1; in NSP, the result redundancy is  J+1 [20]. The second NSCT step provides directional information via the nonsubsampled filter, which combines two-channel quincunx sampling filters and a resampling operation for 2D frequency division on directional edges [21]. More accurate directional details can be sampled discretely on a sample stage. Sampling uses a quincunx matrix Q and considers image direction alignment. The process is shown in Figure 4.Figure 4
Nonsubsampled directional filter banks.


## 2.1. Multidimensional Image Decomposition

The target of image multidimensional representation is to provide a description of image with less characteristic information. The wavelet transform is a classic image multidimensional representation algorithm that has a good effect on image edge points [7, 8]. However, the wavelet transform can capture only limited direction information in the horizontal, vertical, and diagonal directions, as shown in the left side of Figure 1. It is difficult to express image smoothness contours; a better image representation is shown in the right side of Figure 1.Figure 1
Multidimensional image decomposition.Other well-known multidimensional image decomposition algorithms include bandlets, brushlets, edge multidimensional transform, complex wavelets, and wedgelet. However, these algorithms require image edge detection and then summarize a representative adaptive coefficient. A decomposition algorithm that can transform an image into fixed decomposition coefficients is desirable. These coefficients can then be used in a broader context that does rely on edge detection alone but also includes better directional image decomposition.In 2004, Candès and Donoho proposed a curvelet transform that uses a value approximation algorithm for a continuous 2D spatial domain and adds a smooth signal on the basis of a 1D Fourier transform [9]. The best approximation deviation is O((logM)3M-2) for curvelet and O(M-1) for wavelet transforms. The curvelet transform is first applied to a continuous signal and then combines a multidimensional filter and ridgelet transformation. A second curvelet transform is based on frequency segments and extreme judgment. The curvelet transform is universally applicable to continuous signals, but there will be parallel noise in discrete fields [10]. It is also biased in directional image decomposition. The reason is that the typical rectangular sampling mode leads to a priori geometric deviation in decomposition of discrete image signals, especially in the horizontal and vertical directions. This limitation prompted researchers to develop a new multiscale decomposition algorithm that does not depend on edge detection and can decompose images in cross-scale multidimensions.


## 2.2. Contourlet Transform

The contourlet transform is a multidimensional decomposition algorithm proposed by Do and Vetterli in 2005 [11]. The transform can be directly used for multidimensional decomposition of discrete image signals. It has a dual filter for image decomposition and yields a smoother sparse representation of the original image. The two filters are a Laplacian pyramid (LP) filter [12, 13] and a directional filter bank (DFB) [14]. The LP yields nonconsecutive image points, and then the DFB connects consecutive points into a nonlinear structure. The process is shown in Figure 2.Figure 2
Contourlet transform.A subsample contourlet transform uses a relevance factorM for image subsampling at each decomposition level. A 2D filter is evolved from the 1D filter. For complete image reconstruction, the following relationship holds for the 1D filter:(2.1)M0(z)N0(z)+M1(z)N1(z)=2,Sm={D(2l-1,2),for0≤m≤2l-1-1,D(2,2l-1),for2l-1≤m≤2l-1,
where M0(z) and M1(z) represent low- and high-pass analysis filters, and N0(z) and N1(z) represent low- and high-pass synthesis filters, respectively. Downsampling matrices Sm are shown above. For 2D complete decomposition,(2.2)M0(M′(z))N0(M′(z))+M1(M′(z))N1(M′(z))=2.M′(z) represents discrete properties of the heterogeneous domain, which can help to reduce the filter complexity from O(N2) to O(N). In two dimensions, the first DFB step is construction of the spectrum in the frequency domain using two-channel quincunx filter banks to decompose an image into horizontal and vertical directions [15]. The DFB equivalent parallel family is(2.3){Dk(l)[n-Sm(l)]}0≤m≤2l,m∈L2.Cutting operations on both directions for the decomposition spectrum provide 2D directional and segmental image decomposition. Like the discrete wavelet transform, the discrete downsampling contourlet transform is shift-invariant [16].


## 2.3. Nonsubsampled Contourlet Transform

NSCT is a fast implementation of the contourlet transform that provides a shift-invariant and multidimensional image representation [17]. Compared with subsampled contourlet transforms, NSCT is closer to the nonredundant wavelet transform [18]. NSCT uses a 2D nonsubsampled filter bank and can be expressed as(2.4)M0(z)N0(z)+M1(z)N1(z)=1,
where M(z) represents a 2D filter of the z transform, M0(z) and M1(z) represent 2D low- and high-pass analysis filters, and N0(z) and N1(z) represent 2D low- and high-pass synthesis filters, respectively. There are also other limitations for the filter design.NSCT involves two steps: multidimensional representation and directional decomposition. Multidimensional representation is achieved by nonsubsampled pyramid decomposition. This step is similar to the 1D discrete nonsubsampled wavelet transform (NSWT), which uses the à trous method [19]. Compared to NSWT, NSCT uses a nonsubsampled 2D filter. The frame bound of an NSCT directional decomposition is(2.5)P1≤|M0(eε)|+|M1(eε)|︸t(eε)≤P2,P1=essinf⁡t(eε),P2=esssup⁡t(eε),ε∈[-π,π]2.After decomposition of the first layer, the sampling filter banks provide multiscale decomposition of the underlying properties. The process for two-layer nonsubsampled pyramid decomposition is shown in Figure3.Figure 3
Non-subsampled pyramid decomposition.The frequency domain for layerj supported by a low-pass filter is [-(π/2j),-(π/2j)]; the replacement domain is from [-(π/2j-1),-(π/2j-1)] to [-(π/2j),-(π/2j)], which is supported by a high-pass filter.Each step in NSWT image decomposition involves three directions. The total image redundancy is3J+1; in NSP, the result redundancy is  J+1 [20]. The second NSCT step provides directional information via the nonsubsampled filter, which combines two-channel quincunx sampling filters and a resampling operation for 2D frequency division on directional edges [21]. More accurate directional details can be sampled discretely on a sample stage. Sampling uses a quincunx matrix Q and considers image direction alignment. The process is shown in Figure 4.Figure 4
Nonsubsampled directional filter banks.


## 3. NLM Filter Based on a Bayesian Approach

### 3.1. NLM Filter

Different frequency components play different roles in an image structure. Low-frequency components account for most image energy, forming basic local gradation areas, but play a small role in image content or structure. High-frequency components form the main image edges and determine its basic content or structure and are thus the most important components. Changes in high-frequency information lead to changes in the basic image content or structure, and information extracted from the image by the human eye will thus be subject to major changes. Thus, high-frequency components play the most important role in image perception by the human eye.At present, many image filters only consider adjacent pixels; some filters take into account information for neighboring pixels, such as Yaroslavsky neighborhood filters [22] and bilateral filters [23]. A nonlinear filter involves additive white noise and can effectively handle image redundancy [24].The NLM algorithm takes advantage of grayscale image redundancy and structural redundancy through a weighted average of pixel values to estimate the current pixel value. The value of each pixel is calculated using the Gaussian-weighted Euclidean distance between subblocks; a pixel as taken as the right center of the corresponding subblock. This ensures that pixels with a similar structure are assigned greater weight. For an original imagev={v(i)∣i∈I}, the expression for the image processed using the NLM filter, NL(v), is(3.1)NL(v,i)=∑j∈Iw(i,j)v(j),
where w(i,j) is the Gaussian-weighted Euclidean distance between pixels i and j, which represents the similarity of the image subblocks with i and j as centers:(3.2)w(i,j)=1Z(i)exp⁡[-‖v(Ni)d×d-v(Nj)d×d‖2h2],Z(i)=∑jexp⁡[-‖v(Ni)d×d-v(Nj)d×d‖2a2h2],h2=1σ2(1+r2)∑iN∑jN(rNi-Nj)2,
where Z(i) is a normalization factor, v(Ni) is a set of subblocks with pixel i as the center, ∥·∥2a is a Gaussian-weighted Euclidean distance function, a is Gaussian kernel standard deviation, and h is a filter parameter that controls the degree of smoothing. The weight w(i,j) meets 0≤w(i,j)≤1, if neighborhood pixels are more similar with v(Ni), the weight of center pixel is greater.


### 3.2. NLM Filter Combined with a Bayesian Method

The traditional NLM algorithm is very similar in both grayscale and structure content for smooth neighborhood areas. The algorithm yields the best results in flat areas, where better denoising effects can be obtained. At image edges and in texture-rich regions, the algorithm performs poorly because these regions have many repeat structures, the difference in grayscale content is greater, and the larger Euclidean distance makes the weights very small and reduces denoising capability, especially the ability to retain image detail [25]. To improve the edge retention capacity of the NLM algorithm, a Bayesian algorithm was added to make use of image edge information and adjust the similarity of the neighborhood structure so that the center pixel of edge contents that are similar can be given greater weight. This provides a more effective approach for protecting image detail.The Bayesian NLM filter is expressed as(3.3)n(x)=∑y∈Δ(x)p(v(x)∣u(y))p(u(y))u(y)∑y∈Δ(x)p(v(x)∣u(y))p(u(y)),
where v(x) represents noise data, u(x) represents nonnoise image data, n(x) is the average pixel precision weight for gray value u(y) at a rate of change Δx, and p(v(x)∣u(y))p(u(y)) is the similarity between v(x) and u(y). Equation (5) can then be rewritten as(3.4)w(i,j)=1Z(i)exp⁡[-‖v(Ni)d×d-v(Nj)d×d‖2a2h2-‖n(xi)d×d-n(xj)d×d‖2a2h′].


### 3.3. Improved Bayesian NLM Filter

We propose an improved Bayesian NLM filter in which a sigma range is added to the prior image operation for more effective protection of image details [26]. The first step is analysis of the probability density for pixel levels, which takes the average variance for the improved Bayesian filter. Considering the independence and integrity of an image, its conditional probability distribution can be expressed as(3.5)p(v(x)∣u(y))=∑m=1M×Mp(vm(x)∣um(y)),
where vm(x) and um(y) are the image probability densities at pixel m and vm(y) is a subset of um(y).For L-level image decomposition, the conditional probability density function is(3.6)p(vm(x)∣um(y))=vm(x)L-1F(L)(Lum(y))exp⁡(-Lvm(x)um(y)),(3.7)F(L)=Ga(L)∑lLexp⁡(um(y)l-1-vm(x)l-1).Because of multiscale features, we assume that the prior probabilityp(u(y)) is continuous and uniform, p(u(y))=1/|Δx|. The proposed algorithm uses an iterative technique, which takes the observed value v(x) as the initial value u(y). This treatment can process data directly but takes longer, and details can become fuzzy. If the frequency window is too large, the result will have too much edge details and point targets become even more blurred. Experimental results confirmed that 3 × 3 window is an appropriate choice. The algorithm uses an a priori estimate mean u’(y) to replace u(y) to reduce image noise bias and Δx is replaced by N(x). The new Bayesian filter can be expressed as(3.8)n(x)=∑y∈Δ(x)p(v(x)∣u′(y))p(u′(y))u′(y)∑y∈Δ(x)p(v(x)∣u′(y))p(u′(y)).N(x) can be expressed as  N(x)=Δx∩N1(x)∩N2(x), where N1(x) and N2(x) are a prioriregional image characteristics and pixel features, respectively. The a priori regional characteristic is image region Δx, and unrelated points are removed using a region similarity algorithm. The a priori pixel feature is the set obtained by comparing the similarity of adjacent pixels [27]. A priori pixel characteristics are generally always overlooked in NLM filter processes. In fact, a priori pixel characteristics are good for excluding pixel noise [28].The sigma range between pixelx and the a priori mean u’(x) can be defined as (u’(x)I1,u’(x)I2) and the range (I1,I2) meets ξ=∫I1I2p(s)ds, where p(s) is the image probability density function. For different sigma values ξ∈{0.1,0.2,…,0.9}, the range can be calculated by pixel search [29].It is desirable to have a greater sigma weight; however, under conditional probability, the sigma range cannot be greater than the maximum upper boundary,u’(x)I2<Vmax, where Vmax is the maximum image density. It has been demonstrated that a priori pixel characteristics can have a good effect on retention of image edges, but there will be some situations in which isolated pixels are ignored. To solve this problem, the proposed algorithm uses a threshold T=Vmax/2 to separate two pixels [30]. For a priori pixels, only u’(x)<T are retained.


## 3.1. NLM Filter

Different frequency components play different roles in an image structure. Low-frequency components account for most image energy, forming basic local gradation areas, but play a small role in image content or structure. High-frequency components form the main image edges and determine its basic content or structure and are thus the most important components. Changes in high-frequency information lead to changes in the basic image content or structure, and information extracted from the image by the human eye will thus be subject to major changes. Thus, high-frequency components play the most important role in image perception by the human eye.At present, many image filters only consider adjacent pixels; some filters take into account information for neighboring pixels, such as Yaroslavsky neighborhood filters [22] and bilateral filters [23]. A nonlinear filter involves additive white noise and can effectively handle image redundancy [24].The NLM algorithm takes advantage of grayscale image redundancy and structural redundancy through a weighted average of pixel values to estimate the current pixel value. The value of each pixel is calculated using the Gaussian-weighted Euclidean distance between subblocks; a pixel as taken as the right center of the corresponding subblock. This ensures that pixels with a similar structure are assigned greater weight. For an original imagev={v(i)∣i∈I}, the expression for the image processed using the NLM filter, NL(v), is(3.1)NL(v,i)=∑j∈Iw(i,j)v(j),
where w(i,j) is the Gaussian-weighted Euclidean distance between pixels i and j, which represents the similarity of the image subblocks with i and j as centers:(3.2)w(i,j)=1Z(i)exp⁡[-‖v(Ni)d×d-v(Nj)d×d‖2h2],Z(i)=∑jexp⁡[-‖v(Ni)d×d-v(Nj)d×d‖2a2h2],h2=1σ2(1+r2)∑iN∑jN(rNi-Nj)2,
where Z(i) is a normalization factor, v(Ni) is a set of subblocks with pixel i as the center, ∥·∥2a is a Gaussian-weighted Euclidean distance function, a is Gaussian kernel standard deviation, and h is a filter parameter that controls the degree of smoothing. The weight w(i,j) meets 0≤w(i,j)≤1, if neighborhood pixels are more similar with v(Ni), the weight of center pixel is greater.


## 3.2. NLM Filter Combined with a Bayesian Method

The traditional NLM algorithm is very similar in both grayscale and structure content for smooth neighborhood areas. The algorithm yields the best results in flat areas, where better denoising effects can be obtained. At image edges and in texture-rich regions, the algorithm performs poorly because these regions have many repeat structures, the difference in grayscale content is greater, and the larger Euclidean distance makes the weights very small and reduces denoising capability, especially the ability to retain image detail [25]. To improve the edge retention capacity of the NLM algorithm, a Bayesian algorithm was added to make use of image edge information and adjust the similarity of the neighborhood structure so that the center pixel of edge contents that are similar can be given greater weight. This provides a more effective approach for protecting image detail.The Bayesian NLM filter is expressed as(3.3)n(x)=∑y∈Δ(x)p(v(x)∣u(y))p(u(y))u(y)∑y∈Δ(x)p(v(x)∣u(y))p(u(y)),
where v(x) represents noise data, u(x) represents nonnoise image data, n(x) is the average pixel precision weight for gray value u(y) at a rate of change Δx, and p(v(x)∣u(y))p(u(y)) is the similarity between v(x) and u(y). Equation (5) can then be rewritten as(3.4)w(i,j)=1Z(i)exp⁡[-‖v(Ni)d×d-v(Nj)d×d‖2a2h2-‖n(xi)d×d-n(xj)d×d‖2a2h′].


## 3.3. Improved Bayesian NLM Filter

We propose an improved Bayesian NLM filter in which a sigma range is added to the prior image operation for more effective protection of image details [26]. The first step is analysis of the probability density for pixel levels, which takes the average variance for the improved Bayesian filter. Considering the independence and integrity of an image, its conditional probability distribution can be expressed as(3.5)p(v(x)∣u(y))=∑m=1M×Mp(vm(x)∣um(y)),
where vm(x) and um(y) are the image probability densities at pixel m and vm(y) is a subset of um(y).For L-level image decomposition, the conditional probability density function is(3.6)p(vm(x)∣um(y))=vm(x)L-1F(L)(Lum(y))exp⁡(-Lvm(x)um(y)),(3.7)F(L)=Ga(L)∑lLexp⁡(um(y)l-1-vm(x)l-1).Because of multiscale features, we assume that the prior probabilityp(u(y)) is continuous and uniform, p(u(y))=1/|Δx|. The proposed algorithm uses an iterative technique, which takes the observed value v(x) as the initial value u(y). This treatment can process data directly but takes longer, and details can become fuzzy. If the frequency window is too large, the result will have too much edge details and point targets become even more blurred. Experimental results confirmed that 3 × 3 window is an appropriate choice. The algorithm uses an a priori estimate mean u’(y) to replace u(y) to reduce image noise bias and Δx is replaced by N(x). The new Bayesian filter can be expressed as(3.8)n(x)=∑y∈Δ(x)p(v(x)∣u′(y))p(u′(y))u′(y)∑y∈Δ(x)p(v(x)∣u′(y))p(u′(y)).N(x) can be expressed as  N(x)=Δx∩N1(x)∩N2(x), where N1(x) and N2(x) are a prioriregional image characteristics and pixel features, respectively. The a priori regional characteristic is image region Δx, and unrelated points are removed using a region similarity algorithm. The a priori pixel feature is the set obtained by comparing the similarity of adjacent pixels [27]. A priori pixel characteristics are generally always overlooked in NLM filter processes. In fact, a priori pixel characteristics are good for excluding pixel noise [28].The sigma range between pixelx and the a priori mean u’(x) can be defined as (u’(x)I1,u’(x)I2) and the range (I1,I2) meets ξ=∫I1I2p(s)ds, where p(s) is the image probability density function. For different sigma values ξ∈{0.1,0.2,…,0.9}, the range can be calculated by pixel search [29].It is desirable to have a greater sigma weight; however, under conditional probability, the sigma range cannot be greater than the maximum upper boundary,u’(x)I2<Vmax, where Vmax is the maximum image density. It has been demonstrated that a priori pixel characteristics can have a good effect on retention of image edges, but there will be some situations in which isolated pixels are ignored. To solve this problem, the proposed algorithm uses a threshold T=Vmax/2 to separate two pixels [30]. For a priori pixels, only u’(x)<T are retained.


## 4. Image Feature Extraction Based on the Contourlet Transform

The proposed algorithm can improve the accuracy and completeness of image feature extraction based on direct contourlet decomposition. An image is first processed by the contourlet transform to yield a multidimensional domain, with multiple-resolution decomposition coefficients for large-scale details (low-frequency signal) and finer image details (high-frequency signal). Next, the algorithm applies deeper decomposition to the large-scale approximation. The whole process can be repeated until the algorithm yields the detail required.Figure5 shows the two-layer decomposition, where I is the observed image, I’ is the image processed using the contourlet transform, LP is a Laplacian pyramid decomposition filter, DFB is a direction filter bank, LF is the low-frequency signal and HF is the high-frequency signal.Figure 5
Image feature extraction based on the contourlet transform.The decomposition coefficients for different frequencies are processed using the Bayesian-based NLM filter with a decomposition threshold. In particular, we use the wavelet threshold approach for the low-frequency part and the NLM approach for the high-frequency part. The specific steps in the algorithm are as follows.Step 1.
Decompose image I using the nonsubsampled contourlet transform.Step 2.
Apply the decomposition threshold method to the low-frequency part for noise suppression and feature extraction.Step 3.
Apply the improved Bayesian NLM filter to the high-frequency part for feature extraction.Step 4.
Reconstruct the high-frequency and low-frequency parts of the image processed using the contourlet transform.For the low-frequency part, threshold decomposition is used to remove image noise. First, the threshold valueT is set. Decomposition coefficients smaller than T are considered to be noise and thus are set to zero; coefficients greater than T are reserved. The decomposition threshold is(4.1)T=θn2k2log⁡N,θn=M(d1)0.6745,M(d1)=∑i=0N∑n∈z2d1[2n+ki],
where k is the number of layers for wavelet decomposition, θn is a function for estimating mean absolute deviation, where d1 is the high-frequency coefficient for first-layer contourlet decomposition. The high-frequency part of the first layer usually contains few signal components and comprises mainly noise. For the high-frequency part, G[Dk(x)] and G[Dk(y)] are two blocks for image I.G[Dk(p)] represents the rectangular neighborhood around p as center. The proposed algorithm uses an improved Bayesian NLM image filter and a Euclidean distance to represent similarity between high-frequency image blocks. The similarity is represented by w(x,y), which is defined as(4.2)w(x,y)=1Z(x)e-(‖G[Dk(x)]-G[Dk(y)]‖2/h2),Z(x)=∑x,y∈Ie-(‖G[Dk(x)]-G[Dk(y)]‖2/h2),
where h is a constant used to control the exponential decay rate. Compared with the original NLM filter, the contourlet transform decomposes the image at different resolutions and the proposed approach uses different algorithm to process an image: threshold decomposition is used for the low-frequency part, and an improved Bayesian NLM filter is used for the high-frequency part. The NLM filter involves time-consuming calculations. If an image is decomposed by the contourlet transform for k levels, the improved algorithm only has to process 1/2k of the original size. This not only reduces the computational complexity but also greatly improves the accuracy of feature extraction.


## 5. Experimental Results

Spatial image can directly influence the success of spatial rendezvous and docking. We must ensure that the most accurate image feature can be extracted from spatial images. The proposed algorithm can also be applied to general image processing. To verify the performance of the algorithm, we carried out experiments both on Spatial images and standard images.


### 5.1. Performance Evaluation Based on Spatial Images

Tests were carried out on image I and image II for an image size of  512 × 512. Multidimensional contourlet and curvelet decomposition algorithms were used for comparison. The same parameters were used for all algorithms.Image feature extraction results were evaluated according to subjective and objective standards. Figure6 shows the processing results for image I, and Figure 7 shows the results for image II. It is evident that the proposed algorithm yields a better subjective visual effect compared with the curvelet and contourlet algorithms.Process results for image I.
(a)
Original image(b)
Curvelet transform(c)
Contourlet transform(d)
Contourlet + NL-meansProcess results for image II.
(a)
Original image(b)
Curvelet transform(c)
Contourlet transform(d)
Contourlet + NL-meansThe results show that feature extraction with the curvelet transform leads to confusion for some background information, and the contourlet transform yields a blurry image. By contrast, the proposed algorithm effectively suppresses noise and displays the main features of the image. Finer image details are shown in Figure8.Process results for image II.
(a)
Original(b)
Noisy(c)
Curvelet(d)
Contourlet(e)
Contourlet + NL-means


### 5.2. Performance Evaluation Based on Standard Images

Figures9 and 10 show the processing results for image as III, IV. Image as III and IV are standard images for an image size of 512 × 512. It is shown that our algorithm also provides better performance than that of Curvelet transform and Contourlet transform.Process results for image III.
(a)
Original image(b)
Curvelet transform(c)
Contourlet transform(d)
Contourlet + NL-meansProcess results for image IV.
(a)
Original image(b)
Curvelet transform(c)
Contourlet transform(d)
Contourlet + NL-meansResults are objectively evaluated using the peak signal-to-noise ratio (PSNR) and structural similarity (SSIM), defined as follows:(5.1)PSNR=10log10(L-1)2∑i=1M∑j=1N[R(i,j)-F(i,j)]2,SSIM(x,y)=(2uxuy+C1)(2σxy+C2)(ux2+uy2+C1)(σx2+σy2+C2),MSSIM(X,Y)=1W∑r=1WSSIM(Xr,Yr),
where ux and uy are the mean and σx and σy are the standard deviation for the original and processed images, respectively, σxy is the covariance for the original and processed images, and C1 and C2 are constants. MSSIM is mean SSIM, and W is the number of image subblocks. For greater PSNR and MSSIM (0≤MSSIM≤1), the processed image is closer to the original.Table1 compares the feature extraction performance for images I and II for several algorithms and different noise levels. The PSNR results show that the contourlet transform is superior to the curvelet transform for image I by almost 0.5 dB. Both the contourlet and curvelet transforms provide good edge detection. The proposed  CT+NLM  algorithm showed even better performance (~1.1 dB) compared with the contourlet transform but retained the good edge detection of the latter method (Table 1).The proposed algorithm uses an NLM filter for adaptive image expression. The MSSIM results show that the proposed algorithm yields the best performance for Gaussian, Poisson, Salt and Pepper and Speckle noise (Table 2).Table 1
PSNR results.
NoiseσPSNR (dB)NoisyOriginalCurveletContourletCT+NLM10Image I31.2232.9333.0133.5834.62Image II32.9833.5333.8934.2135.4020Image I31.1631.8932.7533.0134.19Image II32.4933.1433.1933.8234.8930Image I30.7831.5632.1332.7833.76Image II31.9832.7532.8233.2134.1740Image I30.5631.0131.7432.5333.61Image II31.2032.1132.4232.9133.9850Image I29.8830.5431.2931.7432.83Image II30.7731.5631.9832.1733.28Table 2
MSSIM results.
NoiseMSSIMNoisyOriginalCurveletContourletCT+NLMGaussianImage I0.29790.31280.35240.42770.5185Image II0.30140.32560.36880.43160.5220PoissonImage I0.57440.62310.68270.71730.8219Image II0.58910.63490.70280.73390.8551Salt and pepperImage I0.29480.32390.37400.41570.5255Image II0.30010.33980.38290.42540.5345SpeckleImage I0.29530.33210.39740.43090.5312Image II0.30210.34760.40120.44520.5422Our algorithm uses a nonsubsampled key point filter for whichJ+1 redundancy is the most efficient. In pyramid decomposition, a lesser extent of image loss can be considered as an effective means to reduce redundancy. The proposed algorithm, which uses a nonsubsampled pyramid filter and a directional filter, leads to some image loss in reducing redundancy. Search windows of 16×16, 32×32, and 64×64 were applied to images I and II. The size of the search window can affect the computational complexity of the NLM filter.Comparison of the experimental results for different window sizes reveals that the proposed algorithm delivers better noise suppression and feature extraction than the other algorithms. It provides a maximum PSNR value and a minimum MSSIM value for all windows (Table3).Table 3
Experimental results for different search windows.
SchemeSearch window16×16Search window32×32Search window64×64Image IImage IIImage IImage IIImage IImage IIPSNRSSIMPSNRSSIMPSNRSSIMPSNRSSIMPSNRSSIMPSNRSSIMOriginal30.530.65831.450.67930.090.66931.210.68529.890.67531.090.703Curvelet30.780.69731.660.72230.120.70131.390.73729.910.68831.150.749Contourlet31.230.75931.980.78430.350.76231.650.79130.170.74431.230.806CT+NLM32.380.78332.890.80231.370.78832.680.81131.230.79632.220.835


## 5.1. Performance Evaluation Based on Spatial Images

Tests were carried out on image I and image II for an image size of  512 × 512. Multidimensional contourlet and curvelet decomposition algorithms were used for comparison. The same parameters were used for all algorithms.Image feature extraction results were evaluated according to subjective and objective standards. Figure6 shows the processing results for image I, and Figure 7 shows the results for image II. It is evident that the proposed algorithm yields a better subjective visual effect compared with the curvelet and contourlet algorithms.Process results for image I.
(a)
Original image(b)
Curvelet transform(c)
Contourlet transform(d)
Contourlet + NL-meansProcess results for image II.
(a)
Original image(b)
Curvelet transform(c)
Contourlet transform(d)
Contourlet + NL-meansThe results show that feature extraction with the curvelet transform leads to confusion for some background information, and the contourlet transform yields a blurry image. By contrast, the proposed algorithm effectively suppresses noise and displays the main features of the image. Finer image details are shown in Figure8.Process results for image II.
(a)
Original(b)
Noisy(c)
Curvelet(d)
Contourlet(e)
Contourlet + NL-means


## 5.2. Performance Evaluation Based on Standard Images

Figures9 and 10 show the processing results for image as III, IV. Image as III and IV are standard images for an image size of 512 × 512. It is shown that our algorithm also provides better performance than that of Curvelet transform and Contourlet transform.Process results for image III.
(a)
Original image(b)
Curvelet transform(c)
Contourlet transform(d)
Contourlet + NL-meansProcess results for image IV.
(a)
Original image(b)
Curvelet transform(c)
Contourlet transform(d)
Contourlet + NL-meansResults are objectively evaluated using the peak signal-to-noise ratio (PSNR) and structural similarity (SSIM), defined as follows:(5.1)PSNR=10log10(L-1)2∑i=1M∑j=1N[R(i,j)-F(i,j)]2,SSIM(x,y)=(2uxuy+C1)(2σxy+C2)(ux2+uy2+C1)(σx2+σy2+C2),MSSIM(X,Y)=1W∑r=1WSSIM(Xr,Yr),
where ux and uy are the mean and σx and σy are the standard deviation for the original and processed images, respectively, σxy is the covariance for the original and processed images, and C1 and C2 are constants. MSSIM is mean SSIM, and W is the number of image subblocks. For greater PSNR and MSSIM (0≤MSSIM≤1), the processed image is closer to the original.Table1 compares the feature extraction performance for images I and II for several algorithms and different noise levels. The PSNR results show that the contourlet transform is superior to the curvelet transform for image I by almost 0.5 dB. Both the contourlet and curvelet transforms provide good edge detection. The proposed  CT+NLM  algorithm showed even better performance (~1.1 dB) compared with the contourlet transform but retained the good edge detection of the latter method (Table 1).The proposed algorithm uses an NLM filter for adaptive image expression. The MSSIM results show that the proposed algorithm yields the best performance for Gaussian, Poisson, Salt and Pepper and Speckle noise (Table 2).Table 1
PSNR results.
NoiseσPSNR (dB)NoisyOriginalCurveletContourletCT+NLM10Image I31.2232.9333.0133.5834.62Image II32.9833.5333.8934.2135.4020Image I31.1631.8932.7533.0134.19Image II32.4933.1433.1933.8234.8930Image I30.7831.5632.1332.7833.76Image II31.9832.7532.8233.2134.1740Image I30.5631.0131.7432.5333.61Image II31.2032.1132.4232.9133.9850Image I29.8830.5431.2931.7432.83Image II30.7731.5631.9832.1733.28Table 2
MSSIM results.
NoiseMSSIMNoisyOriginalCurveletContourletCT+NLMGaussianImage I0.29790.31280.35240.42770.5185Image II0.30140.32560.36880.43160.5220PoissonImage I0.57440.62310.68270.71730.8219Image II0.58910.63490.70280.73390.8551Salt and pepperImage I0.29480.32390.37400.41570.5255Image II0.30010.33980.38290.42540.5345SpeckleImage I0.29530.33210.39740.43090.5312Image II0.30210.34760.40120.44520.5422Our algorithm uses a nonsubsampled key point filter for whichJ+1 redundancy is the most efficient. In pyramid decomposition, a lesser extent of image loss can be considered as an effective means to reduce redundancy. The proposed algorithm, which uses a nonsubsampled pyramid filter and a directional filter, leads to some image loss in reducing redundancy. Search windows of 16×16, 32×32, and 64×64 were applied to images I and II. The size of the search window can affect the computational complexity of the NLM filter.Comparison of the experimental results for different window sizes reveals that the proposed algorithm delivers better noise suppression and feature extraction than the other algorithms. It provides a maximum PSNR value and a minimum MSSIM value for all windows (Table3).Table 3
Experimental results for different search windows.
SchemeSearch window16×16Search window32×32Search window64×64Image IImage IIImage IImage IIImage IImage IIPSNRSSIMPSNRSSIMPSNRSSIMPSNRSSIMPSNRSSIMPSNRSSIMOriginal30.530.65831.450.67930.090.66931.210.68529.890.67531.090.703Curvelet30.780.69731.660.72230.120.70131.390.73729.910.68831.150.749Contourlet31.230.75931.980.78430.350.76231.650.79130.170.74431.230.806CT+NLM32.380.78332.890.80231.370.78832.680.81131.230.79632.220.835


## 6. Conclusions

Focusing on the actual needs of spatial images analysis, an improved contourlet transform, consisting of a nonsubsampled pyramid transform and nonsubsampled directional filter banks, was used to reduce the filter design problem of spatial images. The improved contourlet transform uses a mapping approach to solve the 2D filter design problem. The algorithm uses a Bayesian NLM filter for high-frequency information to suppress noise and improve the accuracy of image feature extraction. Experimental results confirm that the NLM filter can effectively retain structural information and reduce the residual structure. In the NSCT domain, the proposed algorithm showed better denoising and enhancement effects compared with the contourlet transform. Moreover, in comparison with NSWT, the algorithm is a more mature and sophisticated image-processing method.



---

*Source: 467412-2012-06-04.xml*