\section{Implementation
	\label{section:implementation} }
	

\subsection{Reminder of BM3D}


The algorithm developed in this paper is very similar to BM3D and also has two
successive steps. In order to use this similarity as much as possible our notation
and exposition order will be as close as possible to those used in
\cite{bm3d_ipol}. Here is a brief overview of BM3D:

\begin{enumerate}
	\item {\bf Step 1: first denoising\\} Loop on the image. We denote by
$\tilde P$ the reference  current noisy patch.
	\begin{enumerate}
		\item \underline{Grouping:\\} Stacking up similar patches to the reference one, using a similarity threshold applied to the distance between patches in order to build the 3D block $\mP(\tilde P)$;
		\item \underline{Collaborative Filtering:\\} A 3D linear transform
is applied to the 3D block, then a hard thresholding is applied to the
coefficients, and finally the inverse 3D transform is applied. A weight is
associated with the whole 3D block, depending to its sparsity;
		\item \underline{Aggregation:\\} A first basic estimate of the
denoised image is obtained by doing a weighted aggregation of every
estimate obtained in the preceding step for each pixel. This basic estimate
will be denoted by $u^{\basic}$.
	\end{enumerate}
	\item {\bf Step 2: second denoising step using the result of the first as
``oracle''\\}
	\begin{enumerate}
		\item \underline{Grouping: \\}
		The distance between patches is computed on the basic estimate. Two 3D blocks are formed :
		\begin{itemize}
			\item $\mP^{\basic}(P^{\basic})$ by stacking up patches
from the basic estimation $u^{\basic}$ and,
			\item $\mP^{\basic}(\tilde P)$ by stacking up patches in
the same order  from the original noisy image $\tilde u$.
		\end{itemize}
		\item \underline{Collaborative Filtering: \\}
		A 3D transform is applied on both 3D blocks, followed by a Wiener
filtering of the group $\mP^{\basic}(\tilde P)$ using the empirical
oracular coefficients obtained from  the group $\mP^{\basic}(P^{\basic})$,
and finally by the inverse 3D transform. A weight is computed for the whole
3D block. It depends on  the norm of the empirical Wiener coefficients;
		\item \underline{Aggregation:\\}
		A final estimate of the denoised image is obtained by using a
weighted aggregation of every estimate obtained for each pixel. This final
estimate will be denoted by $u^{\final}$.
		\end{enumerate}
\end{enumerate}


\subsection{Comparison of the structure of NL-Bayes with BM3D}


All algorithms will be described for color images. They can also be applied on grey
level images; the changes in that case will be indicated. Like BM3D, NL-Bayes
 is applied in two successive steps, the result of the first one serving as oracle for the second one:
\begin{enumerate}
\vspace{-3mm}\item the first step provides a basic estimate $u^{\basic}$ by
    using \eqref{eq_step1} during the {\it collaborative filtering}. Parameters
    in this step are denoted by the exponent \hard; \vspace{-3mm}\item the
second step is based both on the original noisy image $\tilde u$ and
    on the basic estimate obtained during the first step $u^{\basic}$ in order
    to apply \eqref{eq_step2} during the {\it collaborative filtering}.
    Parameters in this step are denoted by the exponent \wien.
\end{enumerate}




The following table allows to compare steps between BM3D and NL-Bayes for color
images.

