# A Correlation Based Strategy for the Acceleration of Nonlocal Means Filtering Algorithm

Authors: Junfeng Zhang; Jiasong Wu; Jean-Louis Coatrieux; Limin Luo; Yang Chen
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Although the nonlocal means (NLM) algorithm takes a significant step forward in image filtering field, it suffers from a high computational complexity. To deal with this drawback, this paper proposes an acceleration strategy based on a correlation operation. Instead of per-pixel processing, this approach performs a simultaneous calculation of all the image pixels with the help of correlation operators. Complexity analysis and experimental results are reported and show the advantage of the proposed algorithm in terms of computation and time cost.


---

## Body

## 1. Introduction

Digital image denoising has been a fundamental and challenging issue for several decades [1, 2]. Many contributions have been devoted to recover the image degraded by Gaussian noise. Much attention has been paid to the Partial Differential Equation (PDE) approaches among which the total variation (TV) and the Perona and Malik (PM) model are well known [3–5]. Other methods rely on the image transform from the spatial domain to another domain (such as the Fourier domain, the wavelet domain, and the DCT domain). After adjusting the transform coefficients, the image is restituted by applying the inverse transform [6, 7].The algorithm of nonlocal means (NLM) filtering was proposed by Buades et al. [8]. They suggested that a denoised pixel is equivalent to the weighted average of its neighboring pixels, with the weights calculated by the normalized Gaussian weighted Euclidean distance between the blocks centred at those pixels. The NLM algorithm has demonstrated better performance than other main-stream filter methods, such as bilateral filter and TV model in both visual performance and objective measure [9]. Unfortunately, the computational cost of computing the weights is too expensive in many applications.Some solutions have been proposed to alleviate this high computation burden for weights’ calculation. Liu et al. proposed an approximation to the similarity of neighborhood windows by employing an efficient Summed Square Image (SSI) further combined with Fast Fourier Transform (FFT) [10, 11]. In [12, 13], the NLM algorithm is accelerated by eliminating some computation of weights through a preclassification step based on a hard threshold of local block measures (average intensities, gradients, and first- and second-order moments). Also, some Probabilistic Early Termination (PET) schemes, such as cluster tree, Singular Value Decomposition (SVD), or dictionaries for image blocks and image edges, were also employed to speed up the weights’ calculation [14–17].In this paper, compared with the traditional NLM algorithm, in which the calculations of weights are by the pixel by pixel way, our proposed fast strategy was performed on the whole image. Specifically, correlation operators are applied to compute the differential image and lead to a straightforward shortcut to achieve all the weights. By doing this, a lot of redundancy computation can be successfully avoided. Thereby, the computation speed of NLM is improved.This paper is organized as follows. In Section2, we briefly review the NLM algorithm. The repetitive computation causing the original per-pixel NLM slowness is analyzed and our proposed fast NLM algorithm is described in Section 3. In Section 4, some experiments are conducted and they show that a significant improvement over competitive approaches is brought by our method. Section 5 concludes this paper and some relative discussions are given.


## 2. The Non-Local Means Algorithm

### 2.1. The Principle of NLM

