\pagestyle{fancy}
\cleardoublepage\chapter{Non-Local Means (NL-Means)\index{non-local means|see{NL-means}}\index{NL-means}}\label{nlmeans}
\Quote{Science is always wrong. It never solves a problem without creating ten more.}{George Bernard Shaw}

\noindent NL-means\index{NL-means} was proposed by Buades et al.\ \cite{buades05}. The algorithm tries to replace a noisy pixel with a weighted average of other pixels with similar neighborhoods. The basic idea is that images contain self-similarities and many structures show up more than once in the image where we usually assume that noise in similar neighborhoods is uncorrelated. Thus, taking the average of similar neighborhoods yields a version of the image with less noise.

\section{Basic NL-Means}\label{nlmeans.basic}

Given a discrete noisy image $u = \{u(\textbf{\em i}) \; | \; \textbf{\em i} \in I\}$, the estimated value $NL[u](\textbf{\em i})$ for a pixel $\textbf{\em i} = (i,j)$ is computed as a weighted average of all the pixels in the image:
\begin{equation}
	NL[u](\textbf{\em i}) := \sum_{\textbf{\em j} \in I} w(\textbf{\em i},\textbf{\em j})u(\textbf{\em j}),
\end{equation}
where $\textbf{\em j} = (i',j') \in I$ and the usual conditions $0 \leq w(\textbf{\em i},\textbf{\em j}) \leq 1$ and $\sum_\textbf{\em j} w(\textbf{\em i},\textbf{\em j}) = 1$ are satisfied.

The similarity between two pixels $\textbf{\em i}$ and $\textbf{\em j}$ depends on the similarity of the intensity gray level vectors $u(\mathcal{N}_\textbf{\em i})$ and $u(\mathcal{N}_\textbf{\em j})$, where $\mathcal{N}_\textbf{\em i}$ and $\mathcal{N}_\textbf{\em j}$ denote the neighborhood of fixed size centered at pixel $\textbf{\em i}$ and $\textbf{\em j}$, respectively. Pixels with a similar gray level neighborhood to $u(\mathcal{N}_\textbf{\em i})$ have larger weights in the average. The weights are defined as follows:
\begin{equation}\label{nlmeans_equation}
	w(\textbf{\em i},\textbf{\em j}) := \frac{1}{Z(\textbf{\em i})} \exp\left(- \frac{\|u(\mathcal{N}_\textbf{\em i}) - u(\mathcal{N}_\textbf{\em j})\|^2_{2,a}}{\lambda^2}\right),
\end{equation}
where $Z(\textbf{\em i})$ is a normalization constant
\begin{equation}
	Z(\textbf{\em i}) := \sum_\textbf{\em j} \exp\left(- \frac{\|u(\mathcal{N}_\textbf{\em i}) - u(\mathcal{N}_\textbf{\em j})\|^2_{2,a}}{\lambda^2}\right)
\end{equation}
and $\lambda$ is the decay parameter of the exponential function. The expression $\|u(\mathcal{N}_\textbf{\em i}) - u(\mathcal{N}_\textbf{\em j})\|^2_{2,a}$ is an Euclidean distance, where $a > 0$ is the standard deviation of the Gaussian kernel. In other words, the parameter $\lambda$ in equation \eqref{nlmeans_equation} controls the weighting of the distances.

This weighting function penalizes structures which have less similarity in the geometrical configuration of their neighborhood, see Figure \ref{figure_nlmeansscheme}.

\begin{figure}[htb]
	\centering
	\includegraphics[scale=0.5]{figures/buades_morel_non-local_algorithm-1.eps}
	\caption[Scheme of NL-means strategy.]{Scheme of NL-means strategy. Similar pixel neighborhoods give a large weight, $w(\textbf{\em p},\textbf{\em q}_1)$ and $w(\textbf{\em p},\textbf{\em q}_2)$, while much different neighborhoods give a small weight, $w(\textbf{\em p},\textbf{\em q}_3)$. \textbf{Authors:} A. Buades, B. Coll and J.-M. Morel \cite{buades05}.}
	\label{figure_nlmeansscheme}
\end{figure}

\section{Modifications and Refinements to NL-means}\label{nlmeans.refinements}

The main drawback of the basic NL-means algorithm is its high complexity, a fact that was also noted by Buades et al.\ in \cite{buades05}. To address this issue, many ideas were proposed to reduce the complexity of the computation and thus also to reduce the runtime of the algorithm. In the following, we introduce some of them.

\subsection{Restricting the Search Space of NL-Means\index{NL-means}}\label{nlmeans.search_space}
In its very basic form, NL-means\index{NL-means} compares every pixel of the image with the pixel that should be replaced to find a possible denoising candidate. This is computationally very expensive and thus not very practical. To overcome this problem, Buades et al.\ \cite{buades05} introduced a search window (a neighborhood) which reduces the search space dramatically.

Let us assume that $N^2$ is the number of pixels of the image. Let us further assume that the similarity window has the size $W \times W$ and the search window has size $S \times S$. In its basic form, the algorithm has complexity $W^2 \times N^2 \times N^2$. With the introduction of the search window, the algorithm has complexity $W^2 \times S^2 \times N^2$ where $S^2 \ll N^2$.