\begin{changemargin}{-2cm}{0cm}
\begin{center}
	\begin{tabular}{||c||c|c||}
		\hline
		\hline
		\multicolumn{3}{||c||}{{\bf Step 1}}\\
		\hline
		 & BM3D & NL-Bayes\\
		 \hline
		 \hline
		 {\bf Preprocessing} & Transformation to the $Y_0U_0V_0$ color space & Transformation to the $Y_0U_0V_0$ color space \\
		 \hline
		 Scanning step between two patches & $p_{\hard}=3$ & $p_{\hard}=1$\\
		 Processing  already used patches & Yes & No\\
		 \hline
		 {\bf Grouping} & & \\
		 Distance between patches & Channel $Y_0$ & Channel $Y_0$\\
		  & Normalized quadratic distance & Normalized quadratic distance\\
		 Similarity threshold & Fixed, tabulated according to $\sigma$ & -\\
		 Patches kept & $N_{\hard}$ best & $N_{\hard}$ best\\
		 Patches ordered & No & No\\
		 3D groups formed & One for each channel $Y_0$, $U_0$ et $V_0$ &
 one for each channel $Y_0$, $U_0$ et $V_0$\\
		 \hline
		 {\bf Collaborative Filtering} & & \\
		 3D transform & 2D Bior1.5 on each patch followed by & -\\
		  & 1D Hadamard transform along the  & \\
		  & third dimension & \\
		 Filter & Hard thresholding on the coefficients & Bayesian, based on \eqref{eq_estimate1}\\
		  & of the DCT & \\
		 Weighting & Depending on the number of non-zero & -\\
		  & coefficients after Hard thresholding & \\
		 \hline
		 {\bf Aggregation} & \multicolumn{2}{c||}{Identical part}\\
		 \hline
		 {\bf Post processing} & - & Transform to the RGB color space\\
		 \hline
		 \hline
		 \multicolumn{3}{||c||}{{\bf Step 2}}\\
		 \hline
		  & BM3D & NL-Bayes\\
		 \hline
		 \hline
		 Step between two patches & $p_{\wien}=3$ & $p_{\wien}=1$\\
		 Process of an already used patch & Yes & No\\
		 \hline
		 {\bf Grouping} & & \\
		 Distance & Channel $Y_0$ of $u^{\basic}$ & All channels of $u^{\basic}$\\
		  & Normalized quadratic distance & Normalized quadratic distance\\
		 Similarity threshold & fixed, tabulated according to the $\sigma$ & Adaptive according to the distance\\
		  & & of the $N_{\wien}$-th best patch\\
		 Patches kept & $N_{\hard}$ best & All\\
		 Patches ordered & No & No\\
		 3D groups formed & Two for each channel & Two\\
		 \hline
		 {\bf Collaborative Filtering} & & \\
		 3D transform & 2D DCT then 1D Hadamard & -\\
		  & transform on both groups & \\
		 Filter & Wiener filter using $u^{\basic}$ as oracle & Bayesian, based on \eqref{eq_estimate2}\\
		 Weighting & Depending on the norm of the & -\\
		  & empirical Wiener coefficients & \\
		 \hline
		 {\bf Aggregation} & \multicolumn{2}{c||}{Identical part}\\
		 \hline
		 {\bf Post processing} & Transformation to the RGB color space & - \\
		 \hline
	\end{tabular}
\end{center}
\end{changemargin}


\subsection{The first step of the NL-Bayes}


Only for the first step, the noisy image $\tilde u$ in the usual RGB color space is
converted in a different color space where the independent denoising of each
channel will not create noticeable color artifacts. Most algorithms use the $YUV$
system which separates the geometric and chromatic parts of the image. As for BM3D,
it is a linear transform multiplying the $RGB$ vector by the matrix
$$
Y_oU_oV_o = \begin{pmatrix}
		\frac{1}{3} &\frac{1}{3} & \frac{1}{3}\\
		\frac{1}{2} & 0 & -\frac{1}{2} \\
		\frac{1}{4} & -\frac{1}{2} & \frac{1}{4}\\
		\end{pmatrix}
$$
We wrote the matrix above without normalization for readability, but  the matrix is
normalized to become orthonormal. In that way, $\sigma$ still is the value of the
standard deviation of the noise on each channel $Y_0$, $U_0$ and $V_0$. The
transform increases the SNR of the geometric component, the $Y_0$ component being
an average of the three colors. The geometric component is perceptually more
important than the chromatic ones, and the presence of less noise permits a better
performance of the algorithm in this component. The components $U_0$ and $V_0$ are
differences of channels, which cancel or attenuate the signal. Thus a higher noise
reduction on the chromatic components $U_0$ and $V_0$ is possible, due to their
observable regularity.
%
We denote by $\tilde P$ the current reference patch with size $\psize_{\hard}\times
\psize_{\hard}$ (seen as a column vector) of the noisy image $\tilde u$.