The image denoising model for an imageF
∈
R
M
×
M, where R is the real value domain, can be formulated by(1)
F
i
=
U
i
+
N
i
,
1
≤
i
≤
M
,where F
(
i
) denotes the noisy image value, U
(
i
) represents the original noise-free image value, and N
(
i
) is the noise value at the pixel i. The idea behind the NLM filter is to consider that the denoised pixel value of i is equivalent to the weighted mean of all pixels’ values of the noisy image (indexed in image F) [1, 8]. However, considering the high computational cost, it was suggested in [1, 8] that we can estimate the denoised pixel value F
(
i
) using the pixels within a larger search neighborhood S
i
(
2
L
+
1
)
×
(
2
L
+
1
)
⊆
R
M
×
M centred at the pixeli, where L denotes the radius of S
i, whose pixel values are scanned column by column and then concatenated to a vector:(2)
F
S
i
=
F
s
i
,
1
,
F
s
i
,
2
,
…
,
F
s
i
,
q
,
…
,
F
s
i
,
K
T
,where K
=
(
2
L
+
1
)
×
(
2
L
+
1
) is the number of pixels within S
i. The superscript “T” denotes the transpose and F
(
s
i
,
q
), 1
≤
q
≤
K, is the value of the qth pixel in the search neighborhood S
i. In this situation, the NLM algorithm can be expressed by(3)
NLM
F
i
=
∑
s
i
,
q
∈
S
i
W
i
,
s
i
,
q
F
s
i
,
q
,where N
L
M
(
F
(
i
)
) is the filtered image value of the pixel i and W
(
i
,
s
i
,
q
) represents the weight between F
(
i
) and F
(
s
i
,
q
). The weight is acquired by(4)
W
i
,
s
i
,
q
=
1
Z
i
W
~
i
,
s
i
,
q
,where Z
(
i
) is a normalization factor:(5)
Z
i
=
∑
s
i
,
q
∈
S
i
W
~
i
,
s
i
,
q
,and W
~
i
,
s
i
,
q denotes the exponential Gaussian weighted Euclidean distance:(6)
W
~
i
,
s
i
,
q
=
e
-
D
a
i
,
s
i
,
q
/
h
2
.In (6), h acts as a filtering parameter and D
a
(
i
,
s
i
,
q
) is the Gaussian Euclidean distance between blocks B
i
∈
R
(
2
l
+
1
)
×
(
2
l
+
1
) and B
s
i
,
q
∈
R
(
2
l
+
1
)
×
(
2
l
+
1
) centred at the pixels i and s
i
,
q, respectively, and is given by (7)
D
a
i
,
s
i
,
q
=
f
B
i
-
f
B
s
i
,
q
2
,
a
2
.In (7),(8)
f
B
i
=
F
b
i
,
1
,
F
b
i
,
2
,
…
,
F
b
i
,
k
T
,obtained by scanning the pixel values of the block B
i column by column and then concatenating them and k is the number of pixels within the square block B
i whose radius is set as l. That is to say, k
=
(
2
l
+
1
)
×
(
2
l
+
1
). Similarly,(9)
f
B
s
i
,
q
=
F
b
s
i
,
q
,
1
,
F
b
s
i
,
q
,
2
,
…
,
F
b
s
i
,
q
,
k
Tis achieved by scanning the pixel values of the block B
s
i
,
qcolumn by column and then concatenating them. Using (8) and (9), (7) can be deduced step by step as follows:(10)
D
a
i
,
s
i
,
q
=
∑
m
=
1
k
G
m
F
¯
s
i
,
q
m
,
(11)
F
¯
s
i
,
q
m
=
F
~
s
i
,
q
2
m
,
(12)
F
~
s
i
,
q
m
=
F
b
i
,
m
-
F
b
s
i
,
q
,
m
,
1
≤
m
≤
k
,where G
(
m
) is the mth component of the discrete Gaussian kernel G
∈
R
(
2
l
+
1
)
×
(
2
l
+
1
), whose standard deviation is a. Apparently it has the same size as the block B.


### 2.2. Computation Cost Using Per-Pixel Algorithm

From (10)–(12), we can clearly see the cause of the NLM computational complexity. For every pixel i, the NLM algorithm must compute the dissimilarity between f
B
i and all f
B
s
i
,
q in the search neighborhood S
i. For instance, to filter a noisy image F of size M
×
M, it needs 2
×
k
×
K
×
M multiplications and (2
k
-
1
)
×
K
×
M additions. By cutting down the search neighborhood S
i into 15 × 15, with the block B
i set as 3 × 3, for a general 512 × 512 size image F, it needs 2 × 92  × 152  × 5122 multiplications and 17 × 152  × 5122 additions. Hence, a fast NLM algorithm is required. This motivated us to propose a correlation based accelerated nonlocal denoising algorithm.


## 2.1. The Principle of NLM

The image denoising model for an imageF
∈
R
M
×
M, where R is the real value domain, can be formulated by(1)
F
i
=
U
i
+
N
i
,
1
≤
i
≤
M
,where F
(
i
) denotes the noisy image value, U
(
i
) represents the original noise-free image value, and N
(
i
) is the noise value at the pixel i. The idea behind the NLM filter is to consider that the denoised pixel value of i is equivalent to the weighted mean of all pixels’ values of the noisy image (indexed in image F) [1, 8]. However, considering the high computational cost, it was suggested in [1, 8] that we can estimate the denoised pixel value F
(
i
) using the pixels within a larger search neighborhood S
i
(
2
L
+
1
)
×
(
2
L
+
1
)
⊆
R
M
×
M centred at the pixeli, where L denotes the radius of S
i, whose pixel values are scanned column by column and then concatenated to a vector:(2)
F
S
i
=
F
s
i
,
1
,
F
s
i
,
2
,
…
,
F
s
i
,
q
,
…
,
F
s
i
,
K
T
,where K
=
(
2
L
+
1
)
×
(
2
L
+
1
) is the number of pixels within S
i. The superscript “T” denotes the transpose and F
(
s
i
,
q
), 1
≤
q
≤
K, is the value of the qth pixel in the search neighborhood S
i. In this situation, the NLM algorithm can be expressed by(3)
NLM
F
i
=
∑
s
i
,
q
∈
S
i
W
i
,
s
i
,
q
F
s
i
,
q
,where N
L
M
(
F
(
i
)
) is the filtered image value of the pixel i and W
(
i
,
s
i
,
q
) represents the weight between F
(
i
) and F
(
s
i
,
q
). The weight is acquired by(4)
W
i
,
s
i
,
q
=
1
Z
i
W
~
i
,
s
i
,
q
,where Z
(
i
) is a normalization factor:(5)
Z
i
=
∑
s
i
,
q
∈
S
i
W
~
i
,
s
i
,
q
,and W
~
i
,
s
i
,
q denotes the exponential Gaussian weighted Euclidean distance:(6)
W
~
i
,
s
i
,
q
=
e
-
D
a
i
,
s
i
,
q
/
h
2
.In (6), h acts as a filtering parameter and D
a
(
i
,
s
i
,
q
) is the Gaussian Euclidean distance between blocks B
i
∈
R
(
2
l
+
1
)
×
(
2
l
+
1
) and B
s
i
,
q
∈
R
(
2
l
+
1
)
×
(
2
l
+
1
) centred at the pixels i and s
i
,
q, respectively, and is given by (7)
D
a
i
,
s
i
,
q
=
f
B
i
-
f
B
s
i
,
q
2
,
a
2
.In (7),(8)
f
B
i
=
F
b
i
,
1
,
F
b
i
,
2
,
…
,
F
b
i
,
k
T
,obtained by scanning the pixel values of the block B
i column by column and then concatenating them and k is the number of pixels within the square block B
i whose radius is set as l. That is to say, k
=
(
2
l
+
1
)
×
(
2
l
+
1
). Similarly,(9)
f
B
s
i
,
q
=
F
b
s
i
,
q
,
1
,
F
b
s
i
,
q
,
2
,
…
,
F
b
s
i
,
q
,
k
Tis achieved by scanning the pixel values of the block B
s
i
,
qcolumn by column and then concatenating them. Using (8) and (9), (7) can be deduced step by step as follows:(10)
D
a
i
,
s
i
,
q
=
∑
m
=
1
k
G
m
F
¯
s
i
,
q
m
,
(11)
F
¯
s
i
,
q
m
=
F
~
s
i
,
q
2
m
,
(12)
F
~
s
i
,
q
m
=
F
b
i
,
m
-
F
b
s
i
,
q
,
m
,
1
≤
m
≤
k
,where G
(
m
) is the mth component of the discrete Gaussian kernel G
∈
R
(
2
l
+
1
)
×
(
2
l
+
1
), whose standard deviation is a. Apparently it has the same size as the block B.


