%\documentclass[english,a4paper,11pt]{article}
\documentclass[english,a4paper,11pt]{book}
%\usepackage[english]{babel}
\usepackage[T1]{fontenc}
\usepackage{ucs}
\usepackage[utf8x]{inputenc}
\setcounter{secnumdepth}{3}
\setcounter{tocdepth}{1}
\usepackage{float}
\usepackage{graphicx}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{amsthm}
\usepackage{bbm}
\usepackage{listings}
\usepackage{subfigure}
\usepackage{pst-node}
\usepackage{pstricks}
\usepackage{mathrsfs}
\usepackage[active]{srcltx}
\usepackage[ruled,vlined,linesnumbered]{algorithm2e}
\usepackage{color}
\bibliographystyle{amsalpha}
\lstloadlanguages{C++}
%\setlength{\oddsidemargin}{-0.5 cm}
%\setlength{\textwidth}{17 cm}
\newcommand{\pd}[2]{\frac{\partial#1}{\partial#2}}
\newcommand{\pdd}[2]{\frac{\partial^2#1}{\partial#2 ^2}}
\newcommand{\pddd}[3]{\frac{\partial^2#1}{\partial#2\partial#3}}
\newcommand{\bb}[1]{\mathbb #1}
\newcommand{\cl}[1]{\mathcal #1}
\newcommand{\mb}[1]{\mathbf #1}
\newcommand{\epsi}{\varepsilon}
\newcommand{\phii}{\varphi}
\newcommand{\ee}{\text{  }}
\newcommand{\cotan}{\operatorname{cotan}}
\graphicspath{{images/},{figures/}}%{../21_Juillet_2009/figures/},{../22_Aout_2009/figures/}}
\newtheorem{theorem}{Theorem}
\newtheorem{lemma}{Lemma}
\newtheorem{acknowledgement}{Acknowledegment}
\newtheorem{julie}{Julie's remark}
\newcommand{\E}{\operatorname{E}}
\newtheorem{Jean-Michel}{Jean-Michel:}
\title{Neighborhood filters and the recovery of 3D information
\label{Neighborhood_filters_and_the_recovery_of_3D_information}}

 

\author{Julie Digne$^{1}$, Mariella Dimicoli$^{2}$, Neus Sabater$^{3}$ and  Philippe Salembier$^{3}$
       \\$^1$  CMLA, Ecole Normale Sup\'{e}rieure de Cachan,\\
 julie.digne@cmla.ens-cachan.fr\\
         $^2$Laboratory of Physiology of Perception and Action
\\ Collège de France\\
maria.dimiccoli@college-de-france.fr\\
 $^3$  CMLA, Ecole Normale Sup\'{e}rieure de Cachan,\\
 neus.sabater@cmla.ens-cachan.fr\\
F-94235 Cachan\\
 $^4$ Image and Video Processing Group\\
Technical University of Catalonia (UPC)\\
philippe@gps.tsc.upc.edu}
 

\begin{document}
 \maketitle
 
 \paragraph{Abstract}{\it Following their  success in image processing (see Chapter \ref{neighborhoodfilterschapter}), neighborhood filters have been extended to 3D surface processing. This adaptation is not straightforward. It has led to  several variants for surfaces depending on whether the surface is defined as a mesh, or as a raw data set point. The image grey level in the bilateral similarity measure is replaced by a geometric information such as the normal or the curvature. The first section of this chapter reviews the variants of 3D mesh bilateral filters and compares them to the simplest possible isotropic filter, the mean curvature motion.

In a second part, this chapter reviews applications of the bilateral filter to a data composed of a sparse depth map (or of depth clues), and of the image on which they have been computed. Such sparse depth clues can be obtained by stereo vision or by psychophysical techniques.
 The  underlying assumption to these applications is that pixels with similar intensity around a region are likely to have similar depths. Therefore, when diffusing depth information  with a bilateral filter  based  on locality and color similarity, the discontinuities in depth are assured to be consistent with the color discontinuities, which is generally a desirable property. In the reviewed applications, this ends up with the reconstruction of a dense perceptual depth map from the joint data of an image and of depth cues.  }
 \tableofcontents

 \medskip
The above abstract is enough of an introduction, and our plan follows from it. Section \ref{meshesbilateral} reviews bilateral filters applied to 3D data set points, often organized in a triangulation (a mesh). It  ends up with comparative simulations illustrating the advantage of bilateral filters on isotropic filtering. Section \ref{sec:relief} considers the various cases where in an image depth values or depth cues are available, and shows that the bilateral filter used as a diffusion tool makes a good job of restoring a dense depth map.
\section{Bilateral filters processing of meshed 3D surfaces\label{meshesbilateral}}
  This section proceeds by first examining the various adaptations of bilateral filtering on meshes (triangulated 3D surfaces)  and discussing their implementation, which can depend on the surface triangulation. Then it will be shown that the main bilateral filter is actually consistent on smooth surfaces with the mean curvature motion (although it behaves differently on edges).  Finally several comparative experiments on synthetic and real meshes will be performed.
Since a common notation is needed for all methods, we  begin with a small glossary and notation summary to which the reader may refer in this section.

\paragraph{Glossary and notation}
\begin{itemize}
\vspace{-3mm}
\item $\mathcal M$ : the mesh, namely a set of triangles
\vspace{-3mm}
\item $v$ current mesh vertex to be denoised
\vspace{-3mm}\item $p\in\cl N(v)$ neighborhood of the vertex made of other vertices $p$
\vspace{-3mm}\item $n_v$, $n_p$, etc.  normals at vertex $v$ or point $s$, etc.
\vspace{-3mm}\item $w_1(||p-v||)$, $w_2(<n_v, p-v>)$, etc.: 1D centered Gaussians  with various variances, used as weighting functions applied to the distance of neighbors  to the current vertex and to the distance along the normal direction at $v$.
\vspace{-3mm}\item $H_v$, $H_p$, etc. curvatures at $v$, $p$, etc.
    \vspace{-3mm}\item $f$ triangle of a mesh
    \vspace{-3mm}\item $a_f$ area of triangle $f$
    \vspace{-3mm}\item $c_f$ barycenter of the triangle $f$
    \vspace{-3mm}\item $n_f$ normal to triangle $f$
    \vspace{-3mm}\item $\Pi_f$ projection on the plane of triangle $f$
    \vspace{-3mm}\item $V$ voxel containing points of the data set point
    \vspace{-3mm}\item $s'$, $v'$, $p'$, $n_v'$ processed versions of $s$, $v$, $p$, $n_v$, ...
    \end{itemize}

\subsection*{Main adaptations proposed for meshes}
The {\it neighborhood filter or sigma filter} is  attributed to J.S. Lee \cite{lee1983digital} in 1983 but goes back to L. Yaroslavsky and
the Sovietic image processing theory (see the book summarizing these works \cite{yaroslavsky1985digital}) in 2D image analysis. A recent variant names it bilateral filter (\cite{bilatimage}). The bilateral filter denoises a pixel by using a weighted mean of its similar neighbors gray levels. In the original article the similarity measure was the difference of pixel gray levels yielding for a pixel $v$ of an image $I$ with neighborhood $\cl N(v)$:

$$\hat I(v)=\frac{\sum_{p\in \cl N(v)} w_1(\|p-v\|)w_2(|I(v)-I(p)|)I(p)}{\sum_{p\in \cl N(v)} w_1(\|p-v\|)w_2(|I(v)-I(p)|)}$$
where $w_1$ et $w_2$ are decreasing functions on $\bb R^+$ (e.g. Gaussian). Thus $\hat I(v)$ is an average of pixel values for pixels that are similar in position, but also in value. Hence the ``bilaterality''.
 \medskip
Filtering without losing the  sharp features  is as critical for surfaces as it is for images, and a first adaptation of the bilateral filter to surface meshes was proposed in \cite{bilateral}.
Consider a meshed surface $\cl M$ with known normals $n_v$ at each vertex position $v$. Let $\cl N(v)$ be the $1$-ring neighborhood of $v$ (ie the set of vertices sharing an edge with $v$). Then  the filtered position of $v$ writes $v'=v+\delta v \vec n(v)$, where:

\begin{equation}\delta v=\frac{\sum_{p\in \cl N(v)} w_1(\|p-v\|)w_2(<n_v,p-v>)<n_v,p-v>}{\sum_{p\in \cl N(v)} w_1(\|p-v\|)w_2(<n_v,p-v>)|)}\label{eq:bilatfleishman}\end{equation}
In a nutshell, this means that the normal component of the  vertex $v$ is replaced by an weighted average of the normal components of its neighboring points which also are close to the plane tangent to the surface at $v$.  The distance to the tangent plane  takes for meshes the role that was taken for images by the distance between grey levels. If $v$ belongs to a sharp edge, then the only points close to the tangent plane at $v$ are the points on the edge. Thus,  the edge sharpness will not be smoothed away.  One of the drawbacks of the above filter is clearly the use of a mesh-dependent neighborhood. Yet, this is easily fixed by defining an intrinsic Euclidean neighborhood.
\medskip

Another adaptation of the 2D bilateral filter to surface meshes is introduced in \cite{non-iterative_feature_preserving}. This approach considers the bilateral filtering problem as a robust estimation problem for the vertex position. A set of surface predictors are linked to the mesh $\cl M$: for each triangle $f$ the position estimator $\Pi_f$ projects a point to the plane defined by $f$. Let $a_f$ be the surface area and $c_f$ be the center of $f$ then for each vertex $v$, the denoised vertex is

\begin{equation}v'=\frac1{C(v)}\sum_{f\in \cl M}\Pi_f(v)a_fw_1(\|c_f-v\|)w_2(\|\Pi_f(v)-v\|) \label{eq:bilatjones}\end{equation}
where $C(v)=\sum_{f\in \cl M}a_fw_1(\|c_f-v\|)w_2(\|\Pi_f(v)-v\|)$ is the weight normalizing factor and $w_1$ and $w_2$ are two Gaussians.

Thus, $w_1(||c_f-v||)$ is a weight which is small if the triangle $f$ is close to $v$. Thus, this term is the classic locality-in-space term of the bilateral. Similarly, $w_2(\|\Pi_f(v)-v\|)$ measures how far the projection of $v$ onto the plane of the triangle is from $v$. This weight favors the triangles $f$ whose  plane is coherent with $v$.

The projection on the tangent planes operator $\Pi_f$ depending on the normal to $f$, the normals must be robustly estimated. Normals, being first order derivatives, are more subject to noise than vertex positions. Hence the method starts by denoising the normal field. To do so, the mesh is first smoothed by using the same formula as above without the influence weight $w_2$ and with $\Pi_f(v)=c_f$, namely an updated position

$$v'=\frac1{C(v)}\sum_{f\in \cl M}c_fa_fw_1(\|c_f-v\|)$$
where $C(v)=\sum_{f\in \cl M}a_fw_1(\|c_f-v\|)$.
The normal for each face in the denoised mesh is then computed and assigned to the corresponding face of the original noisy mesh. It is with this robust normal field that the bilateral filter of equation  (\ref{eq:bilatjones}) is applied in a second step.
The bilateral filter is not iterated in this method.

In \cite{bilateral_recovering_sharp_edges}, meshes  feature-insensitively sampled  are denoised using a related bilateral approach. (By feature-insensitively we mean that the mesh sampling is independent of the features of the underlying surface, like (e.g.) a uniform sampling.)
The algorithm proceeds as follows: it detects the shape geometry (namely
sharp-regions), denoises the points and finally optimizes the mesh by removing
thin triangles. The bilateral filter is defined in a manner similar to
\cite{non-iterative_feature_preserving}. Let $v$ be a mesh vertex, $\cl N(v)$
the set of triangles within a given range of $v$ and $n_f$, $a_f$, $c_f$ the
normal, area and center of a facet  $f$ (a triangle). Denote by $\Pi_f(v)$ the
projection of $v$ onto the plane of $f$, then the denoised vertex is defined by

$$v'=\frac1{C(v)}\sum_{f\in \cl N(v)}\Pi_f(v)a_fw_1(\|c_f-v\|)w_2(\|\Pi_f(v)-v\|)$$

where $C(v)=\sum_{f\in \cl N(v)}a_fw_1(\|c_f-v\|)w_2(\|\Pi_f(v)-v\|)$ (weight normalizing factor).

The first step is to detect sharp regions. Several steps of bilateral filtering
(as defined in \cite{non-iterative_feature_preserving})
are applied, then a smoothness index is computed by measuring the infimum of angles between normals of faces adjacent to $v$. By thresholding  this measurement, the sharp vertices are selected. Triangles whose three vertices are sharp and whose size does not increase during the bilateral iterations are marked as sharp. This detection done, points are restored to their original positions. Then the bilateral filtering formula is applied to sharp vertices only, and the geometry sharpness is encoded into a Hermite data collection containing normals, centers and areas of filtered triangles. Points are then restored to their original position. Each sharp vertex is moved using the bilateral filtering over the neighboring stored data units, and thin vertices are removed from the mesh (these last two steps are iterated a certain number of times). Finally a post-filtering step consists in applying one step of bilateral filtering on all non sharp edges.
 \medskip

In \cite{fuzzyCmeans} a two-steps denoising method  combines the fuzzy C-means clustering method (see \cite{dunn_fcm}) with a bilateral filtering approach.
Fuzzy C-means is a clustering technique that allows a  piece of data to belong to two different clusters. Each point $p$ gets a parameter $\mu_{p,k}$ which measures the degree of membership of $p$ to a cluster $k$.
Let $m_p$ be the number of points in the spherical neighborhood $S$ of a point $p$. If $m_p<threshold$ the point is deleted. Otherwise a fuzzy C-means clustering center $c_p$ is associated with  $p$. The normal at point $c_p$ is computed as the normal to the regression plane of the data set in a spherical neighborhood of $p$.   Fleishman's bilateral filter \cite{bilateral} is used to filter $c_i$ which yields the denoised point.
This hybrid and complex method  is doubly bilateral. Indeed, the previous C-means clustering selects an adapted neighborhood for each data set point and replaces it by an average which is by itself the result of a first bilateral filter in the wide sense of neighborhood filter. Indeed, the used neighborhood for each point depends on the point.  The second part of the method therefore applies a second classical bilateral method to a cloud that has been filtered by a first bilateral filter.


\medskip

The bilateral filtering idea was also used as a part of a surface reconstruction process. In \cite{reconstruction_with_bilat}, for example, a method for reducing position and sampling noise in point cloud data while reconstructing the surface is proposed. A 3D geometric bilateral filter method for edge-preserving and data reduction is introduced.
 Starting from a point cloud, the points are classified in an octree, the voxel centers are filtered, representative surface points are defined and the mesh is finally reconstructed. A key point is that the denoising depends on the voxel decomposition. Indeed, the filter outputs a result for each voxel. For a voxel $V$, call $v$ its centroid with normal $n_v$. Let $w_1$ and $u_2$ be two functions weighting respectively the distance between a point $p$ position and the centroid location and $\delta(p,v)=\langle n_p,n_v\rangle$ be the scalar product of the normal at $p$ normal and the normal at the  centroid.
 Then the output of the filter for voxel $V$ is
$$v'=\frac1{C(v)}\sum_{p\in V}w_1(\|p-v\|)u_2(\delta(p,v))p$$
where $C(v)=\sum_{p\in V}w_1(\|p-v\|)u_2(\delta(p,v))$.
Here $w_1$  is typically a Gaussian and $u_2$ is an increasing  function on $[0,1]$.
 But this filter proves unable to recover sharp edges, so a modification is introduced: prior to any filtering for each voxel $V$, points of $V$ are projected onto a sphere centered at the centroid $v$. Each mapped point is given a normal $\tilde n_p$ which has direction $p-v$ and is normalized. The geometric filtering is reduced to:
$$v'=\frac1{C(v)}\sum_{p\in V}u_2(\delta(\tilde n_p,n_v))p \text{ with } C(v)=\sum_{p\in V}u_2(\delta(\tilde n_p,n_v)).$$
 Although only the similarity of normals  is taken into account  in the above formula, the filter is bilateral because the average is localized in the  voxel.

\medskip
In \cite{bilateral_new_predictor}, the bilateral filter is interpreted as the association to each vertex $v$  of a weighted average
$$v'=\frac{\sum_{p\in \cl N(v)}w_1(\|p-v\|)w_2(\|\Pi_p(v)-v\|)\Pi_p(v)}{\sum_{p\in \cl N(v)}w_1(\|p-v\|)w_2(\|\Pi_p(v)-v\|)}$$
where  $\Pi_p(v)$ a \emph{predictor} which defines a ``denoised position of $v$ due to $p$'', namely the {\it projection of $v$ on the plane passing by $p$ and having the normal $n_v$}. The used bilateral predictor coming from \cite{bilateral} is $\Pi_p(v)=v+((p-v)\cdot  n_p) n_p$. In \cite{non-iterative_feature_preserving}, the used predictor was $\Pi_p(v)=v+((p-v)\cdot n_p) n_p$ which is the projection of $v$ on the tangent plane passing by $p$. With this last predictor the corners are less smoothed out,  yet there is a tangential drift due to the fact that the motion is not in the normal direction $n_v$ but in a averaged directions of the $n_p$ for $p\in \cl N(v)$. Therefore a new predictor is introduced:

$$\Pi_p(v)=v+\frac{(p-v)\cdot  n_p}{n_v \cdot  n_p}n_v$$
 This predictor tends to preserve better the edges than all other bilateral filters.
\medskip

The question of choosing automatically the parameters for the bilateral filter was raised in \cite{adaptive_bilateral}. It was proposed to choose adaptive parameters. The adaptive bilateral normal smoothing procedure process starts by searching for the set of triangles $(T_i)_i$ whose barycenters are within a given distance of a center triangle $T$. (But this keeps a distance parameter anyway.) Then the influence weight parameter $\sigma_s$ is computed as the standard deviation of the distance between normals $\|n(T_i)-n(T)\|$. The spatial weight parameter is estimated using a minimum length descriptor criterion (for various scales). The estimated parameters are then used to get the smoothed normal. This result is then used for rebuilding the mesh using the smoothed normals using the method described in \cite{adaptive_anisotrop_gauss_filter}.

\medskip
The bilateral filtering of meshes proving an efficient denoising  method preserving sharp features, the idea of a trilateral filtering is a natural extension taking into account still more geometric information.
For example \cite{trilateral_filter} proposes an extension to oriented meshes of the trilateral image filter. It is a $2$-pass filter: a first pass filters the normals and a second pass filters the vertex positions. Starting from an oriented mesh, a first pass denoised bilaterally the vertices normals using the following update

$$n'_v=\frac1{C(n_v)}\sum_{p\in\cl N(v)} n_p w_1(\|p-v\|)w_2(\|n_p-n_v\|)$$
where $C(n_v)=\sum_{p\in\cl N(v)} w_1(\|p-v\|)w_2(\|n_p-n_v\|)$.
Then, an adaptive neighborhood $\cl N(v)$ is found by iteratively adding faces near $v$ until the normals $n_f$ of face $f$ differ too much from $n'_v$. A function measuring the similarity between normals is built using a given threshold $R$,

$$f(v,f)=1 \text{ if $\|n'_v-n_f\|<R$; } 0 \text{ otherwise}.$$

The trilateral filter for normals filters a difference between normals. Define $n_\Delta(f)=n_f-n'_v$. Then the trilaterally filtered normal $n_v$ is

$$n''_v=n'_v+\frac1{C(v)}\sum_{f\in \cl N(v)} n_\Delta(f)w_1(\|c_f-v\|)w_2(n_\Delta(f))f(v,f)$$
where $C(v)=\sum_{f\in \cl N(v)} w_1(\|c_f-v\|)w_2(n_\Delta(f))f(v,f)$.
Finally, the same trilateral filter can be applied to vertices. Call $P_v$ the plane passing through $v$ and orthogonal to $n'_v$. Call $\tilde c_f$ the projection of $c_f$ onto $P_v$ and $c_\Delta(f)=\|\tilde c_f-c_f\|$. Then the trilateral filter for vertices, using the trilaterally filtered normal $n''_v$ writes

$$v'=v+n''_v \frac1{C(v)}\sum_{p\in\cl N(v)}c_\Delta(f) w_1(\|\tilde c_f-v\|)w_2(n_\Delta(f))f(v,f)$$
where $C(v)=\sum_{p\in\cl N(v)}w_1(\|\tilde c_f-v\|)w_2(c_\Delta(f))f(v,f)$
 % Starting from a point $V$, the algorithm finds a set $\cl N_F$ of connected faces with normals similar to the normal at the point \begin{Jean-Michel}
% How do you define the normal to the  point?\end{Jean-Michel}
% \begin{julie}
%  it is supposed an oriented mesh, I added it
% \end{julie}
%
%  (``similar'' means that the angle between the point normal and the face normal is less than a given threshold). It also allows to define a filtered plane passing through the center point and orthogonal to the point's filtered normal. Denote by $V$ the center point with normal $N_V$, $F$ a face, $X_F$ the face center, and $N_F$ the face normal. Set $f_\theta(V,F)=1$ if $N_v\cdot N_F<R $ and $0$ otherwise ($R$ is a given threshold), let $s_N$ and $c_N$ be two Gaussians.  Then the filtered normal is
%
% $$N_{v,out}=\frac1C\sum_{F\in \cl N_F}N_Fc_N(\|V-X_F\|)s_N(\|N_V-N_F\|)f_\theta(V,F)$$
%
% where $C=\sum_{F\in \cl N_F}c_N(\|V-X_F\|)s_N(\|N_V-N_F\|)f_\theta(V,F)$
% \begin{Jean-Michel}  The formula you wrote was
% $$N_{v,out}=N_v+\frac1C\sum_{F\in \cl N_F}(N_V-N_F)c_N(\|V-X_F\|)s_N(\|N_V-N_F\|)f_\theta(V,F).$$
% I tried to adjust it because it seemed incoherent. Indeed, it is not an average of normals. It  is a normal ($N_v$) + an average of differences of normals. But the differences  of normals seem to me to have the wrong sign because we have basically $N_v+averageof(N_v-N_f)$. Please, confirm or correct.
% \end{Jean-Michel}
% \begin{julie}
%  No the first expression had the right form. I will rewrite this whole part, because the article itself needs to be made clearer.
% \end{julie}
%
%
% Once all normals are smoothed point positions are filtered. Call $p_F$ the projection of $X_F$ onto the filtered plane,  $c_V$ and $s_V$ two gaussian functions, then the updated position $V$ writes
%
% $$X_{v,out}=X_v+N_{v,out}\frac1C\sum_{F\in \cl N_F}\|X_F-p_F\|c_V(\|X_V-p_F\|)s_V(\|X_F-p_F\|)f_\theta(V,F)$$
%
% where $C=\sum_{F\in \cl N_F}c_V(\|X_V-p_F\|)s_V(\|X_F-p_F\|)f_\theta(V,F)$
%
% \begin{Jean-Michel} Here again I was unable to  understand. Indeed, the formula adds to $X_v$ a term which has the dimension and direction of $N_v$. Thus it  all depends on the direction of $N_v$: is $N_v$ pointing towards the concavity or convexity of the  surface?
% I tried to make a sense of these formulas but failed. Depending on the orientation of the normals to the facets, these definitions can have completely different results!
% \end{Jean-Michel}
The results are similar to \cite{non-iterative_feature_preserving} though slightly better. They are comparable to the results of \cite{bilateral} since both methods use the distance to the tangent plane as a similarity between points.
\medskip

In \cite{similarityPSS} the authors proposed a trilateral filter with slightly different principles.
A {\it geometric intensity} of each sampled point is first defined as depending of the neighborhood of the point
$$\delta(p)=\frac{\sum_{q\in\cl N(p)}w_{pq}<n_p,q-p>}{\sum_{q\in\cl N(p)}w_{pq}}$$
with $$w_{pq}=w_1(\|q-p\|)w_2(\|<n_p,q-p>\|)w_h(\|H_{q}-H_p\|).$$

This type of filter is a trilateral filter, which means that it depends on three variables: distance between the point $p$ and its neighbors $q$, distance along the normal $n_p$ between the point $p$ and its neighbors $q$, and the difference of their mean curvatures $H_p$ and $H_q$)

At each point, a local grid is built on the local tangent plane (obtained by local covariance analysis), at each point of this grid, the geometry intensity is defined by interpolation. Thus, neighborhoods of same geometry are defined for each pair of distinct points and the similarity can be computed as a decreasing function of the $L^2$ distance between these neighborhoods.

Since the goal is to denoise one point with  similar points, the algorithm proposes to cluster the points into various classes by the mean shift algorithm. To denoise a point, only points of the same class are used. This gives a denoised geometry intensity $\delta'_i$ and the final denoised position $p'=p+\delta(p) n_p$.
\medskip

More recently the NL-means (\cite{nonlocalmeans}) method which proved very powerful in image denoising was adapted to meshes and point clouds in
  \cite{meshNLmeans}. Recall that for an image $I(x)$, the NL-means filter computes a filtered value $J(x)$ as
$J(x)=\frac1{C(x)}\int_\Omega w(x,y)I(y)dy$, an adaptive average  with  weights $w(x,y)=\exp{-\frac1{h^2}\int G_a(|t|)|I(x-t)-I(y-t)|^2dt}$ and $C(x)=\int_\Omega w(x,y)dy$.
Here $G_a$ is a gaussian or a compactly supported function, so that it defines a patch.
Thus, the denoised point is a mean of pixel values with weights measuring the local image similarity of patches around other pixels with the patch around the current pixel.

Consider now  the adaptation to a mesh $\cl M$. Let $\Omega_\sigma(x)=\{y\in\cl M |\quad |x-y|\leq 2\sigma\}$. The smoothing is done by changing $x$ at each step: $x^{n+1}=x^n+k(x^n)n_x^n$ with $n_x$ the normal to $\cl M$ at $x$. Let $S_y$ be the surface associated to node $y$. The following definitions are directly adapted from the image case (we adopt here a continous formalism for clarity):
$$k(x)=\frac1{C(x)}\int_{\Omega_{\sigma_2}}w(x,y)I(y)dS_y$$
$$C(x)=\int_{\Omega_{\sigma_2}}w(x,y)dS_y$$
$$I(y)=<n_x,y-x>$$
$$w(x,y)=\exp{-\frac{D(x,y)}{h^2}}$$

The problem is to define the similarity kernel $D$. Let $\sigma_2>\sigma_3$ be the half radius of the comparison domain: $\sigma_3$ is the half radius of the neighborhood to define the geometric similarity between two points. $\sigma_2$ is the half radius of the domain where similar points are looked for.
The local tangent plane at $y$ is parameterized by $t_1$ et $t_2$. For all $z$ of $\Omega_{\sigma_2}(y)$ the translation $t$ is defined as $t=-(<t_1,z-y>,<t_2,z-y>)$.

A local approximation by RBF functions is built:
$F_x(u,v)=p(u,v)+\sum_{w\in\Omega_{\sigma_3}(x)}\lambda_w\psi(\sqrt{u^2+v^2})$.

RBF coefficients $\lambda_\omega$ are found by a linear system resolution. And the similarity kernel finally yields:
$$D(x,y)=\int_{\Omega_{\sigma_3}(y)}G_{\sigma_3}(|t|)|F_x(u_z,v_z)-I(y-t)|^2dt$$ with $I(x-t)=<n_x,z-x>$ et $G_\sigma$ a gaussian kernel.


Thus each node is compared with nodes in a limited domain around it and the weighted mean over all these nodes yields the denoised position. This results in a better feature preserving mesh denoising method, but at the cost of a considerably higher computation time.
\medskip

To improve the computation time when denoising datas using neighborhood filters, Gaussian kd-trees were introduced in \cite{kdtree_filtering}. The method proposed was designed to compute efficiently a class of $n$-dimensional filers which replace a value by a linear combination of other values. The basic idea is to consider those filters as nearest neighbors search in a higher dimensional space, for example $(r,g,b,x,y)$ in case of a 2D color image and of a bilateral filter. To accelerate this neighbor search, a Gaussian k-d tree is introduced.
 The non local means filter in its naive implementation has a $O(n^2f^2)$ complexity. To apply Gaussian k-d tree, the position of a pixel is set to be the patch and the value is set to be the color value of the pixel. A simple PCA on patches helps to capture the dimensions that best describe the patches.
The authors also applied the Gaussian k-d tree to perform 3D NL-means on meshes or point clouds. To produce a meaningful value to describe geometry, they use the idea of spin images; At each point sample, a regression plane can be estimated and the coordinates of the point in this coordinate system are used as values; the position of the patch is set to be a couple of values (namely the spin image which is a histogram of cylindrical coordinates).
\medskip

\paragraph{Summary of 3D mesh bilateral filters}
The filters we reviewed in this section are almost all defined for meshes. Yet, with very little effort almost all of them can be adapted to unstructured point clouds by simply  redefining the neighborhoods as the set of points within a given distance from the center point (spherical neighborhood).  We have examined several classic variants of bilateral filters, but their main principle is to perform an average of neighboring vertices pondered by the distance of these vertices to an estimated tangent plane of the current vertex. This distance takes the role played by the grey level in image  bilateral filters. This distance can be implemented in several ways by either projecting the current vertex to the neighboring triangles, or by projecting the neighboring vertices on the current triangle, or by using an estimate of the  normal at the current vertex which has been itself previously filtered. An interesting and simple possibility is to directly combine distance  of vertices and of their normals or even distances of vertices, normals, and curvatures (but this requires a previous smoothing to get denoised normals and curvatures). Notice that position, normal and curvature characterize the shape of cloud in a larger neighborhood. Thus at this point the obvious generalization of bilateral filters is NL-means, which directly compares point-wise the  shape of the neighborhood of a vertex with the overall shape of the neighborhoods of others before performing an  average of the most similar neighborhoods to deliver a filtered neighborhood.

Sticking to the simplicity of comparisons and to the essentials of bilateral filter, we shall be contented in the comparative section to illustrate the  gains of the bilateral filter with respect to a (good) implementation of its unilateral counterpart, the mean curvature motion, performed by projection of each vertex on a local regression plane.
The remainder of this section is divided as follows: section \ref{sec:artificial} presents experiments and comparisons on artificial shapes, and section \ref{sec:real} present results on some real shapes.


\subsection{Experiments on artificial shapes \label{sec:artificial}}

In the following experiments, we will compare the denoising of the bilateral filter as introduced in \cite{bilateral} with the mean curvature motion. Recall that \cite{bilateral} defined the update of a point as:

$$\delta v=\frac{\sum_{p\in \cl N(v)} w_1(\|p-v\|)w_2(<n_v,p-v>)<n_v,p-v>}{\sum_{p\in \cl N(v)} w_1(\|p-v\|)w_2(<n_v,p-v>)|)}$$
(see first section, equation \ref{eq:bilatfleishman} for notations).
The mean curvature motion used here is the projection on the regression plane: a vertex $v$ with normal $n_v$ and spherical neighborhood $\cl N(v)$ is projected on the regression plane of $\cl N(v)$. In \cite{preprint1}, this operator was shown to be an approximation of the mean curvature motion:
$$\pd v t=H n_v$$
Let us first observe the effects of bilateral denoising on some artificial shapes. A $5$ unit side-long cube is  created with added gaussian noise with standard deviation $0.02$ (fig \ref{fig:noisy_cube}).
\begin{figure}
 \begin{center}
   \includegraphics[scale=0.5]{cube_noisy.eps}
 \end{center}
 \vspace{-8mm}
\caption{A noisy cube with gaussian noise \label{fig:noisy_cube}}
\end{figure}


%\begin{figure}[H]
% \begin{center}
%  \subfigure[$1$ bilateral iteration]{\includegraphics[scale=0.5]{cube_b1.eps}}
%  \subfigure[$2$ bilateral iterations]{\includegraphics[scale=0.5]{cube_b2.eps}}
%  \subfigure[$5$ bilateral iterations]{\includegraphics[scale=0.5]{cube_b5.eps}}
% \end{center}
%\caption{Bilateral iterations. Notice how the sharpness is preserved \label{fig:cube_bilat}. The bilateral variant uses the distance between points and the distance to the point's regression plane as weights.}
%\end{figure}
%
%The experiments of fig \ref{fig:cubed} show the denoising power of the bilateral filter in term of preserving edges and should be compared with the standard mean curvature motion filter (fig \ref{fig:cubed}).
%
%%\begin{figure}[H]
% \begin{center}
%  \subfigure[$1$ iteration]{\includegraphics[scale=0.5]{cube_p1.eps}}
%  \subfigure[$2$ iterations]{\includegraphics[scale=0.5]{cube_p2.eps}}
%  \subfigure[$5$ iterations]{\includegraphics[scale=0.5]{cube_p5.eps}}
% \end{center}
%\caption{Mean curvature motion iterations. Notice how the sharpness is smoothed away \label{fig:cube_mcm}}
%\end{figure}

\begin{figure}[H]
 \begin{center}
  \subfigure[$1$ bilateral iteration]{\includegraphics[scale=0.4]{cubed_b1.eps}}
  \subfigure[$2$ bilateral iterations]{\includegraphics[scale=0.4]{cubed_b2.eps}}
  \subfigure[$5$ bilateral iterations]{\includegraphics[scale=0.4]{cubed_b5.eps}}
  \subfigure[$1$ MCM iteration]{\includegraphics[scale=0.4]{cubed_p1.eps}}
  \subfigure[$2$ MCM iterations]{\includegraphics[scale=0.4]{cubed_p2.eps}}
  \subfigure[$5$ MCM iterations]{\includegraphics[scale=0.4]{cubed_p5.eps}}
 \end{center}
 \vspace{-8mm}
\caption{Bilateral and MCM iterations on the cube corner. Notice how the sharpness is much better preserved by the bilateral filter than by the mean curvature equation. \label{fig:cubed} }
\end{figure}
The experiments of fig \ref{fig:cubed}(a, b, c) show the denoising power of the bilateral filter in term of preserving edges and should be compared with the standard mean curvature motion filter (fig \ref{fig:cubed} d, e f), the comparison is particularly interesting in the corner areas.
The bilateral filter implies an anisotropic curvature motion leading to a diffusion only in smooth parts while preserving the sharp areas.
Let us now see how those filters perform in case of a sharp edge. An estimation of the noise for each of the denoising methods is shown on tabular \ref{tab:sharp_noise}. These numbers tend to prove that Mean Curvature Motion, although it smoothes well the noisy flat parts also smoothes away the sharpness whereas the bilateral filter tends to preserve the sharp edges better. With few iterations the noisy parts are smoothed out, which decreases the root mean square error, then when iterating the operator, the sharpness tends to be smoothed, which increases the RMSE again. This phenomenon is of course far quicker with the mean curvature motion since this filter does not preserve the edge at all.

%\begin{figure}[H]
% \begin{center}
%\includegraphics[scale=0.5]{sharp_noisy.eps}
% \end{center}
%\caption{A noisy sharp edge with added gaussian noise \label{fig:sharp_noisy}}
%\end{figure}
%
%\begin{figure}[H]
% \begin{center}
%  \subfigure[$1$ iteration]{\includegraphics[scale=0.4]{sharp_b1.eps}}
%  \subfigure[$2$ iterations]{\includegraphics[scale=0.4]{sharp_b2.eps}}
%  \subfigure[$5$ iterations]{\includegraphics[scale=0.4]{sharp_b5.eps}}
% \end{center}
%\caption{Bilateral filter iterations. Notice how sharpness is preserved \label{fig:sharp_bilat}}
%\end{figure}
%
%\begin{figure}[H]
% \begin{center}
%  \subfigure[$1$ iteration]{\includegraphics[scale=0.4]{sharp_p1.eps}}
%  \subfigure[$2$ iterations]{\includegraphics[scale=0.4]{sharp_p2.eps}}
%  \subfigure[$5$ iterations]{\includegraphics[scale=0.4]{sharp_p5.eps}}
% \end{center}
%\caption{Mean curvature motion iterations. Notice how the sharpness is smoothed away \label{fig:sharp_mcm}}
%\end{figure}

\begin{figure}[H]
\begin{center}
\begin{tabular}{|c|c|c|c|c|}
\hline
& Input & Iteration 1 & Iteration 2 & Iteration 5 \\ \hline
RMSE (bilateral)& $0.01$ & $0.0031$ & $0.0019$ & $0.0035$\\ \hline
RMSE (mcm)& $0.01$ & $0.0051$ & $0.0085$ & $0.0164$\\ \hline
\end{tabular}
\end{center}
\vspace{-8mm}
\caption{Noise estimation for the Sharp edge denoising \label{tab:sharp_noise}}
\end{figure}

\subsection{Experiments on real shapes \label{sec:real}}
This section runs some experiments on a real shape, namely Michelangelo's David.
The experiments were run on the point cloud. At each step an interpolating mesh
was built for visualization.
\begin{figure}[H]
 \begin{center}
 \subfigure[Initial noisy David]{\includegraphics[scale=0.25]{david_noisy_back.eps}}
 \subfigure[Bilateral denoising]{\includegraphics[scale=0.25]{david_noisy_back_bilat1.eps}}
 \subfigure[MCM]{\includegraphics[scale=0.25]{david_noisy_back_mcm1.eps}}
 \end{center}
 \vspace{-1cm}
\caption{Denoising of the David (back) \label{fig:david_back_denoising}}
\end{figure}

\begin{figure}[H]
 \begin{center}
 \subfigure[Original of the David] {\includegraphics[scale=0.24]{david.eps}}
  \subfigure[Noisy David]{\includegraphics[scale=0.2]{david_noisy2.eps}}\\
  \subfigure[Bilateral denoising]{\includegraphics[scale=0.2]{david_noisy2_bilat.eps}}
  \subfigure[MCM]{\includegraphics[scale=0.2]{david_noisy2_mcm.eps}}
 \end{center}
 \vspace{-1cm}
\caption{Denoising of the David's face \label{fig:david_denoising}}
\end{figure}

\begin{figure}[H]
 \begin{center}
 \subfigure[Iteration $1$]{\includegraphics[scale=0.25]{david_noisy_back_bilat1.eps}}
 \subfigure[Iteration $2$]{\includegraphics[scale=0.25]{david_noisy_back_bilat2.eps}}
 \end{center}
 \vspace{-8mm}
\caption{Iterating the bilateral filter on the David (back) \label{fig:david_back_denoising_iter}}
\end{figure}

\begin{figure}[H]
 \begin{center}
 \subfigure[Bilateral Filtering]{\includegraphics[scale=0.15]{david_detail_top_bilat.eps}}
 \subfigure[MCM]{\includegraphics[scale=0.15]{david_detail_top_mcm.eps}}
 \end{center}
 \vspace{-8mm}
\caption{Detail of the David \label{fig:david_hair_denoising}}
\end{figure}


\begin{figure}[H]
 \begin{center}
 \subfigure[Initial scan]{\includegraphics[scale=0.15]{screw_noisy.eps}}
 \subfigure[Bilateral Denoising]{\includegraphics[scale=0.15]{screw_bilat.eps}}
 \subfigure[MCM]{\includegraphics[scale=0.15]{screw_mcm.eps}}
 \end{center}
 \vspace{-8mm}
\caption{Denoising of a screw nut driver scan.   \label{fig:screw}}
\end{figure}



\begin{figure}[H]
 \begin{center}
 \subfigure[Initial fragment]{\includegraphics[scale=0.4]{31u_input.eps}}\\
 \subfigure[Bilateral Denoising]{\includegraphics[scale=0.4]{31u_bilat.eps}}
 \subfigure[MCM]{\includegraphics[scale=0.4]{31u_mcm.eps}}
 \end{center}
 \vspace{-8mm}
\caption{Denoising of fragment '31u' of Stanford Forma Urbis Romae \label{fig:31u}}
\end{figure}
On fig \ref{fig:david_denoising}, denoising
artefacts created by the bilateral filter can be seen. They appear as
oscillations (on David's cheek for example). These artifacts can be explained by
the fact that the bilateral filter enhances structures.  Added noise structures
can be randomly enhanced by the bilateral filter.
On fig \ref{fig:david_back_denoising}, we can see that some noise remain after one iteration of bilateral denoising. We therefore iterate the bilateral filter with same parameters. Then, obviously, the remaining noise disappears at the cost of some sharpness loss (see fig \ref{fig:david_back_denoising_iter}). This can also be seen on a fragment of the Stanford Forma Urbis Romae Project fig \ref{fig:31u} and on a noisy simple scan of a screw nut driver.


\newpage
\section{Depth-oriented applications}
\label{sec:relief}

This section focuses on the applications of the bilateral filter and its generalized version to depth-oriented image processing tasks.
The common idea to all these applications is of constraining the diffusion of depth information to the intensity similarity between pixels.
The  underlying assumption is that pixels with similar intensity around a region are likely to have similar depths. Therefore, when diffusing depth information  basing on intensity similarity, the discontinuities in depth are assured to be consistent with the color discontinuities, which is often a desirable property \cite{Gamble:87,kellman:91}.


The remainder of this section is organized as follows.
Section \ref{subsec:search_correspondence} reviews the applications of the bilateral filter to stereo matching algorithms, while section \ref{subsec:range_images} describes an application to the resolution enhancement of range images.
Section \ref{subsec:monocular_depth_by_nonlinear_diffusion}  review applications to the estimation of depth in single images.








\subsection{Bilateral filter for improving the depth map provided by stereo matching algorithms}
\label{subsec:search_correspondence}


Stereo matching algorithms address the problem of recovering the depth map of 3D scene from two images captured from different viewpoints. This is achieved by finding a set of points in one image which can be identified as the same points in the other image. In fact, the knowledge of the point-to-point correspondences permits to compute the relative disparities, which are directly related to the distance of the object normal to the image plane.
The search of correspondences is usually based on a similarity measure between pixels belonging to two different images of the same scene.
Due to the presence of noise and repetitive texture, these correspondences are extremely difficult to find without global reasoning.
In addition, occluded and textureless regions are inherently ambiguous since the relative disparity cannot be found reliably by image matching alone.

In applications demanding  accurate disparities such as structure reconstruction or multiview video coding, stereo matching algorithms provide a sparse disparity map.
Improvements can be obtained through filtering or interpolation, for instance by using median or morphological filters. However, their ability to do so is limited. In \cite{Yin04}, the authors  have proposed a post-processing step to improve dense depth maps produced by any stereo matching algorithm. The proposed method consists in applying an iterated bilateral filter, which diffuses the depth values by relying on the gradient of the original image instead of that of the depth image. This allows to incorporate edge information into the depth map, assuring discontinuities in depth to be consistent with intensities discontinuities.

In \cite{Yoon:06}, the color-weighted correlation idea underlying the bilateral filter has been exploited to reduce the ambiguity of the  correspondence search problem.  Classically, this problem has been addressed by area-based methods relying on the use of local support windows.
In this approach,  all pixels in a support window are assumed to have similar depth in the scene and, therefore, similar
disparities. Accordingly, pixels in homogeneous regions get assigned the disparities inferred from the disparities of neighboring pixels.
However, when the support windows are located in correspondence of depth discontinuities, the same disparity is assigned to pixels having different depths, resulting in a foreground-fattening phenomenon.
To obtain accurate results, an appropriate support window should be selected for each pixel adaptively.
This problem is addressed in \cite{Yoon:06} by weighting the pixels in a given support window taking into account their color similarity and geometric proximity to the reference pixel.
The similarity measure between two pixels is then computed using the support-weights in the support windows of both pixels, allowing to incorporate the edge information into the disparity map.
Experimental results show that the use of adaptive support-weights produces noise-removed accurate piecewise smooth disparity maps while preserving depth discontinuities.

The idea of exploiting the color-weighted correlation to reduce the ambiguity of the correspondence problem has been implemented in a parallel architecture \cite{Yang:06} \cite{wang:06}, allowing its use in  real-time applications and integrated in more complex stereo systems \cite{Yang:06} \cite{wang:06}, which nowadays achieve the top rank in  the Middlebury benchmark \cite{Scharstein}.




The interpolation of disparity maps and in particular of Digital Elevation Models (DEMs) has been considered in several recent works.  \cite{Facciolo06} proposes to interpolate unknown areas by constraining a diffusion anisotropic process to the geometry imposed by a reference image, and coupling the process with a data fitting term which tries to adjust the reconstructed surface to the known data. More recently, \cite{Facciolo09} has proposed a new interpolation method which defines a geodesic neighborhood and fits an affine model at each point. The geodesic distance is used to find the set of points that are used to interpolate a piecewise affine model in the current sample. This interpolation is refined by merging the obtained affine patches with a Mumford-Shah like algorithm.  The \textit{a contrario} methodology has been used in a merging procedure.  In the urban context, \cite{Lafarge08} uses a dictionary of complex building models to fit the disparity map. However, the applicability of such a method is less evident because of the initial delineation of buildings by a rectangle fitting.



The bilateral filter averages the pixel colors, based on both their geometric closeness and their photometric similarity,  preferring of course near values to distant values in space and color.    \cite{Yoon06} have used it to weight the correlation windows before the stereo correspondence search.
 We shall illustrate this process with experiments from  \cite{Neusthesis}. Let $q$ be a point in $I$. Consider $L_q\subset I$ the subimage where the weight is learned. For each $p \in L_q$ the weight due to color similarity and proximity are computed:


\begin{description}
\item[color similarity:] We consider the color distance $$d_c(u_q,u_p) = \left((R_u(q) -  R_u(p))^{2} + (G_u(q) -  G_u(p))^{2} + (B_u(q) -  B_u(p))^{2}\right)^{1/2}\,,$$ where $R_u$,$G_u$ and $B_u$ are the red, green and blue channels of $u$. Then the weight corresponding to the color similarity between $p$ and $q$ is


$$w_c(p,q) = exp \left(-\dfrac{{d_c(u_q,u_p)}^2}{h_1^{2}}\right)\,.$$

\item[proximity:] We consider the Euclidean distance between the points positions in the image plane $$d(q,p)=\left((q_1-p_1)^{2} + (q_2-p_2)^{2}\right)^{1/2}\,,$$  where $p=(p_1,p_2)$ and $q=(q_1,q_2)$. Then the weight corresponding to proximity is


$$w_d(p,q) = exp \left(-\dfrac{d(q,p)^{2}}{h_2^{2}}\right)\,.$$

 \end{description}


Therefore, the total associated weight between the two points $q$ and $p$ is

$$
W(p,q) =  \dfrac{1}{Z_q} \, w_c(p,q)w_d(p,q) =
 \dfrac{1}{Z_q} \,  exp \left(-\left(\dfrac{d_c(u_q,u_p)^{2}}{h_1^{2}}+ \dfrac{d(q,p)^{2}}{h_2^{2}} \right)\right)\,,
$$
where $Z_q$ is the normalizing factor $Z_q =\displaystyle { \sum_{p \in L_q} w_c(p,q)w_d(p,q)}\,.$
The interpolated disparity map $\mu_I$ is computed via an iterative schema


$$
\mu_I(q,k)=\sum_{p \in L_q} W(p,q) \mu_I(p,k-1)\,,
$$
where $k$ is the current iteration and the initialization $\mu_I(\cdot,0)=\mu_M(\cdot)$.




Figures \ref{interpolated_tsukuba_venus} and \ref{interpolated_teddy_cones} show the interpolated Middlebury results (100\% density).  The experiments demonstrate that, starting from a disparity map which is very sparse near image boundaries, the bilateral diffusion process can recover a reasonable depth map.

\begin{figure}
\begin{center}
 \begin{tabular}{ccc}
\includegraphics[width=4.5cm]{color_imL_tsukuba.eps}  &
\includegraphics[width=4.5cm]{color_imR_tsukuba.eps}  &
\includegraphics[width=4.5cm]{tsukuba_finaldisp.eps} \\
\includegraphics[width=4.5cm]{interpol_tsukuba.eps}   &
\includegraphics[width=4.5cm]{groundtruth_tsukuba.eps} \\
%\includegraphics[width=4.5cm]{error_tsukuba.eps}\\

\includegraphics[width=4.5cm]{color_imL_venus.eps}&
\includegraphics[width=4.5cm]{color_imR_venus.eps} &
\includegraphics[width=4.5cm]{venus_finaldisp.eps}\\
\includegraphics[width=4.5cm]{interpol_venus.eps}&
\includegraphics[width=4.5cm]{groundtruth_venus.eps}\\
%\includegraphics[width=4.5cm]{error_venus.eps}\\

\end{tabular}
\end{center}
\caption{Tsukuba and Venus results. For each couple of images: stereo pair of images, output of a sparse algorithm retaining only sure points, points (in red in the .pdf file) are the rejected correspondences, interpolated version of these results, ground truth and signed error between the interpolated image and the ground truth.}
 \label{interpolated_tsukuba_venus}
\end{figure}



\begin{figure}
\begin{center}
 \begin{tabular}{ccc}

\includegraphics[width=4.5cm]{color_imL_teddy.eps}&
\includegraphics[width=4.5cm]{color_imR_teddy.eps}&
\includegraphics[width=4.5cm]{teddy_finalcolor.eps}\\
\includegraphics[width=4.5cm]{interpol_teddy.eps}&
\includegraphics[width=4.5cm]{groundtruth_teddy.eps}\\
%\includegraphics[width=4.5cm]{error_teddy.eps}\\

\includegraphics[width=4.5cm]{color_imL_cones.eps}&
\includegraphics[width=4.5cm]{color_imR_cones.eps}&
\includegraphics[width=4.5cm]{cones_finalcolor.eps}\\
\includegraphics[width=4.5cm]{interpol_cones.eps}&
\includegraphics[width=4.5cm]{groundtruth_cones.eps}
%\includegraphics[width=4.5cm]{error_cones.eps}\\
\end{tabular}
\end{center}
\caption{Teddy and Cones results. For each couple of images: stereo pair of images, output of a sparse algorithm retaining only sure points, points (in red in the .pdf file) are the rejected correspondences, interpolated version of these results, ground truth and signed error between the interpolated image and the ground truth.}
 \label{interpolated_teddy_cones}
\end{figure}











\subsection{Bilateral filter for enhancing the resolution of low-quality range images}
\label{subsec:range_images}



In \cite{yang:07}, the authors have proposed a post-processing step to enhance the resolution of low-quality range images.
Contrary to intensity images, each pixel of a range image expresses the distance between a known reference frame and a visible point in the scene.
Range images are acquired by range sensors that, when acquired at video rate, are either very expensive or very limited in terms of resolution.
To increase the resolution of low-quality range images acquired at video rate, the authors of \cite{yang:07} have proposed a post-processing step relying on an iterated bilateral  filter.
The filter diffuses the depth values of the low-quality range image, steering the diffusion by the color information provided by a registered high-quality camera image.
The input range image $z_0$ is first up-sampled to the same size of the camera image and then used as the initial depth map hypothesis according to which an initial cost volume is built, $c_0$.  An iterative process is then applied that consists in first performing a bilateral filter on each slide of the current cost volume $c_i$ and then in generating a new depth map hypothesis $z_{i+1}$ based on the filtered cost volume. At each filtering step, an aggregated  cost volume is computed for each pixel based on the cost volume of the pixels in its neighborhood, which is defined  based on the color similarity.
The new depth map hypothesis $z_{i+1}$ is then obtained by first selecting the depth hypothesis having the minimal cost volume and then by applying a sub-pixel estimation using a quadratic polynomial interpolation.



\subsection{Bilateral filter for the global integration of local depth information}
\label{subsec:monocular_depth_by_nonlinear_diffusion}

In \cite{dimiccoli:08}, an iterated neighborhood filter has been proposed as a depth diffusion mechanism to globally integrate a set of initial local depth hypothesis derived from different monocular depth cues.




\begin{figure}
\begin{center}
\begin{tabular}{ccrrrrr}
 \includegraphics[width=0.45\linewidth]{diff2.eps}&
  \includegraphics[width=0.45\linewidth]{diff3.eps}
  \\(a) & (b)
  \\

  \includegraphics[width=0.45\linewidth]{diff4.eps}&
    \includegraphics[width=0.45\linewidth]{diff5.eps}
  \\ (c) & (d)

\end{tabular}
\end{center}


\caption{Example of depth diffusion by using equation (\ref{eq:depth_by_average}). (a) Gray level image, where BSPs and FSPs are marked in blue and red respectively. (b) Depth image, where points corresponding to FSPs are initialized with a positive value (marked in white) and the rest of the image with value zero. (c) and (d) Depth images after an increasing number of iterations of the DDF. }

\label{fig:example_diffusion}
\end{figure}


\begin{figure}
\begin{center}
\begin{tabular}{ccrrrrrrrrrrr}
 \includegraphics[width=0.30\linewidth]{diff33.eps}&
  \includegraphics[width=0.30\linewidth]{diff6.eps}
  \\ (a) & (b)
    \\
   \includegraphics[width=0.30\linewidth]{diff7.eps}&
  \includegraphics[width=0.30\linewidth]{diff8.eps}
  \\ (c) & (d)
      \\

   \includegraphics[width=0.30\linewidth]{diff9.eps}&
  \includegraphics[width=0.30\linewidth]{diff10.eps}
  \\ (e) & (f)
      \\
   \includegraphics[width=0.30\linewidth]{diff11.eps}&
  \includegraphics[width=0.30\linewidth]{diff12.eps}
  \\ (g) & (h)
\end{tabular}
\end{center}


\caption{Example of depth diffusion by using equation \ref{eq:depth_by_sup}. (a) Depth image, where FSPs have been initialized with a positive value (marked in gray) and the rest of the image with value zero. (b), (c), (d), (e), (f), and (g) depth images corresponding to an increasing number of iterations. After each iteration, the depth difference between corresponding FSPs and BSPs is forced to be at least equal to the initial depth difference $\Delta$, by adding $\Delta$ to FSPs when the difference  between corresponding FSPs and BSPs is less than $\Delta$.(h) Final depth image. }

\label{fig:example_diffusion_with_sup}
\end{figure}






With this strategy, the occlusion boundaries and the relative distances from the viewpoint of depicted objects are simultaneously recovered from local depth information without the need of any explicit segmentation and possible conflicting depth relationships are automatically solved by the diffusion process itself.

Once monocular depth cues are detected, each region involved in a depth relationship is marked by one or few points, called \emph{source points} (see Figure \ref{fig:example_diffusion} (a)). Source points marking the regions closer to the viewpoint are called Foreground Source Points (FSPs), whereas source points marking the regions  more distant to the viewpoint are called Background Source Points (BSPs).
The depth image $z$ is initialized by assigning value $1$ to FSPs and value $0$ to BSPs. The rest of the image is initialized with value zero (see Figure \ref{fig:example_diffusion} (b)).
The diffusion process is applied to the depth image $z$ by using the gradient of the original image  $u$ rather than that of the depth image to define the neighborhood.
By diffusing the depth information using the gray level (color) gradient of the original image instead of that of the depth map itself, the edge information is incorporated into the depth map, assuring that discontinuities in depth are consistent with gray level (color) discontinuities.


The  depth diffusion filter (DDF) proposed in \cite{{dimiccoli:08}}  is
\begin{equation}
%DDF_{h,\rho} z(x) = \frac{1}{C(x)}\sup_{y \in B_\rho(x)} z(y) e^\frac{-|u(x)-u(y)|^2}{h^2},
DDF_{h,\rho} z(x) = \frac{1}{C(x)}\int_{B_\rho(x)}z(y)e^\frac{-|u(x)-u(y)|^2}{h^2}  dy,
\label{eq:depth_by_average}
\end{equation}
where $B_{\rho}(x)$ is a ball of center $x$ and radius $\rho$, $h$ is the filtering parameter which controls the decay of the exponential function, and
\begin{equation}
C(x)=\int_{B_\rho(x)}e^\frac{-|u(x)-u(y)|^2}{h^2}  dy
\end{equation} is the normalization factor.
%$C(x) = \sup_{y \in B_\rho(x)} e^\frac{-|u(x)-u(y)|^2}{h^2}$.
We set $\rho=1$ and $h=10$.



Equation (\ref{eq:depth_by_average}) is applied iteratively until the stability is attained. After each iteration, the values of FSPs and BSPs are modified in order to hold at least the initial depth gradient. In practice, when the difference between the values of a FSP and the corresponding BSP become less than $1$, the value $1$ is added to the value of the FSP. This constraint corresponds to Neumann internal boundary conditions which are understood as a pre-specified jump on the $c\frac{Dz}{Dn}$ as the boundary is crossed, where $c$ is a positive constant and $n$ is the normal to the boundary. This allows one to handle simple sorting when objects are located in multiple layers.
Figure \ref{fig:example_diffusion} is an example of the diffusion through the DDF. As can be sensed looking at Figures \ref{fig:example_diffusion} (c) and (d), by taking an average of the depth values in a neighborhood, a very large number of iterations would be needed to attain the stability.
To make the diffusion process faster, the following equation is used to initialize the process
\begin{equation}
%DDF_{h,\rho} z(x) = \frac{1}{C(x)}\sup_{y \in B_\rho(x)} z(y) e^\frac{-|u(x)-u(y)|^2}{h^2},
DDF_{h,\rho} z(x) = \sup_{B_\rho(x)}z(y)e^\frac{-|u(x)-u(y)|^2}{h^2}  dy,
\label{eq:depth_by_sup}
\end{equation}
while equation (\ref{eq:depth_by_average}) is used only in the last iterations (see Figure \ref{fig:example_diffusion_with_sup}).



In the case of occlusion  there is also a depth order between the two regions separated respectively by the stem of the T. However,  occlusion does not carry any information about the partial order between the objects respectively in partial occlusion and the background. This depth order can be inferred by other cues, such as convexity. When
information about this partial order is present, the depth gradient between one of the BSPs and the
FSPs increases. This is the reason for which force source points are forced to hold "at least" the initial
depth gradient.



\begin{figure}
\begin{center}
\begin{tabular}{cccrrrrrrr}
%\includegraphics[width=0.2\linewidth]{convexity_ori.eps}&
\includegraphics[width=0.25\linewidth]{ori_natural.eps}&
\includegraphics[width=0.25\linewidth]{cues_natural.eps}&
\includegraphics[width=0.25\linewidth]{res_natural.eps}




\\

%\includegraphics[width=0.2\linewidth]{junctions_convexity.eps}&

\includegraphics[width=0.25\linewidth]{ori.eps}&
\includegraphics[width=0.25\linewidth]{local_cues.eps}&
\includegraphics[width=0.25\linewidth]{res.eps}



\\
%\includegraphics[width=0.2\linewidth]{conv_thr_4.eps}\\



\includegraphics[width=0.3\linewidth]{horse.eps}&
\includegraphics[width=0.31\linewidth]{horse_vectors_conv.eps}&
\includegraphics[width=0.31\linewidth]{horse_depth.eps}


\\
\includegraphics[width=0.25\linewidth]{self_occlusion_ori.eps}&
\includegraphics[width=0.25\linewidth]{mod_junctions_self_occlusion.eps}&
\includegraphics[width=0.25\linewidth]{res_self.eps}
\\
($a$) & ($b$) &($c$)

\end{tabular}
\end{center}
  \caption{($a$) Original image. ($b$) Local depth cues are represented through vectors that point to the region closer to the viewpoint.  ($c$) Depth image. }
\label{fig:experiments3}
\end{figure}

% \includegraphics[width=0.2\linewidth]{failure_ori.eps}&
% \includegraphics[width=0.2\linewidth]{failure_lsd.eps}&
% \includegraphics[width=0.2\linewidth]{vectors_failure.eps}&
% \includegraphics[width=0.2\linewidth]{failure_res.eps}\\

Experimental results on real images (see Fig. \ref{fig:experiments3}) have proved that this simple formulation turns out to be very effective for the integration of several monocular depth cue such as occlusion, transparency, convexity, visual completion (both a-modal and modal) and self-occlusion. In particular, contradictory information given by conflicting depth cues is dealt with correctly by the bilateral diffusion mechanism, which permits two regions to invert harmoniously their depths, in full agreement with the phenomenology.

For each  experiment  on Fig. \ref{fig:experiments3}  are shown  four images: the original image (Fig. \ref{fig:experiments3}($a$)); the image where the initial depth gradient at depth cue points is represented through vectors pointing to the region closer to the viewpoint (red vectors arise from T-junctions, green vectors arise from local convexity and each of them represents the point having the biggest curvature value of the connected components obtained by thresholding the curvature) (see Fig.\ref{fig:experiments3} ($b$)); the depth image obtained by performing the bilateral diffusion method (see \ref{fig:experiments3} ($c$)). The depth map is rendered through gray level values (high values indicate regions that are close to the camera).
On the first and the second rows of Figure \ref{fig:experiments3} there are examples of indoor scenes, for which a proper solution is obtained. Instead,  on the third row there is an example of outdoor scene involving a conflict. The T-junction detected on the back of the horse is due to a reflectance discontinuity and its local depth interpretation is incorrect. However, on the depth map, the shape of the horse appear clearly on the foreground since the diffusion process allowed to  overcome the local inconsistency.
In the last row there is an example involving self-occlusion: occluding
contours have different depth relationships at different points along its continuum.
However, the bilateral diffusion method performs well also in this ambiguous situation.


%
%\subsection{Generalized bilateral filter for segmenting an image and diffusing the depth information.}
%\label{subsec:segmentation_with_depth}
%This section will describe another use of bilateral filters to actually provide an accurate image segmentation by diffusing according to local pixel statistics \cite{dimiccoli:09a}. Interestingly, the bilateral strategy serves twice,  first for
%image clue classification and by the depth diffusion strategy described in the preceding section.  This reference develops an automatic algorithm that uses only the cue of occlusion to infer a global, consistent depth ordering.
%The   strategy consists in first detecting T-junctions, then in segmenting the image preserving the T-junctions previously detected, and finally in depth ordering the regions of the final partition by relying on the depth information provided by T-junctions. The global depth ordering is achieved through a graph formalization, which allows to easily detect and solve possible conflicting interpretations.
%
%The region-merging process in \cite{dimiccoli:09a} is based on the work in \cite{calderero:10}, which has recently proposed a new region model based on color histogram and a new family of statistical similarity measures between the region models based on information theory. The region modeling based on empirical distribution has demonstrated a noticeable improvement with respect to first order statistical models where mean or median color values are used as region model since they do not assume that regions are homogeneous in color nor texture. However, the merging process starts by considering that each pixel is a single region, which is modeled still deterministically by its color value and therefore the effect of the statistical modeling become really important only in the late stages of the merging process.
%In the context of T-junction detection, since the region to be segmented is a small neighborhood of a  given candidate point, the importance of a good modeling in early stages of the merging process becomes crucial.
%
%In \cite{dimiccoli:09b}, the same authors have proposed  to solve this problem by modeling each pixel statistically by a probability distribution, instead of deterministically by its color value.  The probability distribution of a given pixel is obtained by exploiting self-similar structures, which can be
%detected by patch comparisons. By \emph{self-similarity}, we mean that every small patch in a natural image has many similar patches in the same image. This assumption is intuitively true for natural images since most objects in the real world have a self-similar or periodic structure: different parts of the same object show the same statistical properties at many different locations, as for instance in correspondence of objects contours.
%The fact that natural images have such a self-similarity property is a kind of stationarity assumption, actually more general and more accurate than any existing image statistics since it does not rely on an underlying model but directly on the data itself.
%
%Under the stationarity assumption, the probability distribution of a single pixel can be computed as follows.
%Let $U$ be an image, $x$ a pixel of the image domain and let $\mathcal{N}(x)$ be a square image patch centered at $x$ which does not include $x$.
%Let us assume that the probability distribution of $U(x)$  depends only on the values of the pixels in $\mathcal{N}(x)$ and it is  independent of the rest of the image (Markovian model). Then, the probability distribution of $U(x)$ given the pixel values of its neighborhood $\mathcal{N}(x)$, can be estimated by computing the set
%
%
%$\Gamma(x)= \lbrace y: \frac{d(\mathcal{N}(x),\mathcal{N}(y))}{d(\mathcal{N}(x),\mathcal{N}_{best})\rbrace}<(1 + \epsilon) \rbrace$,
%where $d(\mathcal{N}(x),\mathcal{N}(y))$ is a distance between a patch centered at $x$ and a patch centered at another pixel location $y$ of the image domain, $\mathcal{N}_{best}$ is the patch that gives the best patch match and $\epsilon$ is a small constant. The histogram of all pixel values in $\Gamma(x)$ gives an estimation of the probability distribution of the value of $x$ given the values of its neighborhood $\mathcal{N}(x)$ \cite{Efros99}. More precisely, assuming that the pixel values range from $1$ to $L$, the histogram is constructed by adding one to the value $U(y)$ for each $y \in \Gamma(x)$ and then normalizing so that the histogram integral is equal to one.
%However, setting a hard threshold ($1 + \epsilon$) to defines the set  $\Gamma(x)$ leads to be able to estimate the probability distribution only of pixels for which similar patches can be found.
%Indeed, in the absence of similar patches, the hard threshold strategy would leave the set $\Gamma(x)$  empty.
%To overcome such a problem, Buades et al. \cite{buades:05} proposed to use an exponential function, which allows a more continuous  distribution.
%More precisely, the probability distribution of a pixel $x$ conditioned to its neighborhood $\mathcal{N}(x)$, can be computed by computing for each pixel $y$ the quantity:
%
%\begin{equation}
%w(x,y)= \frac{1}{Z(x)} e^\frac{-d(\mathcal{N}(x),\mathcal{N}(y))}{h},
%\end{equation}
%where  $Z(x)$ is the normalizing factor:
%\begin{equation}
%Z(x) = \sum_{y  \in U} e^\frac{-d(\mathcal{N}(x),\mathcal{N}(y))}{h},
%\end{equation}
%and
%\begin{equation}
%d(\mathcal{N}(x),\mathcal{N}(y)) = \sum_{z  \in \lbrace -1,1,... \rbrace} \frac{(U(x - z) - U(y - z))^2}{K(z)},
%\end{equation}
%is the similarity between pixel values of a patch centered at $x$ and a patch centered at $y$.
%The variable $z$ indicates the displacement on the patch with respect to $x$, and $K$ is a Gaussian-like function decaying from the center of the patch to its boundary. More precisely,  $K(z)=(2*\lVert z \rVert +1)^2$ acts as a weight function of the euclidean distance between two patches. The goal of the function $K$ is to give more importance on the patch to pixels closer to the reference pixel. Indeed, since we would like to compare local structure as accurately as possible, the error for nearby pixels in the patch should be greater than that of distant pixels.
%The parameter $h$ controls the decay of the exponential function, and therefore of the function $w$. Due to the fast decay of the exponential term, large euclidean distances lead to nearly zero weights acting as an automatic threshold.
%To reduce the computation cost, the search for similar patches is restricted to a search window of size $S\times S$.
%The histogram corresponding to the probability distribution of $x$ is obtained by adding, for each pixel $y$ the value of the function $w(x,y)$ to $U(y)$.
%
%Summarizing, the function $w(x,y)$ associated to the pixel $x$, depends on the similarity $d(\mathcal{N}(x),\mathcal{N}(y))$ between a patch centered at $x$ and a patch centered at $y \in S \times S$ and satisfies the conditions $0< w(x,y) <1$ and $\sum_{y \in \mathcal x} w(x,y) =1$.
%As the size of the image grows, for every pixel $x$, $w(x,y)$ converges to the conditional expectation of $x$ given its neighborhood $\mathcal{N}(x)$. Indeed, as demonstrated in \cite{buades:05},  $w(x,y)$ can be seen as an instance for the exponential operator of the Naradaya-Watson estimator  \cite{watson:64,nadaraya:64}, which estimates the conditional expectation of a random variable.
%
%Once pixels, which are considered as initial regions have been modeled by their pdf, the region merging algorithm starts to iteratively merge pairs of neighboring regions in $W-\Omega$ following the statistical \emph{Kullback-Leibler merging criterion} (KL) \cite{Calderero08}, until three regions are obtained. The KL merging criterion for color images is defined as follows.
%If $P_a$ and $P_{a \cup b}$ are the color histograms in the color space $(YUV)$, then $ D_{KL}(P_a \parallel P_{a \cup b})$ is given by
%
% $ D_{KL}(P_a \parallel P_{a \cup b}) =   \alpha \cdot D_{KL}(P_{Ya} \parallel P_{Ya \cup Yb}) + (1- \alpha) \cdot (D_{KL}(P_{Ua} \parallel P_{Ua \cup Ub}) + D_{KL}(P_{Va} \parallel P_{Va \cup Vb})). $
%We set $\alpha = \frac{1}{2}$.
%
%
%
%In Fig.\ref{fig:statistic_versus_deterministic_local},  we illustrate the advantage of the above-described pixel modeling in the context of the branch extraction.
%In Fig.\ref{fig:statistic_versus_deterministic_local} (a) are shown the neighborhoods of candidate
%points those branches have to be extracted. In Fig.\ref{fig:statistic_versus_deterministic_local} (b), \ref{fig:statistic_versus_deterministic_local} (c), and \ref{fig:statistic_versus_deterministic_local} (d) are shown the results of the branch extraction when using respectively a region model of order zero, the region model proposed in \cite{Calderero08}, and the statistical
%region modeling proposed in \cite{dimiccoli:09b}. As  can be observed, the proposed pixel modeling yields a significant improvement in terms of  contour definition with respect a deterministic pixel modeling.
%
%\begin{figure}[t]
%\begin{center}
%\begin{tabular}{cccc}
%      \includegraphics[width=21mm]{ori1.eps}&
%      \includegraphics[width=20mm]{model0_1.eps}&
%      \includegraphics[width=25mm]{local_segmentation1_labels_complete.eps}&
%       \hspace{-5mm}\includegraphics[width=21mm]{model1_1.eps}
%\\[-10mm]
%      \includegraphics[width=21mm]{ori2.eps}&
%      \includegraphics[width=20mm]{model0_2.eps}&
%      \includegraphics[width=25mm]{local_segmentation2_labels_complete.eps}&
%      \hspace{-5mm}\includegraphics[width=21mm]{model1_2.eps}
%\\[-10mm]
%     \includegraphics[width=21mm]{ori3.eps}&
%      \includegraphics[width=20mm]{model0_3.eps}&
%      \includegraphics[width=25mm]{local_segmentation3_complete.eps}&
%      \hspace{-5mm}\includegraphics[width=21mm]{model1_3.eps}
%     \\ (a) & (b) & (c) & (d)
%\end{tabular}
%\end{center}
%\vspace{-5mm}
%\caption{\small{ Example of comparison: (a) Neighborhood of a candidate point to be segmented. (b) Local segmentation by using a region model of order zero. (c) Local segmentation by using the region model proposed in \cite{Calderero08}. (d) Local segmentation by using the statistical region modeling proposed in \cite{dimiccoli:09b}.} }
%\label{fig:statistic_versus_deterministic_local}
%\end{figure}
%
%In Fig.\ref{fig:statistic_versus_deterministic_global} we show the results of the segmentation obtained by modeling each pixel deterministically (see Fig.\ref{fig:statistic_versus_deterministic_global} (b)) and by modeling each pixel statistically (see Fig.\ref{fig:statistic_versus_deterministic_global} (c)). As  can be observed, the pixel modeling described above gives much accurate results.
%
%
%
%\begin{figure}
%\begin{center}
%\begin{tabular}{cccccrrrr}
%      \includegraphics[width=30mm]{tigre.eps}&
%  \includegraphics[width=30mm]{tigre_vectors.eps}&
%      \includegraphics[width=30mm]{bhat_noarea_tigre_without_nlmeans.eps}&
%         \includegraphics[width=30mm]{merseq_yuv_iid_kl_noarea_nbins_50_alpha_15_tigre_h_50nwin1_nreg5.eps}&
% \includegraphics[width=30mm]{tigre_depth2.eps}
%\\
%      \includegraphics[width=30mm]{paysaje.eps}&
%  \includegraphics[width=30mm]{paysaje_junction_vectors.eps}&
%      \includegraphics[width=30mm]{bhat_noarea_paysaje_without_nlmeans.eps}&
%      \includegraphics[width=30mm]{merseq_yuv_iid_bhat_noarea_nbins_50_alpha_15_paysaje_h_50nwin1_nreg7.eps}&
%\includegraphics[width=30mm]{paysaje_depth.eps}
%\\
%      \includegraphics[width=30mm]{horse.eps}&
% \includegraphics[width=30mm]{horse_junction_vectors.eps}&
%      \includegraphics[width=30mm]{bhat_noarea_horse_without_nlmeans.eps}&
%      \includegraphics[width=30mm]{labeled_horse.eps}&
%\includegraphics[width=30mm]{horse_depth_region_merging.eps}
%     \\ (a) & (b) & (c)  & (d)
%\end{tabular}
%\end{center}
%\vspace{-5mm}
%\caption{\small{ Example of comparison: (a) Original image to be segmented. (b) T-junction detection. (c) Segmentation which preserves TJs by using the region model proposed in \cite{Calderero08}. (d) Segmentation which preserves TJs by using the statistical region modeling proposed in \cite{dimiccoli:09b}. (e) Depth map.} }
%\label{fig:statistic_versus_deterministic_global}
%\end{figure}


\paragraph{Acknowledgements}
 The David raw point set is courtesy of the the \emph{Digital Michelangelo Project, Stanford University}. Fragment '31u' of the \emph{Stanford Fragment Urbis Romae, Stanford University}, the Screw Nut point set is provided by the AIM@SHAPE repository and is courtesy of Laurent Saboret, INRIA. Research  partially financed by Institut Farman, ENS Cachan,  the Centre National d'Etudes Spatiales
(MISS Project),  the European Research Council (advanced grant Twelve Labours) and  the
Office of Naval research (grant N00014-97-1-0839).

\bibliography{bib_review,biblio} 
 %%\bibliography{biblio}
\end{document}