\begin{enumerate}
	\vspace{-3mm}\item \underline{Grouping:}\\
	  The original noisy image $\tilde u$ is searched in a $\tilde P$-centered
$\zsize_{\hard}\times \zsize_{\hard}$ neighbourhood for patches $\tilde Q$
similar to the reference patch $\tilde P$. The normalized quadratic distance
between each patch $\tilde Q$ of the neighbourhood and the reference patch
$\tilde P$ is computed as
$$
d^2\left(\tilde P, \tilde Q\right) = \frac{\|\tilde P - \tilde Q\|_2^2}{\left(\psize_{\hard}\right)^2}.
$$
This distance is computed  on the luminance channel $Y_0$ only. All patches
$\tilde Q$ of the neighbourhood are sorted according to their distance to the
reference patch $\tilde P$, and the $N_{\hard}$ closest patches to $\tilde P$
are kept. Then three sets of similar patches -one for each channel- are built:
$\mathcal{P}_{Y_0}(\tilde P)$, $\mathcal{P}_{U_0}(\tilde P)$ and
$\mathcal{P}_{V_0}(\tilde P)$. But the sets of similar patches for the
chromatic channels are built with patches whose index are the same as for
$\tilde Q\in \mathcal{P}_{Y_0}(\tilde P)$, and in the same order. After this
step, the same procedure is applied on each channel, but separately. For a sake
of simplicity, the distinction between channels will not be done any more, and
we will work for a generic channel.
	\vspace{-3mm}\item \underline{Collaborative Filtering:}\\
Let $\mathcal{P}(\tilde P)$ be the set of patches $\tilde Q$ similar to the
reference patch $\tilde P$ obtained at the grouping step. We  start by
detecting if $\tilde P$ belongs to a homogeneous area by processing the square
of the standard deviation of $\mathcal{P}(\tilde P)$:
	\begin{equation}
	\sigma_{\tilde P}^2 = \frac{M_{\hard}}{M_{\hard} - 1}\left(\frac{1}{M_{\hard}} \sum_{\tilde Q\in \mathcal{P}(\tilde P)}\sum_{x\in \tilde Q}\left(\tilde{Q}(x)\right)^2 - \left(\frac{1}{M_{\hard}}\sum_{\tilde Q\in\mathcal{P}(\tilde P)}\sum_{x\in \tilde Q}\tilde{Q}(x)\right)^2\right)
	\label{eq_gamma}
	\end{equation}
	where $M_{\hard} = N_{\hard}(\psize_{\hard})^2$. Since  a huge number
($M_{\hard}$) of realizations of the variable $u(i)$ is taken  into account, in
a homogeneous area this random variable should be very  concentrated around its
mean. Thus, fixing a threshold $\gamma$ close to 1,
	\begin{itemize}
		\vspace{-3mm}\item if $\sigma_{\tilde P}\leq \gamma \sigma$,  we
can assume that with high probability $\tilde P$  belongs to a homogeneous
area. In this case, the better result that can be obtained for the group is
the average. Therefore, the estimate of all patches in the set of similar
patches $\mathcal{P}(\tilde P)$ is
%	$$
%	\forall \tilde Q\in \mathcal{P}(\tilde P), \forall i\in (k_{\hard})^2, Q^{\basic}_i = \frac{1}{M_{\hard}}\sum_{\tilde Q\in \mathcal{P}(\tilde P)}\sum_{j=1}^{(k_{\hard})^2}\tilde{Q}_j
%	$$
	$$
	\forall \tilde Q\in \mathcal{P}(\tilde P), \forall x\in Q^{\basic}, Q^{\basic}(x) = \frac{1}{M_{\hard}}\sum_{\tilde Q\in \mathcal{P}(\tilde P)}\sum_{y\in \tilde Q}\tilde{Q}(y)
	$$
		\vspace{-3mm}\item else, $\tilde P$ is assumed to contain some