## 2.2. Computation Cost Using Per-Pixel Algorithm

From (10)–(12), we can clearly see the cause of the NLM computational complexity. For every pixel i, the NLM algorithm must compute the dissimilarity between f
B
i and all f
B
s
i
,
q in the search neighborhood S
i. For instance, to filter a noisy image F of size M
×
M, it needs 2
×
k
×
K
×
M multiplications and (2
k
-
1
)
×
K
×
M additions. By cutting down the search neighborhood S
i into 15 × 15, with the block B
i set as 3 × 3, for a general 512 × 512 size image F, it needs 2 × 92  × 152  × 5122 multiplications and 17 × 152  × 5122 additions. Hence, a fast NLM algorithm is required. This motivated us to propose a correlation based accelerated nonlocal denoising algorithm.


## 3. The Fast Nonlocal Means Algorithm

### 3.1. Repetitive Computation of the Per-Pixel Algorithm

From the above observations, we can see that, for one specific central pixeli, the Gaussian weighted Euclidean distance d
a
(
i
,
s
i
,
q
) between the block B
i centred at i and block B
s
i
,
q centred at s
i
,
q in the neighboring needs to be computed. However, for any pixel j adjoining the pixel i, the Gaussian weighted Euclidean distance d
a
(
j
,
s
j
,
q
) between the block B
j centred at j and block B
s
j
,
q centred at s
j
,
q in the neighboring also needs to be computed. Because the two pixels are adjacent, B
i and B
j have lots of overlap and are the same as B
s
i
,
q and B
s
j
,
q. Therefore, by theper-pixel way to implement the NLM algorithm, the blocks based Gaussian weighted Euclidean calculation is highly repetitive.More precisely, for a blockB with size k
=
(
2
×
l
+
1
)
×
(
2
×
l
+
1
), the repetitive calculated square term is 2
×
l
2. The repetitive rate is (
(
2
×
l
+
1
)
×
(
2
×
l
)
)
/
2
×
l
+
1
2
=
(
2
×
l
)
/
(
2
×
l
+
1
). To get an intuitive observation, a diagram is depicted in Figure 1 in which we set the size of block B
i as k =3 × 3. It can be observed that there are 2/3 repetitive operators of intensities differences when calculating the weights for the adjacent centered values of pixels i and j.Figure 1
Two adjacent pixelsi and j with their blocks B
i and B
j, and their search neighborhoods S
i and S
j. The pixels within the 11 × 11 blue block diagram represent the search neighborhood S
i of the pixel i, the 5 × 5 blue block diagram is the block B
i of the pixel i, and the centered blue block diagram is the pixel i. The red block diagrams (from small to big) correspond to the pixel j, block B
j, and search neighbourhood S
j, respectively.


### 3.2. Our Proposed Fast Algorithm