\subsection{Neighborhood Classification to Speedup NL-Means\index{NL-means}}\label{nlmeans.neighborhood_classification}
Let us now explore how neighborhood classification can help to lower the algorithmic complexity of NL-means\index{NL-means}. We look at a neighborhood classification that works with two components. The first one is based on \emph{average neighborhood gray values} and the second one is based on \emph{gradient directions}. We follow here the work proposed by Mahmoudi and Sapiro \cite{mahmoudi05}.

Intuitively, similar neighborhoods should have similar average gray values. Let $\overline{u}(\textbf{\em i})$ and $\overline{u}(\textbf{\em j})$ be the average gray value in the neighborhoods of pixels $\textbf{\em i}$ and $\textbf{\em j}$, respectively. Let furthermore $\eta_1 < 1$ and $\eta_2 > 1$ be two constants close to 1. The weight $w(\textbf{\em i},\textbf{\em j})$ has a non-zero value, i.e.\ the neighborhood is only considered, if
\begin{equation}\label{nlmeans_avg_grayvalue_frac}
	\eta_1 < \frac{\overline{u}(\textbf{\em i})}{\overline{u}(\textbf{\em j})} < \eta_2
\end{equation}
holds.

A second method to approximate the similarity between two neighborhoods is their \emph{average gradient orientation}. If $\nabla u(\textbf{\em i}) = (u_x(\textbf{\em i}),u_y(\textbf{\em i}))$ is the image gradient then the average gradient orientation around pixel $\textbf{\em i}$ is defined as follows:
\begin{equation}
	\overline{\nabla u}(\textbf{\em i}) = \left(\overline{u_x}(\textbf{\em i}),\overline{u_y}(\textbf{\em i})\right),
\end{equation}
where $\overline{u_x}(\textbf{\em i})$ and $\overline{u_y}(\textbf{\em i})$ are the average $x$ and $y$ derivatives in the neighborhood of pixel $\textbf{\em i}$. Having this, we can formulate the difference of the average gradient orientations at pixels $\textbf{\em i}$ and $\textbf{\em j}$:
\begin{equation}
	\theta(\textbf{\em i},\textbf{\em j}) = \angle (\overline{\nabla u}(\textbf{\em i}),\overline{\nabla u}(\textbf{\em j})).
\end{equation}

To find a threshold above which $\theta$ is considered as outlier, Mahmoudi et al.\ propose to use robust statistics, such as
\begin{equation}\label{nlmeans_sigma_theta}
	\sigma_\theta = 1.4826 \cdot \text{median}_{I \times I}\Big[\big||\theta(\textbf{\em i},\textbf{\em j})| - \text{median}_{I \times I}(|\theta(\textbf{\em i},\textbf{\em j})|)\big|\Big].
\end{equation}
To compute a threshold for the gradient, we do a similar computation:
\begin{equation}\label{nlmeans_sigma_gradient}
	\sigma_\nabla = 1.4826 \cdot \text{median}_{I \times I}\Big[\big||\nabla u|^2 - \text{median}_{I \times I}(|\nabla u|^2)\big|\Big].
\end{equation}
Plugging everything together, we end up with the following weighting
\begin{equation}\label{nlmeans_fast_weighting_function}
	w(\textbf{\em i},\textbf{\em j}) := \left\{\begin{array}{ll}\frac{1}{Z(\textbf{\em i})} \exp\left(- \frac{\|u(\mathcal{N}_\textbf{\em i}) - u(\mathcal{N}_\textbf{\em j})\|^2_{2,a}}{\lambda^2}\right), & \big[(\|\nabla u(\textbf{\em i})\| < \sigma_\nabla) \\ & \text{or } (\|\nabla u(\textbf{\em j})\| < \sigma_\nabla) \\ & \text{or }(|\theta(\textbf{\em i},\textbf{\em j})| < \sigma_\theta)\big] \\ & \text{and } (\eta_1 < \frac{\overline{u}(\textbf{\em i})}{\overline{u}(\textbf{\em j})} < \eta_2) \\ 0, & \text{otherwise}\end{array}\right. .
\end{equation}

\subsection{Locally Adaptive Weighting for NL-Means\index{NL-means}}\label{nlmeans.adaptive_weighting}
Since the parameter $\lambda$ in equation \eqref{nlmeans_equation} controls the weighting of the distances, it is the most important and crucial parameter of NL-means\index{NL-means}: Choosing it too low can lead to a noisy result, choosing it too high can blur fine structures. Brox and Cremers \cite{brox07} showed that there is in general no global value for $\lambda$ that is well suited for the whole image. In the following we look at an idea proposed by Zimmer et al.\ \cite{zimmer08}.

Instead of using a fixed $\lambda$, a \emph{locally adaptive function} $\Lambda(\lambda,B)$ is introduced. The authors' first proposal is to use the standard deviation $s_B$ of a block $B$. The standard deviation depends on the block structure and on the amount of noise contained in $B$. Sharp structures usually mean that there are few similarities, since sharp structures lead to a high variance. Thus, we have to smooth more there compared to blocks where the variance is low. To avoid a too high influence on the smoothing, Zimmer et al.\ propose to use a sub-linear function $\Psi$. Finally, we end up with
\begin{equation}\label{nlmeans_lambda_adaptive_weighting}
	\Lambda(\lambda,B) := \lambda \cdot \Psi(s_B),
\end{equation}
where $\Psi(s_B) = \sqrt{s_B}$. Since this locally adaptive function replaces $\lambda$ in equation \eqref{nlmeans_equation}, we refer to it by a formulation like \emph{adaptively chosen (distance weight) $\lambda$} or \emph{adaptive (distance weight) $\lambda$}.