signal, for which a Gaussian model is built. By the law of large numbers we
have
		\begin{equation}
\bC_{\tilde P}\! \simeq \! \frac 1{\# \mathcal{P}(\tilde P) - 1}\!\!
\sum_{\tilde Q\in\mathcal{P}(\tilde P)}\!\!\!\!\left(\tilde Q - \overline{\tilde P}\right)\!\!
\left(\tilde Q-\overline{\tilde P}\right)^t\!\!, \quad \;\;
\overline {\tilde P}\!\simeq \! \frac 1{\# \mathcal{P}(\tilde P)}
\!\!\!\sum_{\tilde Q\in\mathcal{P}(\tilde P)}\tilde Q.				\label{eq_cov1}
	\end{equation} Once $\bC_{\tilde P}$ and $\overline{\tilde P}$ have
been computed, \eqref{eq_step1} yields an estimate for every patch in the
set of similar patches,
	\begin{equation}
	\forall \tilde Q\in \mathcal{P}(\tilde P), Q^{\basic} = \overline{\tilde P} + \left[\bC_{\tilde P} - \beta_{\hard} \sigma^2\bI\right]\bC_{\tilde P}^{-1}\left(\tilde Q - \overline{\tilde P}\right).
	\label{eq_estimate1}
	\end{equation}
	\end{itemize}
	where $\beta_{\hard}$ is a parameter of conservative attenuation close to
1.
	\vspace{-3mm}\item \underline{Aggregation:}\\
When the collaborative filtering is achieved, an estimate is associated with
every used patch. This yields a variable number of estimates for each pixel. To
take advantage of these multiple estimates an aggregation must be done.
Contrarily to BM3D, this aggregation is not weighted. The final estimate after
this first step is given by
	$$
	u^{\basic}(x) = \frac{\displaystyle{\sum_{\tilde P}\sum_{\tilde Q \in\mathcal{P}(\tilde P)} \chi_{\tilde Q}(x)Q^{\basic}(x)}}{\displaystyle{\sum_{\tilde P}\sum_{\tilde Q \in\mathcal{P}(\tilde P)}\chi_{\tilde Q}(x)}}
	$$
	with $\chi_{\tilde Q}(x) = 1$ if and only if $x\in\tilde Q$, $0$ otherwise.
	
	
	\underline{Remark:} For practical purposes, this computation is simplified by using two buffers $\bm{\nu}$ and $\bm{\delta}$, where respectively the numerator and the denominator are kept in memory
	$$\forall \tilde Q \in \mathcal{P}(\tilde P), \forall x \in \tilde Q,
\left\{
	\begin{array}{ll}
		\bm{\nu}(x) = \bm{\nu}(x) + Q^{\basic}(x)\\
		\bm{\delta}(x) = \bm{\delta}(x) + 1\\
	\end{array}
\right.$$
	Thus the final estimate is simply obtained by dividing both buffers
element-by-element.

\item Acceleration: To speed up the algorithm, each patch that has been used
    (and therefore denoised at least once) in a 3D group is no more considered
    as reference patch $\tilde P$. Nevertheless, it may be denoised several
    times, being potentially chosen in other groups.
\end{enumerate}
Once $u_{Y_0}^{\basic}$, $u_{U_0}^{\basic}$ and $u_{V_0}^{\basic}$ have been
obtained, inverting the color transform yields back $u^{\basic}$, the first
estimate of the image in the $RGB$ color space.

\subsection{Second step of NL-Bayes}


In this second step of the algorithm a basic estimate $u^{\basic}$ of the denoised
image is available. The second step follows exactly the same scheme as the first,
but performs a Wiener filter of the original noisy image $\tilde u$,  using as
oracle the basic estimate $u^{\basic}$.


\begin{enumerate}
	\item \underline{Grouping:}\\
The patch matching is  processed on the basic estimate only. But this time the
distance involves all channels, which are assumedly denoised by the first step:
	\begin{equation}
	\forall Q^{\basic}, d^2(P^{\basic}, Q^{\basic})=\frac{1}{N_c}\sum_{c=1}^{N_c}\frac{\|P^{\basic}_c-Q^{\basic}_c\|_2^2}{(k_{\wien})^2}
	\label{eq_dist2}
	\end{equation}
	where $N_c$ denotes the number of channels in the image. As a difference
with the first step where only $N_{\hard}$ patches were kept, here a threshold
$\tau$ is used to obtain a set of similar patches
	$$
	\mP^{\basic}(P^{\basic})=\{Q^{\basic}:d^2(P^{\basic}, Q^{\basic})\leq\tau\},
	$$
	with
	\begin{itemize}
		\vspace{-3mm}\item $\tau = \max\left(\tau_0, d_{N_{\wien}}\right)$;
		\vspace{-3mm}\item $\tau_0$ is a fixed parameter;
		\vspace{-3mm}\item $d_{N_{\wien}}$ is the distance between the
reference patch and its $N_{\wien}$-th best similar patches, sorted by
their distance to $P_{\wien}$.
	\end{itemize} Thus, using $\tau$, many more similar patches can be
picked in homogeneous areas. 		A second set of similar patches from
the noisy image $\tilde u$ is built
	$$
	\mP^{\basic}(\tilde P) = \{\tilde Q:d^2(P^{\basic}, Q^{\basic})\leq\tau\},
	$$ by stacking up patches together in the same order as
$\mathcal{P^{\basic}}(\tilde P)$.
	\item \underline{Collaborative Filtering:}\\
	Once 3D-blocks are built the collaborative filtering is applied. Then by
the law of large numbers,

	\begin{equation}
	\begin{array}{c}
	\displaystyle{\bC_{P}^{\basic}\! \simeq \! \frac 1{\# \mP^{\basic}(P^{\basic}) - 1}\!\! \sum_{Q^{\basic} \in \mP^{\basic} (P^{\basic})}\!\!\!\!\! \left(Q^{\basic} - \overline{\tilde P}^{\basic}\right)\!\!
\left(Q^{\basic}-\overline{\tilde P}^{\basic}\right)^t\!\!}, \\
 \\
	\displaystyle{\overline {\tilde P}^{\basic}\!\simeq \! \frac 1{\# \mP^{\basic}(\tilde P)}\!\!\!
\sum_{\tilde Q\in\mP^{\basic}(\tilde P)}\tilde Q}.
	\end{array}
	\label{eq_cov2}
	\end{equation} Once $\bC_{\tilde P}^{\basic}$ and $\overline{\tilde
P}^{\basic}$ are computed,  \eqref{eq_step2} yields an estimate for every patch
in the set of similar patches 	
	\begin{equation}
	\forall \tilde Q\in \mP^{\basic}(\tilde P), Q^{\final} = \overline{\tilde P}^{\basic} + \bC_{\tilde P}^{\basic} \left[\bC_{\tilde P}^{\basic} + \beta_{\wien} \sigma^2\bI\right]^{-1} \left(\tilde Q - \overline{\tilde P}^{\basic}\right)
	\label{eq_estimate2}
	\end{equation}
	
	\item \underline{Aggregation:}\\ When the collaborative filtering is
achieved, an estimate is associated with every used patch and therefore a
variable number of estimates for every pixel. Once again, contrarily to BM3D,
this aggregation is not weighted. The final estimate after this second step is
given by
	$$
	u^{\final}(x) = \frac{\displaystyle{\sum_{\tilde P}\sum_{\tilde Q\in\mP^{\basic}(\tilde P)} \chi_{\tilde Q}(x) Q^{\final}(x)}} {\displaystyle{\sum_{\tilde P}\sum_{\tilde Q\in\mP^{\basic}(\tilde P)} \chi_{\tilde Q}(x)}}
	$$
	with $\chi_{\tilde Q}(x) = 1$ if and only if $x\in\tilde Q$, $0$ otherwise.
	
	
	\underline{Remark:} For practical purposes, this computation is simplified by using two buffers $\bm{\nu}$ and $\bm{\delta}$, where respectively the numerator and the denominator are kept in memory
	$$\forall \tilde Q \in \mP^{\basic}(\tilde P), \forall x \in \tilde Q,
\left\{
	\begin{array}{ll}
		\bm{\nu}(x) = \bm{\nu}(x) + Q^{\final}(x)\\
		\bm{\delta}(x) = \bm{\delta}(x) + 1\\
	\end{array}
\right.$$
	Thus the final estimate is simply obtained by dividing both buffers element-by-element.
	 Once again, in order to speed up the algorithm, each patch used once in a
3D group is no more processed as reference patch. It can be used anyway several
times as secondary patch in other 3D blocks.
\end{enumerate}


\subsection{NL-PCA: a particular case of NL-Bayes
	\label{subsec:nl_pca} }


From this Bayesian approach on which NL-Bayes is based, another  algorithm can be
formulated, which we shall call NL-PCA. The idea is to compute a Principal
Component Analysis (PCA) on the 3D group $\mP(\tilde P)$. This idea was first
proposed in \cite{tsid} and is also studied in detail in \cite{Deledalle2011}.
NL-PCA is obtained by replacing in BM3D the fixed linear transform  (DCT,
bi-orthogonal spline wavelet) by a learnt basis for each patch, obtained by PCA on
the patches of the 3D block.  Indeed, applying a PCA to the 3D group amounts to
diagonalize its covariance matrix, and the eigenvectors give the adaptive basis. In
continuation, the collaborative filtering and the aggregation parts can be applied
exactly like in BM3D.The algorithm developed in \cite{Deledalle2011} is very close
to the first step of the NL-PCA algorithm described hereafter. This article will be
described and commented more in detail in section \ref{section:analysis_biblio}.


Diagonalizing the covariance matrix $\bC_{\tilde P}$, denoting the associated
isometry by $R_{\hard}$ and denoting by $S_{\hard}(j)$ the squares of the
associated eigenvalues, the restoration formula \eqref{eq_step1} becomes on the
eigenfunction basis:
$$
\left(R_{\hard}(P^{\basic} - \overline{\tilde P})\right)_j =
\frac{S_{\hard}(j) - \sigma^2}{S_{\hard}(j)}\left(R_{\hard}(\tilde P - \overline{\tilde P})\right)_j.
$$
The only (classic) variation with respect to this estimate is  that $\bC_{\tilde
P}$ should {\bf be} positive semi-definite. Thus $S_{\hard}(j)-\sigma^2$ is
replaced in the above formula by $\left(S_{\hard}(j) - \sigma^2\right)^+$ and,
instead of $-\sigma^2$, a more conservative attenuation is applied,
$-\beta^2\sigma^2$ where an empirical $\beta$ slightly larger than $1$ accounts for
the error of model. A still more conservative estimate is applied for large noises
(typically if $\sigma\leq 40$), where the estimate becomes
$$
\left(R_{\hard}(P^{\basic} - \overline{\tilde P})\right)_j=
\begin{cases} \left(R_{\hard}(\tilde P-\overline{\tilde P})\right)_j & \hbox{if}\; S_{\hard}(j)\geq \beta^2\sigma^2\\
0 & \hbox{otherwise}.
\end{cases}
$$
In the same way, in the second step, diagonalizing $\bC_{\tilde P}^{\basic}$ and
denoting the associated isometry by $R_{\wien}$ and the squares of the associated
eigenvalues by $S_{\wien}(j)$, the restoration formula \eqref{eq_step2} becomes,
without any alteration to the model,
$$
\left(R_{\wien}(P^{\final}-\overline{\tilde P}^{\basic})\right)_j = \frac{S_{\wien}(j)}{S_{\wien}(j) + \sigma^2} \left(R_{\wien}(\tilde P - \overline{\tilde P}^{\basic}\right)_j
$$
which retrieves exactly a classical Wiener filter based on the PCA of $\mP^{\basic}(\tilde P)$.