To avoid repetitive computations of the pixel differenceF
~
s
i
,
q
(
m
)
(
1
≤
i
≤
M, 1
≤
m
≤
k
) in (12) by the traditional pixel by pixel way, the proposed correlation based accelerated fast NLM algorithm deals with the image F as a whole. In this case, it is necessary to deal with (1)–(9) in matrix form. Mathematically, it is trivial to extend them to the nonscalar matrix case. Then, the variables in (1)–(9) are restricted to all the M pixels in the noisy image F rather than some specific pixel i.To deal with the operation with the whole imageF, our algorithm can be accordingly divided into 8 steps being consistent with the opposite sequence of (1)–(9).Step 1.
Compute the “differential image”F
~
s
q, which is matrix case consistent with F
~
s
i
,
q
(
m
) in (12).This step can be efficiently implemented via a correlation operator:(13)
F
~
s
q
=
P
s
q
∘
F
,where “∘” denotes the correlation operation and F
~
s
q
∈
R
M
×
M denotes the qth differential image of the original image F. P
s
q
∈
R
(
2
L
+
1
)
×
(
2
L
+
1
) is a 2D correlation kernel, which contains only two nonzero elements; that is, the center pixel value is 1 and the qth pixel value is -
1. Note that, for every different q values, there is a different correlation kernel P
s
q since q determines the position of the pixel value -
1. Figure 2 gives an illustration of the q = 5th differential image of the whole image within an 11 × 11 search neighborhood.Figure 2
The introduced correlation kernel matrix and its correspondingv.By means of (13), the achieved differential image contains all the available information in (12). In fact, (13) can be straightforwardly implemented by a translation operation:(14)
F
~
s
q
=
F
v
q
-
F
,where F
v
q represents the image shifted by a coordinate vector v
q, 1
≤
q
≤
K. Here, the qth vector v
q points out the coordinate (x
q, y
q) from the origin (0
,
0) as shown in Figure 2. Note that the correlation kernel P
s
q and v
q, 1
≤
q
≤
K, are one-to-one mapping; that is, when (15)
q
=
c
2
L
+
1
+
r
+
1
,
0
≤
c
≤
2
L
,
0
≤
r
≤
2
L
,where (c
,
r) is the coordinate of “-
1” in P
s
q shown in Figure 2, then, its corresponding v
q is given by(16)
v
q
=
x
q
,
y
q
=
r
-
L
,
c
-
L
.Step 2.
Compute the “square image”F
¯
s
q, which is matrix case consistent with F
¯
s
i
,
q
(
m
) in (11).This step can be easily conducted by an element-wise square operation as follows:(17)
F
¯
s
q
=
F
~
s
q
•
F
~
s
q
,where “•” acts as an element-wise product operator.Step 3.
Compute the “Gaussian weighted Euclidean distance image”D
a
s
q, which is matrix case consistent with D
a
i
,
s
i
,
q in (10).This step can also be implemented by a correlation as follows:(18)
D
a
s
q
=
G
∘
F
¯
s
q
.
D
a
s
q
∈
R
M
×
M is the Gaussian weighted Euclidian distance corresponding to the qth pixel in the search neighborhood S for the whole image F. Equation (18) can be efficiently implemented by FFT as follows:(19)
D
a
s
q
=
I
F
F
T
F
F
T
G
∗
•
F
F
T
F
¯
s
q
.Step 4.
Compute the exponential Gaussian weighted Euclidian distance imageW
~
s
q, which is matrix case consistent with W
~
i
,
s
i
,
q in (6).Step 5.
Then, repeat Steps1
–
4 and compute all the W
~
s
q, 1
≤
s
q
≤
K.Step 6.
Compute the normalization factor imageZ, which is matrix case consistent with Z
(
i
) in (5).Step 7.
Compute the weightsW
(
s
q
), which is matrix case consistent with W
(
i
,
s
i
,
q
) in (4). This step can be implemented by an element-wise division operation.Step 8.
Compute the processedN
L
M
(
F
) in (3).


### 3.3. Symmetry Property

Note that the computational cost of Step5 can be further reduced by exploiting the symmetry property of the v
q; that is, (20)
v
q
=
-
v
K
-
q
+
1
,where v
(
K
-
q
+
1
) denotes the (
K
-
q
+
1
)th translation vector. The symmetry property is shown in Figure 3.Figure 3
The schematic representation of the symmetry property for weight calculation.Thus, when we calculate theqth weights W
~
s
q for the image, the (
K
-
q
+
1
)th weights W
~
s
(
K
-
q
+
1
) of the image are also meanwhile obtained as a secondary product, which saves the computational cost by one half.Specifically, when the kernelP
s
q is employed to yield the qth differential image F
~
s
q according to (14), the (
K
-
q
+
1
)th differential image F
~
s
(
K
-
q
+
1
) is achieved by means of (20) as follows: (21)
F
~
s
K
-
q
+
1
=
F
v
K
-
q
+
1
-
F
=
F
-
v
q
-
F
=
-
F
v
q
-
F
.Next, whenF
¯
s
q is calculated by (17), F
¯
s
(
K
-
q
+
1
) is also meanwhile achieved by(22)
F
¯
s
K
-
q
+
1
=
F
~
s
K
-
q
+
1
•
F
~
s
K
-
q
+
1
.In Steps3
-
4, only once same operation is implemented, W
~
s
q and W
~
s
(
K
-
q
+
1
) are simultaneously obtained.Intuitively, this can be observed in Figure3. By means of v
s
13, for each pixel i within the image F, the 13th differential image F
~
i
,
s
13
m (m
=
1
,
…
,
9) is available which is corresponding to difference between the red and green blue rectangle. At the same time, the negative F
~
i
,
s
109
m (m
=
1
,
…
,
9) is also available which is corresponding to the difference between the red and blue rectangle. Following Steps 3
-
4, only once operation is implemented, W
~
i
,
s
13 and W
~
i
,
s
109 are simultaneously obtained.Therefore, the necessary translation operations in terms ofv
q can be reduced from (K
-
1) to (
K
-
1
)
/
2 according to the symmetry property (see (20)). Thus, the computational cost will be reduced by one half.


