This is (finally) a followup to my latest post!
From the previous episode…
In a previous post, we have talked of a recent paper that was investigating the locality of NLmeans. This paper by Simon Postec pointed out the following paradox:
The more NLmeans is nonlocal, the less it is a good denoising algorithm.
Indeed, this phenomenon was seldom observed and even less investigated, since for computational reasons most NLmeans (and derivatives) are experimented with small learning neighbourhoods only. So, shall we consider that NLmeans is local after all?
Experimental confirmation
In this post, we will confirm experimentally the results from Postec et al.
I will also take advantage of the post to present an implementation of NLmeans in Python, because:
 I want to learn and practice Python, and that’s a good reason per se!
 I didn’t find any code from Simon Postec or his coauthors, so let’s make a bit of Reproducible Research^{1}.
So, let’s dive into the code!
A Python implementation of NLmeans
Extracting patches
NLmeans is a patchoriented processor. Obviously, we’ll need a patch extractor. Here is a simple one, that applies mirroring on the image boundaries:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 

My personal style for implementing NLmeans is a bit unusual: I collect all the patches first and I store them in a matrix. Regular implementations usually build the patches on the fly. Building the patches on the fly has some advantages on lowermemory devices^{2}, but nowadays RAM is cheap and building the patches first allows for easier parallel processing later.
Comparing the patches
Once we have patches, NLmeans will denoise each pixel by comparing the corresponding patch with the patches of all the pixels in a learning neigbourhood centered on this pixel. The denoised value will be the weighted average of all these pixels, where the weights are given by the formula:
$P(\mathbf{x}), P(\mathbf{y})$ are the two patches to be compared, $h$ is a filtering parameter, and $Z$ a normalization constant (the sum of all the weights in the learning neighbourhood).
For simplicity, I have omitted the multiplication of teh patches by a Gaussian window (that has for effect to emphasize the center part of the patch), but in practice this yields only slightly different results.
The only difficulty is how to deal with the central pixel. In the weight formula, the central pixel will always have a weight of 1 and will reduce the efficeincy of the denoising. Thus, I will not compute a weight for this particular location, but will assign it the maximum of the weights computed inside the learning window.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 

Getting the final denoised value
The denoised value for the pixel at hand is obtained by taking the weighted sum of all the pixels in the learning window, adding the value of the central pixel (multiplied by the maximum weight in the neighbourhood) and dividing by the sum of the weights. In other words, a barycenter.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 

Et voilà !
The final routine is really simple: computingthe bounds of the learning windows, computing the weights, assigning the barycenter as output. And this is probably the most stunning part of NLmeans: its simplicity.
Experimental results
I’ve added a Gaussian noise of standard deviation 0.2 times the dynamic range to an image of Lena^{3}. Then, I’ve fixed the value of the $h$ parameter to 2 and run NLmeans for two different patch sizes: $5 \times 5$ and $7 \times 7$ pixels. Finally, for each patch size I let the learning neighbourhood size vary from $3 \times 3$ to $21 \times 21$.
The following figure shows the MSE for each patch and learning neighbourhood size. The neighbourhood size 0 depicts the original noisy image, and the horizontal axis depicts the half learning neighbourhood size. The final learning neighbourhood size is given by twice this value plus 1.
As can be seen from the graph, both patch sizes exhibit the same phenomenon: when the learning neighbourhood gets bigger than $5 \times 5$ pixels, the denoising performance starts to decrease. This is slightly worse than the result in Postec et al., but this can be due to the lack of a Gaussian window, a bad choice of $h$ or simply by chance since I used only 1 test image and 1 noise pattern.
Conclusions and future post
In this post, we have seen an example of an easy, nonoptimized implementation of NLmeans. Using this implementation, we have been able to confirm the observations from the Postec et al. paper: NLmeans efficiency is reduced when the implementation gets more nonlocal!
Show me the code!
I plan to release the full code package when the post series is complete. I would like to make another two or three posts^{4}.
Next post
In the next post, I will make a small experiment to confirm Postec’s intuition about the locality of NLmeans: the nonlocal hypothesis is more easily violated far from the pixel to be denoised.