\section{Theory
	\label{section:theory} }

This section presents a short derivation of the main formulas used in the
algorithm. For a detailed analysis, see \cite{nonlocalbayesSIIMS}. Given $u$ the
noiseless ideal image and $\tilde u$ the noisy image corrupted with Gaussian noise
of standard deviation $\sigma$ so that
\begin{equation}
	\tilde u=u+n,
	\label{eq_noise}
\end{equation}
 the conditional distribution $\pro (\tilde u \mid u)$ reads
\begin{equation}
	\pro(\tilde u\mid u) = \frac{1}{(2\pi\sigma^2)^\frac{N}{2}}e^{-\frac{||u-\tilde{u}||^2}{2\sigma^2}},
\label{eq_prob}
\end{equation}
where $N$ is the total number of pixels in the image. Given a noiseless patch $P$
of $u$ with dimension $\psize\times \psize$ but seen as a column vector, and
$\tilde P$ an observed noisy version of $P$, the same model gives by the
independence of noise pixel values:
\begin{equation}
\pro(\tilde P\mid P)=c.e^{-\frac{\|\tilde P-P\|^2}{2\sigma^2}}
\label{eq_prob2}
\end{equation}
where $P$ and $\tilde P$ are considered as vectors with $\psize^2$ components and
$\|P\|$ denotes the Euclidean norm of $P$. Knowing $\tilde P$, our goal is to
deduce $P$ by maximizing $\pro(P\mid\tilde P)$.  Using Bayes' rule, we can compute
this last conditional probability as
\begin{equation}
\pro(P\mid\tilde P)=\frac{\pro(\tilde P\mid P)\pro(P)}{\pro(\tilde P)}.
\label{eq_bayes}
\end{equation}
$\tilde P$ being observed, this formula could be used to  deduce the patch $P$
maximizing the right term, viewed as a function of $P$. This is unfortunately not
possible, unless we have a probability model for $P$.   We shall now discuss how to
proceed when all observed patches are noisy. Assume that the patches $Q$ similar to
$P$ follow a Gaussian model with covariance matrix $\textbf{C}_P$ and mean
$\overline{P}$. This means that
\begin{equation}
	\pro(Q)= c. e^{-\frac{(Q-\overline{P})^t\textbf{C}_P^{-1}	(Q-\overline{P})}{2}}
\end{equation}
From \eqref{eq_prob} and \eqref{eq_bayes} we obtain for each observed $\tilde P$
the following equivalence of problems:
\begin{eqnarray*}
\max_{P}\ \pro(P\mid\tilde P) & \Leftrightarrow &\max_{P}\ \pro(\tilde P\mid P)\pro(P) \\
& \Leftrightarrow & \max_{P}e^{-\frac{\|P-\tilde P\|^2}{2\sigma^2}}
e^{-\frac{(P-\overline{P})^t\textbf{C}_P^{-1}(P-\overline{P})}{2}}\\
& \Leftrightarrow & \min_P \ \frac{\|P-\tilde P\|^2}{\sigma^2} + (P-\overline{P})^t\textbf{C}_P^{-1} (P-\overline{P})\\
\end{eqnarray*}
where $\bC_P$ denotes the covariance matrix of the patches similar to $P$ and
$\overline{P}$ the expectation of the patches similar to $P$. Nevertheless, this
expression is not satisfactory, the noiseless patch $P$ and the patches similar to
$P$ not being  observable. Since we are observing the noisy version $\tilde P$, the
patch-based algorithms can at least compute the patches $\tilde Q$ similar to
$\tilde P$. An empirical covariance matrix can therefore be obtained with enough
such observable samples $\tilde Q$. In theory we have from\eqref{eq_noise}, $P$ and
$n$ being independent,
\begin{equation}
\bC_{\tilde P} = \textbf{C}_P + \sigma^2\textbf{I};\quad\quad E{\tilde Q}= \overline{P}.
\label{eq_cov}
\end{equation}
Notice that these relations assume that we searched for patches similar to $\tilde
P$ in a neighbourhood large enough to include all patches similar to $P$, but not
too large either, because otherwise it can contain outliers. Thus a safe strategy
is to search similar patches in a distance slightly larger  than the plausible
distance caused by noise. If the above estimates are correct, the MAP (maximum
\emph{a posteriori} estimation) problem finally boils down by \eqref{eq_cov} to the
following feasible  minimization problem:
$$
\max_{P}\ \pro(P\mid\tilde P)  \Leftrightarrow \min_P \ \frac{\|P-\tilde P\|^2}{\sigma^2} + (P-\overline{\tilde P})^t(\textbf{C}_{\tilde P}-\sigma^2\textbf{I})^{-1} (P-\overline{\tilde P}).
$$
Differentiating this quadratic function with respect to $P$ and equating to zero
yields
$$
P-\tilde P+\sigma^2(\textbf{C}_{\tilde P}-\sigma^2\textbf{I})^{-1}(P-\overline{\tilde P})=0.
$$
Taking into account that $\bI+\sigma^2 (\bC_{\tilde P} - \sigma^2 \bI)^{-1} = (\bC_{\tilde P} - \sigma^2 \bI)^{-1} \bC_{\tilde P},$ this yields
$$
(\bC_{\tilde P}-\sigma^2\bI)^{-1}\bC_{\tilde P} P=\tilde P +\sigma^2(C_{\tilde P}-\sigma^2 \bI)^{-1}\overline{\tilde P}.
$$
and therefore

\begin{eqnarray*}
P
&=& \bC_{\tilde P}^{-1}(\bC_{\tilde P}-\sigma^2\bI)\tilde P+\sigma^2\bC_{\tilde P}^{-1}\overline{\tilde P}\\
&=& \overline{\tilde  P}+\left[\bC_{\tilde P}-\sigma^2\bI\right]\bC_{\tilde  P}^{-1}(\tilde P-\overline{\tilde P})
\end{eqnarray*}

Thus  a restored patch $P^{\basic}$ can be obtained from the observed patch $\tilde
P$ by the one step estimation
\begin{equation}
P^{\basic}=\overline{\tilde  P}+\left[\bC_{\tilde P}-\sigma^2\bI\right]\bC_{\tilde  P}^{-1}(\tilde P-\overline{\tilde P}).
\label{eq_step1}
\end{equation}
In a second step, where all patches have been denoised at the first step, all the
denoised patches can be used again to obtain an estimation $\bC_{\tilde
P}^{\basic}$ for $\bC_P$, the covariance of the cluster containing $P$, and
$\overline{\tilde P}^{\basic}$ a new estimation of $\overline{\tilde P}$, the
average of patches similar to $\tilde P$. Indeed, the patch similarity is better
estimated with the denoised patches. Then it follows from \eqref{eq_cov} and
\eqref{eq_step1} that we can obtain a second better denoised patch,
\begin{equation}
P^{\final} = \overline{\tilde P}^{\basic} + \bC_{\tilde P}^{\basic} \left[\bC_{\tilde P}^{\basic} + \sigma^2 \bI \right]^{-1}(\tilde P-\overline{\tilde P}^{\basic}) \label{eq_step2}
\end{equation}
The way of computing $\bC_{\tilde P}$ and $\overline{\tilde P}$ will be discussed
in the section \ref{section:implementation}.