### 3.4. Computational Complexity

Note also that the differences between the proposed fast algorithm and the traditional per-pixel NLM algorithm lie in Steps1
–
3. Steps 4 to 8 in the proposed algorithm include the same operations as those in the original NLM algorithm. So the computational complexity in these steps is not taken into account in the following subsection of complexity analysis.The computational complexity of the proposed algorithm is analyzed in this section. We assume that FFT and Inverse Fourier transform (IFFT) need the same computational complexity. In general, the computational complexity of two-dimensional FFT for anM
×
M image needs approximately (
5
/
9
)
M
log
2
⁡
M multiplications and (
4
/
3
)
M
log
2
⁡
M additions [18].Step 1.
ComputingF
~
s
q in (14) for all 1
≤
s
q
≤
(
K
-
1
)
/
2 requires (
K
-
1
)
/
2 shifts; therefore, M(K
-
1
)
/
2 additions are necessary.Step 2.
ComputingF
¯
s
q in (17) for all 1
≤
q
≤
(
K
-
1
)
/
2 needs (
K
-
1
)
M
/
2 multiplications.Step 3.
ComputingD
a
s
q in (19) for all 1
≤
s
q
≤
(
K
-
1
)
/
2 involves (
(
K
-
1
)
/
2
+
1
) FFTs, (
K
-
1
)
/
2 IFFTs, and (
K
-
1
)
M
/
2 multiplications for element-wise multiplications. Therefore, it needs (
5
/
9
)
M
K
log
2
⁡
M
+
(
K
-
1
)
M
/
2 multiplications and (
(
4
/
3
)
log
2
⁡
M
)
M
K additions.
The computational cost ofD
a
s
q using the proposed algorithm, the original algorithm [8], and the FFT based algorithm [11] are, respectively, given in Table 1. Both the proposed algorithm and the algorithm of [11] are superior to the original NLM algorithm in terms of the computational complexity. According to Table 1, if (23) holds, the complexity relations in (24) and (25) will also hold:(23)
M
≤
K
3
,
(24)
5
9
K
log
2
⁡
M
+
K
-
1
M
≤
15
9
K
log
2
⁡
K
+
K
+
1
M
,
(25)
4
3
K
log
2
⁡
M
+
1
2
K
-
1
M
≤
4
K
log
2
⁡
K
+
2
K
+
4
M
-
1
.Here, to satisfy (23), the minimum pixel number in the search neighborhood K
m
i
n will be M
3 and that can be used to deduce the minimal radius L
m
i
n of search region. Table 2 lists the K
m
i
n and L
m
i
n values according to different image sizes. These requirements on minimal values can be easily met in most practical cases; thus, the proposed method is superior to that of [11].Table 1
The computational complexity of different methods.
Methods
Multiplication operation number
Addition operation number
Proposed method
5
9
K
l
o
g
2
M
+
K
-
1
M
4
3
K
l
o
g
2
M
+
1
2
K
-
1
M
Method in [8]
2
k
K
M
2
k
-
1
K
M
Method in [11]
15
9
K
l
o
g
2
K
+
K
+
1
M
4
K
l
o
g
2
K
+
2
K
+
4
M
-
1Table 2
TheK
m
i
n and L
m
i
n values for different image sizes.
Image sizeM
×
M
M
K
m
i
n
L
m
i
n
2
56
×
256
65536
40.32
3
512
×
512
262144
64
4
1024
×
1024
1048576
101.59
5


## 3.1. Repetitive Computation of the Per-Pixel Algorithm

From the above observations, we can see that, for one specific central pixeli, the Gaussian weighted Euclidean distance d
a
(
i
,
s
i
,
q
) between the block B
i centred at i and block B
s
i
,
q centred at s
i
,
q in the neighboring needs to be computed. However, for any pixel j adjoining the pixel i, the Gaussian weighted Euclidean distance d
a
(
j
,
s
j
,
q
) between the block B
j centred at j and block B
s
j
,
q centred at s
j
,
q in the neighboring also needs to be computed. Because the two pixels are adjacent, B
i and B
j have lots of overlap and are the same as B
s
i
,
q and B
s
j
,
q. Therefore, by theper-pixel way to implement the NLM algorithm, the blocks based Gaussian weighted Euclidean calculation is highly repetitive.More precisely, for a blockB with size k
=
(
2
×
l
+
1
)
×
(
2
×
l
+
1
), the repetitive calculated square term is 2
×
l
2. The repetitive rate is (
(
2
×
l
+
1
)
×
(
2
×
l
)
)
/
2
×
l
+
1
2
=
(
2
×
l
)
/
(
2
×
l
+
1
). To get an intuitive observation, a diagram is depicted in Figure 1 in which we set the size of block B
i as k =3 × 3. It can be observed that there are 2/3 repetitive operators of intensities differences when calculating the weights for the adjacent centered values of pixels i and j.Figure 1
Two adjacent pixelsi and j with their blocks B
i and B
j, and their search neighborhoods S
i and S
j. The pixels within the 11 × 11 blue block diagram represent the search neighborhood S
i of the pixel i, the 5 × 5 blue block diagram is the block B
i of the pixel i, and the centered blue block diagram is the pixel i. The red block diagrams (from small to big) correspond to the pixel j, block B
j, and search neighbourhood S
j, respectively.


## 3.2. Our Proposed Fast Algorithm

To avoid repetitive computations of the pixel differenceF
~
s
i
,
q
(
m
)
(
1
≤
i
≤
M, 1
≤
m
≤
k
) in (12) by the traditional pixel by pixel way, the proposed correlation based accelerated fast NLM algorithm deals with the image F as a whole. In this case, it is necessary to deal with (1)–(9) in matrix form. Mathematically, it is trivial to extend them to the nonscalar matrix case. Then, the variables in (1)–(9) are restricted to all the M pixels in the noisy image F rather than some specific pixel i.To deal with the operation with the whole imageF, our algorithm can be accordingly divided into 8 steps being consistent with the opposite sequence of (1)–(9).Step 1.
Compute the “differential image”F
~
s
q, which is matrix case consistent with F
~
s
i
,
q
(
m
) in (12).This step can be efficiently implemented via a correlation operator:(13)
F
~
s
q
=
P
s
q
∘
F
,where “∘” denotes the correlation operation and F
~
s
q
∈
R
M
×
M denotes the qth differential image of the original image F. P
s
q
∈
R
(
2
L
+
1
)
×
(
2
L
+
1
) is a 2D correlation kernel, which contains only two nonzero elements; that is, the center pixel value is 1 and the qth pixel value is -
1. Note that, for every different q values, there is a different correlation kernel P
s
q since q determines the position of the pixel value -
1. Figure 2 gives an illustration of the q = 5th differential image of the whole image within an 11 × 11 search neighborhood.Figure 2
The introduced correlation kernel matrix and its correspondingv.By means of (13), the achieved differential image contains all the available information in (12). In fact, (13) can be straightforwardly implemented by a translation operation:(14)
F
~
s
q
=
F
v
q
-
F
,where F
v
q represents the image shifted by a coordinate vector v
q, 1
≤
q
≤
K. Here, the qth vector v
q points out the coordinate (x
q, y
q) from the origin (0
,
0) as shown in Figure 2. Note that the correlation kernel P
s
q and v
q, 1
≤
q
≤
K, are one-to-one mapping; that is, when (15)
q
=
c
2
L
+
1
+
r
+
1
,
0
≤
c
≤
2
L
,
0
≤
r
≤
2
L
,where (c
,
r) is the coordinate of “-
1” in P
s
q shown in Figure 2, then, its corresponding v
q is given by(16)
v
q
=
x
q
,
y
q
=
r
-
L
,
c
-
L
.Step 2.
Compute the “square image”F
¯
s
q, which is matrix case consistent with F
¯
s
i
,
q
(
m
) in (11).This step can be easily conducted by an element-wise square operation as follows:(17)
F
¯
s
q
=
F
~
s
q
•
F
~
s
q
,where “•” acts as an element-wise product operator.Step 3.
Compute the “Gaussian weighted Euclidean distance image”D
a
s
q, which is matrix case consistent with D
a
i
,
s
i
,
q in (10).This step can also be implemented by a correlation as follows:(18)
D
a
s
q
=
G
∘
F
¯
s
q
.
D
a
s
q
∈
R
M
×
M is the Gaussian weighted Euclidian distance corresponding to the qth pixel in the search neighborhood S for the whole image F. Equation (18) can be efficiently implemented by FFT as follows:(19)
D
a
s
q
=
I
F
F
T
F
F
T
G
∗
•
F
F
T
F
¯
s
q
.Step 4.
Compute the exponential Gaussian weighted Euclidian distance imageW
~
s
q, which is matrix case consistent with W
~
i
,
s
i
,
q in (6).Step 5.
Then, repeat Steps1
–
4 and compute all the W
~
s
q, 1
≤
s
q
≤
K.Step 6.
Compute the normalization factor imageZ, which is matrix case consistent with Z
(
i
) in (5).Step 7.
Compute the weightsW
(
s
q
), which is matrix case consistent with W
(
i
,
s
i
,
q
) in (4). This step can be implemented by an element-wise division operation.Step 8.
Compute the processedN
L
M
(
F
) in (3).


## 3.3. Symmetry Property

Note that the computational cost of Step5 can be further reduced by exploiting the symmetry property of the v
q; that is, (20)
v
q
=
-
v
K
-
q
+
1
,where v
(
K
-
q
+
1
) denotes the (
K
-
q
+
1
)th translation vector. The symmetry property is shown in Figure 3.Figure 3
The schematic representation of the symmetry property for weight calculation.Thus, when we calculate theqth weights W
~
s
q for the image, the (
K
-
q
+
1
)th weights W
~
s
(
K
-
q
+
1
) of the image are also meanwhile obtained as a secondary product, which saves the computational cost by one half.Specifically, when the kernelP
s
q is employed to yield the qth differential image F
~
s
q according to (14), the (
K
-
q
+
1
)th differential image F
~
s
(
K
-
q
+
1
) is achieved by means of (20) as follows: (21)
F
~
s
K
-
q
+
1
=
F
v
K
-
q
+
1
-
F
=
F
-
v
q
-
F
=
-
F
v
q
-
F
.Next, whenF
¯
s
q is calculated by (17), F
¯
s
(
K
-
q
+
1
) is also meanwhile achieved by(22)
F
¯
s
K
-
q
+
1
=
F
~
s
K
-
q
+
1
•
F
~
s
K
-
q
+
1
.In Steps3
-
4, only once same operation is implemented, W
~
s
q and W
~
s
(
K
-
q
+
1
) are simultaneously obtained.Intuitively, this can be observed in Figure3. By means of v
s
13, for each pixel i within the image F, the 13th differential image F
~
i
,
s
13
m (m
=
1
,
…
,
9) is available which is corresponding to difference between the red and green blue rectangle. At the same time, the negative F
~
i
,
s
109
m (m
=
1
,
…
,
9) is also available which is corresponding to the difference between the red and blue rectangle. Following Steps 3
-
4, only once operation is implemented, W
~
i
,
s
13 and W
~
i
,
s
109 are simultaneously obtained.Therefore, the necessary translation operations in terms ofv
q can be reduced from (K
-
1) to (
K
-
1
)
/
2 according to the symmetry property (see (20)). Thus, the computational cost will be reduced by one half.


## 3.4. Computational Complexity

Note also that the differences between the proposed fast algorithm and the traditional per-pixel NLM algorithm lie in Steps1
–
3. Steps 4 to 8 in the proposed algorithm include the same operations as those in the original NLM algorithm. So the computational complexity in these steps is not taken into account in the following subsection of complexity analysis.The computational complexity of the proposed algorithm is analyzed in this section. We assume that FFT and Inverse Fourier transform (IFFT) need the same computational complexity. In general, the computational complexity of two-dimensional FFT for anM
×
M image needs approximately (
5
/
9
)
M
log
2
⁡
M multiplications and (
4
/
3
)
M
log
2
⁡
M additions [18].Step 1.
ComputingF
~
s
q in (14) for all 1
≤
s
q
≤
(
K
-
1
)
/
2 requires (
K
-
1
)
/
2 shifts; therefore, M(K
-
1
)
/
2 additions are necessary.Step 2.
ComputingF
¯
s
q in (17) for all 1
≤
q
≤
(
K
-
1
)
/
2 needs (
K
-
1
)
M
/
2 multiplications.Step 3.
ComputingD
a
s
q in (19) for all 1
≤
s
q
≤
(
K
-
1
)
/
2 involves (
(
K
-
1
)
/
2
+
1
) FFTs, (
K
-
1
)
/
2 IFFTs, and (
K
-
1
)
M
/
2 multiplications for element-wise multiplications. Therefore, it needs (
5
/
9
)
M
K
log
2
⁡
M
+
(
K
-
1
)
M
/
2 multiplications and (
(
4
/
3
)
log
2
⁡
M
)
M
K additions.
The computational cost ofD
a
s
q using the proposed algorithm, the original algorithm [8], and the FFT based algorithm [11] are, respectively, given in Table 1. Both the proposed algorithm and the algorithm of [11] are superior to the original NLM algorithm in terms of the computational complexity. According to Table 1, if (23) holds, the complexity relations in (24) and (25) will also hold:(23)
M
≤
K
3
,
(24)
5
9
K
log
2
⁡
M
+
K
-
1
M
≤
15
9
K
log
2
⁡
K
+
K
+
1
M
,
(25)
4
3
K
log
2
⁡
M
+
1
2
K
-
1
M
≤
4
K
log
2
⁡
K
+
2
K
+
4
M
-
1
.Here, to satisfy (23), the minimum pixel number in the search neighborhood K
m
i
n will be M
3 and that can be used to deduce the minimal radius L
m
i
n of search region. Table 2 lists the K
m
i
n and L
m
i
n values according to different image sizes. These requirements on minimal values can be easily met in most practical cases; thus, the proposed method is superior to that of [11].Table 1
The computational complexity of different methods.
Methods
Multiplication operation number
Addition operation number
Proposed method
5
9
K
l
o
g
2
M
+
K
-
1
M
4
3
K
l
o
g
2
M
+
1
2
K
-
1
M
Method in [8]
2
k
K
M
2
k
-
1
K
M
Method in [11]
15
9
K
l
o
g
2
K
+
K
+
1
M
4
K
l
o
g
2
K
+
2
K
+
4
M
-
1Table 2
TheK
m
i
n and L
m
i
n values for different image sizes.
Image sizeM
×
M
M
K
m
i
n
L
m
i
n
2
56
×
256
65536
40.32
3
512
×
512
262144
64
4
1024
×
1024
1048576
101.59
5


## 4. Experiments

In this section, we assess the performance of the proposed algorithm by conducting experiments with the “Lena” image with size 512 × 512 (Figure4). The degradation is simulated by adding zero-mean Gaussian noise with standard deviations 10, 20, and 30, respectively. For Gaussian noise with different standard deviations, we use different size settings for search neighborhood S and block B to obtain the optimal denoising effect. For the standard deviations 10, 20, and 30, the sizes of block B are set as 3 × 3, 5 × 5, and 7 × 7 and the search neighborhoods S are set as 15 × 15, 17 × 17, and 21 × 21, respectively. For a fair comparison with the methods in [8, 11], we use the simplified weight calculation strategy based on the Euclidian distance rather than Gaussian weighted one. All the experiments are performed in the Matlab R2015a environment using a PC computer with 16 G ROM and Intel i7 CPU. Our experiments are compared with Buades’s classical per-pixel algorithm [8] and the Wang’s algorithm [11].Figure 4
The “Lena” image test in the experiment.We conduct quality assessment in terms of peak signal to noise ratio (PSNR) and structural similarity index measuring (SSIM) [19]. These quantitative metrics are defined by (26):(26)
P
S
N
R
P
,
I
=
10
l
o
g
10
P
m
a
x
2
1
/
m
n
∑
i
=
1
m
∑
j
=
1
n
P
i
,
j
-
I
i
,
j
2
,
S
S
I
M
P
,
I
=
1
m
n
2
μ
p
μ
i
+
C
1
2
σ
p
i
+
C
2
μ
p
2
+
μ
i
2
+
C
1
σ
p
2
+
σ
i
2
+
C
2
,where μ
p and μ
i are the means of the 8 × 8 square window and m and n are the image size. σ
p and σ
i are the corresponding standard deviations and σ
p
i is the corresponding covariance. C
1 and C
2 are constant and set according to [19].The experiments results are reported in Table3. We can note that the proposed algorithm leads to a very significant decrease in computation cost. Because all the methods are based on the same Euclidean distance, the three different methods yield the same values of PSNR and SSIM in this table.Table 3
The experiment results with different methods.
σ
PSNR
SSIM
Time (seconds)
Method in [8]
Method in [10]
Proposed method
10
33.35
0.89
306.74
31.51
4.86
20
29.73
0.80
373.44
34.72
6.76
30
27.96
0.77
556.84
41.07
10.74


## 5. Conclusion and Discussion

In this paper, we proposed a correlation based strategy to accelerate NLM algorithms. This method filters the whole image simultaneously instead of filtering pixel by pixel as the original NLM one. This new approach shows that both the computational complexity and the running time are greatly reduced when compared to the original NLM method [8] and the method reported in [11]. The acceleration method in [11] has to be applied using the simplified Euclidean distance but the proposed approach can overcome this limit by effectively calculating the Gaussian kernel convolution in (18). Our future work will consider its parallelization using GPU techniques and will address applications in medical image processing such as low-dose CT scanner [20–22] and also color image processing [23–25].Recently, the collaborative filtering algorithms, such as the BM3D algorithm (block matching and 3D filtering algorithm), which are extension of the general concepts of grouping and block matching from the NLM, have shown the superior denoising performance. These are described in [26] and applied in [27] to complex data sets with good results. However, the method also encounters the problem of the large computational complexity costed to group. The application of our proposed correlation based strategy can simply deal with the problem. Actually, the grouping process is consistent with the process of weights’ calculation within NLM algorithm. Thereby, our method can be employed to group. In addition, because of the existence of noise, the tradition grouping strategy is easily disturbed [26, 27]. Our method can overcome the issue to a large extent to enhance grouping in terms of the adjustment of the Gaussian filtering kernel employed in the NLM algorithm (see (18)).Eventually, referring to Table3, it can be observed that the computational time increases versus the augment of noise standard deviation. This is because, for larger noise corrupted image, a larger search neighborhood S and block B are necessary, and vice versa. Fortunately, the concept of random resampling mask has been developed in [28, 29], which can lower the noise standard deviation. Thus, the random mask prefiltering helps to reduce the computational time of these methods, which will be further studied in our future work.



---

*Source: 5483485-2016-07-13.xml*