\documentclass[a4paper,10pt]{article}

\title{Fattening free block matching}
\date{}

\author{G. Blanchet \thanks{CNES, 18 avenue Edouard Belin, Toulouse 31401, France}
\and A. Buades \thanks{Dpt Matematiques Informatica, Universitat Illes Balears, Ctra Valldemossa km 7.5, Palma de Mallorca, Spain} \thanks{MAP5, CNRS - Universit\'{e} Paris Descartes, 45 rue des Saints P\`{e}res, 75270 Paris Cedex 06, France}
 \and B. Coll$^\dag$
 \and J.M Morel \thanks{CMLA, ENS Cachan, 61 av. Pr\'{e}sident Wilson, Cachan 94235, France}
 \and B. Roug\'{e} \thanks{CESBIO, 18 Av. Edouard Belin,  Toulouse 31400, France}
}

\usepackage{graphicx, epsfig}
\usepackage{amsmath}
\usepackage{amssymb,gensymb}
\usepackage[latin1]{inputenc}
\usepackage{color,epic,eepic}
\usepackage{subfigure}
\usepackage{caption}
\usepackage{color}
\usepackage{minitoc}
\usepackage{nicefrac}

\newtheorem{definition}{Definition}
\newtheorem{theorem}{Theorem}
\newtheorem{lemma}{Lemma}
\newtheorem{proposition}{Proposition}
\newtheorem{corollary}{Corollary}
\newenvironment{proof}[1][Proof]{\begin{trivlist}
\item[\hskip \labelsep {\bfseries #1}]}{\end{trivlist}}
\newenvironment{remark}[1][Remark]{\begin{trivlist}
\item[\hskip \labelsep {\bfseries #1}]}{\end{trivlist}}

\def\bc{{\mathbf c}}
\def\bB{{\mathbf B}}
\def\bq{{\mathbf q}}
\def\br{{\mathbf r}}
\def\bx{{\bf x}}
\def\by{{\bf y}}
\def\bn{{\bf n}}
\def\bm{{\bf m}}
\def\bk{{\bf k}}
\def\bs{{\bf s}}
\def\bt{{\bf t}}
\def\a{{a}}
\newcommand{\Var}{\operatorname{Var}}
\newcommand{\Cov}{\operatorname{Cov}}

\newcommand{\sinc}{{\ensuremath{\operatorname{sinc}}}}
\newcommand{\half}{{\ensuremath{\nicefrac{s}{2}}}}
\newcommand{\R}{\mathbb{R}} %real numbers
\newcommand{\Z}{\mathbb{Z}} %integer numbers
\newcommand{\C}{\mathbb{C}} %complex numbers
\newcommand{\E}{\mathbb{E}} %expectation
\newcommand{\w}{{\varphi_{\bq}}} % correlation window around x_0
\newcommand{\scalar}[2]{{\left\langle\,#1,\,#2\right\rangle}} % scalar product
\newcommand{\scalarX}[3]{{\scalar{#1}{#2}_{#3}}} % weighted scalar product
\newcommand{\scalarphi}[2]{\scalarX{#1}{#2}{\w}} % \varphi_x0 weighted scalar product
\newcommand{\normX}[2]{\| #1 \|_{#2}}
\newcommand{\normphi}[1]{\| #1 \|_\w}
\newcommand{\Proba}{\mathbb{P}}

\newcommand{\comment}[1]{{\textbf\textit #1}}

\newcommand{\utilde}{\tilde{u}}




\begin{document}

\maketitle

\begin{abstract}  Block matching along epipolar lines  is the core of most stereovision algorithms in geographic information systems. The usual distances between blocks are the sum of squared distances in the block (SSD) or the correlation. Minimizing these distances causes the {\it fattening} effect, by which the center of the block inherits the disparity of the more contrasted pixels in the block. This fattening error  occurs everywhere in the image, and not just on strong depth discontinuities. The fattening effect at strong depth edges is is a particular case of fattening, called  {\it foreground fattening} effect. A theorem proved in the present paper shows that  a simple and universal adaptive weighting of the SSD  resolves the fattening problem at all smooth disparity points\footnote{A Spanish patent has been applied for by Universitat de Illes Balears  \cite{PatentBCMR09}}.  The optimal SSD weights are nothing but the inverses of the squares of the image gradients in the epipolar direction. With these adaptive weights, it is shown that the optimal  disparity function  is the result of the convolution of the real disparity with a prefixed kernel. Experiments on simulated and real pairs prove that the method does what the theorem predicts, eliminating surface bumps caused by fattening. However, the method does not resolve the foreground fattening.
\end{abstract}
\graphicspath{{./figures/}}


\section{Introduction}

Stereovision consists in finding the depth of a scene from several views of it. This is one of the central problems in computer vision, and it has been an active object of research in the last forty years. Stereovision is based on the fact that differences of depth in a 3D scene create geometrical disparities  between  views of the same scene taken from different points of view.


Given two stereo rectified images $u$ and $v$, the question reduces to finding a disparity function  $\epsilon$ such that
$u(x) = v(x+ \epsilon(x))$.
 Like in motion estimation, the above equation presents the aperture problem, namely the ambiguity of the solution, even when some regularity is demanded for the  disparity.  For this reason, many stereovision algorithms do not look for a function $\epsilon$ matching the grey level intensity of each pixel. They prefer to  compare  the grey levels of an entire block around each pixel. The simplest resulting algorithm is known as block matching by SSD (sum of squared distances).

The most important drawback of SSD is the well known ``fattening effect''. According to Kanade et Okutomi  \cite{Kanade94},

\begin{quotation} A central problem in stereo matching by computing correlation or sum of squared differences (SSD) lies in selecting an appropriate window size. The window size must be large enough to include enough intensity variation for reliable matching, but small enough to avoid the effects of projective distortion. If the window is too small and does not cover enough intensity variation, it gives a poor disparity estimate, because the signal (intensity variation) to noise ratio is low. If, on the other hand, the window is too large and covers a region in which the  depth of scene points (i.e. disparity) varies, then the position of maximum correlation or minimum SSD may not represent correct matching due to different projective distortion in the left and right images.  The fattening effect occurs when the selected window contains pixels at different depth. In that case we cannot find exactly the same window and the obtained disparity depends on the different disparities of the window and not only the central pixel itself.\end{quotation}



% In satellite imaging this phenomenon is mostly noticeable as a dilation of buildings by a width of half the correlation window.
% However, as is observed in Fig. \ref{fig:house_adhesion} the fattening effect has also an internal effect on buildings.


%\begin{figure}[h]
%\centering
%\epsfig{file=scheme_adhesion.eps, height = 2.5 cm}
%\epsfig{file=scheme_house.eps, height = 2.5 cm}
%\caption{Continuous line: cut of the depth map of a house. Dotted line: the estimated depth affected with fattening when the border of the house coincides with a gray level discontinuity. Fattening provokes a dilation of the building, but also a flattening of the roof near the border.}
%\label{fig:house_adhesion}
%\end{figure}


The usual way to cope with the fattening effect is to use adaptive windows that avoid image discontinuities as was first proposed by Kanade et al \cite{Kanade94}. Similar works pre-computing edge points and recursively growing a comparison window avoiding them were proposed by Lotti et al. \cite{Lotti94} and  recently by Wang et al. \cite{Wang06}.  Patricio et al. \cite{patricio2004similarity} and Yoon et al. \cite{Yoon06} select an adaptive window containing only pixels with a  grey level similar to the reference one, like in neighborhood and bilateral filters  \cite{tomasi1998bilateral, yaroslavsky2003fundamentals}.

Other approaches do not try to avoid the discontinuities of the image. They select an adaptive window with a minimum distance criterion.  The subjacent idea is that windows which do not contain discontinuities will be matched with a small window distance.
Fusiello et al. \cite{fusiello2000symmetric} choose among all the windows containing the reference pixel the one which has a minimal distance with its corresponding  one in the second image. Veksler \cite{Veksler03b} applied the same strategy but used in addition square windows of different sizes. A more elaborated version by  Hirschmuller et al \cite{Hirschmuller02} adapts the shape of the window by dividing the comparison window into small sub-windows and taking those which attain the minimum distance. The Delon et al. \cite{Delon07} paper proposes a different strategy, the barycentric correction attributing the disparity of a window to the window barycenter pondered by the image gradients.

Point feature matching methods overcome the fattening problem at the cost of a drastic reduction of the match density. Matched features can also be curvilinear, which also circumvents
the fattening problem to some extent. For instance, Schmid \cite{Schmid00} describes a set of algorithms for automatically matching individual line segments and curves. Robert \cite{Robert91} presents an edge-based stereovision algorithm, where the primitives to be matched are cubic B-splines approximations of the 2-D edges. Mus\'{e} et al. \cite{Muse06} and Cao et al. \cite{Cao07} discuss how to
automatically match pieces of level lines and extract coherent groups of such matches. The Matas et al. \cite{matas04} MSER method solves the problem by matching stable and homogeneous image regions, but their match set is  again sparse. Even if features may seem more local, they depend anyway on a broad neighborhood.  The same remark applies to the SIFT method (Lowe \cite{Lowe04}) and their affine invariant extensions  \cite{morel2009asift}.  Even if the fine scale Laplacian extrema used (e.g.) in the SIFT method are very local, their descriptor around involves anyway a $8\times 8$ window (see \cite{cao2008theory, kimmelmser10, mikolajczyk2005comparison} for comparison on MSER and SIFT). Thus the fattening problem can occur anyway with these methods.


The fattening effect is not the sole obstacle to a correct disparity computation.  Occlusions and moving objects make it a very difficult and sometimes ill-posed problem. Taking simultaneous snapshots with a low baseline avoids partially these drawbacks. However, when using a low baseline  a larger precision in the disparity computation is needed to get the same depth precision. The use of a low B/H (where B is the baseline and H is the altitude) was proposed in satellite imaging by Delon and Rouge \cite{Delon07}.


\section{Mathematical analysis of SSD}

Let us denote by $\bx=(x,y)$  an image point in the continuous image domain, and by $u_1(\bx)=u_1(x,y)$ and $u_2(\bx)$ the images of an ortho-rectified stereo pair. Assume that the epipolar direction is the $x$ axis. The underlying depth map can be deduced from the disparity function $\varepsilon(\bx)$ giving the shift of an observed physical point $\bx$ from the left image $u_1$ in the right image $ u_2$. The physical disparity $\varepsilon(\bx)$ is not well-sampled. Therefore, it cannot be recovered at all points, but only essentially at points $\bx$ around which the depth map is continuous. Following the formulation by Delon and Rouge \cite{Delon07} and Sabater \cite{ThesisNeus}, around such points, the deformation model from an image to the other is
\begin{align}\label{disparityequation}
u_1(\bx) &= u(x+\varepsilon(\bx),y) +n_1(\bx)\; \nonumber \\
u_2(\bx) &= u(\bx)+n_2(\bx),
\end{align}
where $u$ the true scene image and $n_1(\bx)$ and $n_2(\bx)$  independent Gaussian white noises with standard deviation $\sigma$. (The captor noises are independent because the snapshots are different.)
Block matching amounts to finding the disparity at $\bx_0$ minimizing
\begin{equation}\label{minimizationproblem}
e_{\bx_0}(\mu)=\int_{[0,N]^2}\varphi(\bx-\bx_0)\big(u_1(\bx)-u_2(\bx+(\mu,0))\big)^2d\bx.\end{equation}
%
where $\varphi(\bx-\bx_0)$ is a soft window function centered at $\bx_0$. For a sake of compactness in notation, $\varphi_{\bx_0}(\bx)$ stands for $\varphi(\bx-\bx_0)$,  $\int_{\varphi_{\bx_0}} u(\bx)d\bx$ will be an abbreviation for $\int\varphi(\bx-\bx_0)u(\bx)d\bx$; we will write $u(\bx+\mu)$ for $u(\bx+(\mu,0))$ and $\varepsilon$ for $\varepsilon(\bx)$. The  minimization problem  \eqref{minimizationproblem} rewrites
$$\min_\mu\int_{\varphi_{\bx_0}}\big(u(\bx+\varepsilon(\bx) ) +n_1(\bx)-u(\bx+\mu)-n_2(\bx+\mu)\big)^2 d\bx.$$
Differentiating this energy with respect to $\mu$ implies that any local minimum $\mu=\mu(\bx_0)$ satisfies

\begin{equation}\label{derivada}
\int_{\varphi_{\bx_0}}\Big(u(\bx+\varepsilon(\bx))+n_1(\bx)-u(\bx+\mu)-n_2(\bx+\mu)\Big) \times
\Big(u_x(\bx+\mu)+(n_2)_x(\bx+\mu)\Big) \,d\bx  = 0.
\end{equation}
One has by Taylor-Lagrange formula
$u_x(\bx+\mu)=(u_x(\bx+\varepsilon))+O_1(\mu-\varepsilon)$, with \begin{equation}\label{TaylorLagrange}O_1(\mu-\varepsilon)\leq |\mu-\varepsilon|\max|u_{xx}(\bx+\varepsilon)|\end{equation}
 and
$u(\bx+\varepsilon(\bx))-u(\bx+\mu)=u_x(\bx+\varepsilon) (\varepsilon-\mu)+O_2((\varepsilon-\mu)^2),$ where $$|O_2((\varepsilon-\mu)^2)|\leq\frac 12 \max|(u_{xx}(\bx+\varepsilon))|(\varepsilon-\mu)^2\,.$$
Thus equation (\ref{derivada}) yields

\begin{align}
\int_{\varphi_{\bx_0}}\Big(u_x(\bx+\varepsilon)(\varepsilon-\mu)+O_2((\varepsilon-\mu)^2)+n_1(\bx)-n_2(\bx+\mu)\Big) \times &  \nonumber \\
\quad\quad\quad\quad\quad \Big(u_x(\bx+\varepsilon)+O_1(\mu-\varepsilon)+
 (n_2)_x(\bx+\mu)\Big) \, d\bx & = 0.
\end{align}
and therefore

\begin{equation}\label{formuledonnantm}\mu\int_{\varphi_{\bx_0}} (u_x(\bx+\varepsilon))^2 d\bx =\int_{\varphi_{\bx_0}} (u_x(\bx+\varepsilon))^2\,\varepsilon(\bx) \,d\bx +\tilde{\mathcal A} + \tilde{\mathcal B}
+\mathcal O_1 +\mathcal O_2,\end{equation}
where

\begin{eqnarray}
\tilde{\mathcal A} & = & \int_{\varphi_{\bx_0}} u_x(\bx+\varepsilon)  \big(n_1(\bx)-n_2(\bx+\mu)\big)d\bx ;\\
\tilde{\mathcal B} & = & \int_{\varphi_{\bx_0}} \big(n_1(\bx)-n_2(\bx+\mu)\big)(n_2)_x(\bx+\mu) d\bx ; \\
\mathcal O_1       & = & \int_{\varphi_{\bx_0}} u_x(\bx+\varepsilon)(\varepsilon-\mu)(n_2)_x(\bx+\mu)d\bx  \nonumber\\
                   &   & + \int_{\varphi_{\bx_0}} O_1(\mu-\varepsilon)\big(n_1(\bx)-n_2(\bx+\mu)\big) d\bx ;\\	
\mathcal O_2       & = & \int_{\varphi_{\bx_0}} O_2(\varepsilon-\mu)^{2}(u_x(\bx+\varepsilon)) d\bx  \nonumber\\
		   &   & + \int_{\varphi_{\bx_0}} O_2(\varepsilon-\mu)^{2} [O_1(\mu-\varepsilon) + (n_2)_x(\bx+\mu)] d\bx \nonumber\\
                   &   & + \int_{\varphi_{\bx_0}} O_1(\mu-\varepsilon)(u_x(\bx+\varepsilon))(\varepsilon-\mu) d\bx.
\end{eqnarray}
Denote by $\overline{\varepsilon}$ the average of $\varepsilon$ on the support of $\varphi(\bx-\bx_0)$, denoted by $B_{\bx_0}.$
By the Taylor-Lagrange theorem we have
$$\tilde{\mathcal A} = \mathcal A +\mathcal O_\mathcal A$$
where
\begin{equation}\label{defdeA}\mathcal A=\int_{\varphi_{\bx_0}} u_x(\bx+\varepsilon) \big(n_1(\bx)-n_2(\bx+\overline{\varepsilon})\big) d\bx\end{equation} and

\begin{equation}\label{defdeOA}\mathcal O_\mathcal A=(\overline{\varepsilon}-\mu)\int_{\varphi_{\bx_0}} (u_x(\bx+\varepsilon))(n_2)_x(\bx+\tilde\varepsilon(\bx))d\bx,\end{equation}
where $\tilde\varepsilon(\bx)$ satisfies $\tilde\varepsilon(\bx)\in[\min(\mu,\overline{\varepsilon}), \max(\mu,\overline{\varepsilon})].$
In the same way,
$$\tilde{\mathcal B}=\int_{\varphi_{\bx_0}}\big(n_1(\bx)-n_2(\bx+\mu)\big)
(n_2)_x(\bx+\mu) d\bx$$   so that
$\tilde{\mathcal B}=\mathcal B +\mathcal O_{\mathcal B},$ where

\begin{equation}\label{defdeB}\mathcal B = \int_{\varphi_{\bx_0}}\big(n_1(\bx)-n_2(\bx+\overline{\varepsilon})\big)
(n_2)_x(\bx+\overline{\varepsilon}) d\bx\end{equation} and

\begin{equation}\label{defdeOB}\mathcal O_\mathcal B=(\mu-\overline{\varepsilon})\int_{\varphi_{\bx_0}}  n_1(\bx)(n_2)_{xx} (\bx+\tilde\varepsilon(\bx)) - (n_2(n_2)_{x})_x(\bx+\tilde\varepsilon(\bx)) d\bx.\end{equation}
The terms $\mathcal A$ and $\mathcal B $    are stochastic and we must estimate their expectation and variance. The terms $\mathcal O_1$, $\mathcal O_2$, $\mathcal O_{\mathcal A}$,
  $\mathcal O_{\mathcal B}$  are higher order terms with respect to $\varepsilon-\mu$ and are   negligible if $\varepsilon-\mu$ is small, and  the noise samples bounded.

\begin{lemma} \label{estimatevariance} Consider the main error terms
 $$\mathcal A=\int_{\varphi_{\bx_0}} u_x(\bx+\varepsilon(\bx))  \big(n_1(\bx)-n_2(\bx+\overline{\varepsilon})\big) d\bx$$ and
 $$\mathcal B = \int_{\varphi_{\bx_0}}\big(n_1(\bx)-n_2(\bx+\overline{\varepsilon})\big)
(n_2)_x(\bx+\overline{\varepsilon}) d\bx$$ as defined above. One has ${\bf E}\mathcal A={\bf E}\mathcal B=0$ and
\begin{align*}
\Var(\mathcal A) &  =2\sigma^2\int\left[\varphi(\bx-\bx_0)u_x(\bx+\varepsilon)\right]^2_Nd\bx\\
                 &\leq 2\sigma^2\int\varphi(\bx-\bx_0)^2 u_x(\bx+\varepsilon)^2; \\
\Var(\mathcal B) & \leq \frac{2\pi^2 \sigma^4}3 \int \varphi (\bx-\bx_0)^2 d\bx  + \sigma^4 \int \varphi_x (\bx-\bx_0)^2 d\bx.
\end{align*}

\end{lemma}

\noindent{\bf Proof:}
Notice that $n_1(\bx)$ and $n_2(\bx+\overline{\varepsilon})$ are independent Gaussian noises with variance $\sigma^2$. Thus their difference is again a Gaussian noise with variance $2\sigma^2$.  It therefore follows  that
$$\Var(\mathcal A)=2\sigma^2\int\left[\varphi(\bx-\bx_0)u_x(\bx+ \varepsilon))\right]^2_Nd\bx\leq 2\sigma^2 \int\varphi(\bx-\bx_0)^2(u_x(\bx+ \varepsilon))^2d\bx.$$
$$\Var(\mathcal  B)\leq 2\left[ \Var(\int_{\varphi_{\bx_0}} n_1(\bx)(n_2)_x(\bx+\overline{\varepsilon})
+\Var(\int_{\varphi_{\bx_0}} n_2(\bx+\overline{\varepsilon})(n_2)_x(\bx+\overline{\varepsilon}))\right]$$
$$\leq 2\left[\sigma^2\times \frac{\pi^2\sigma^2}3\int\varphi^2(\bx-\bx_0)+\frac{\sigma^4}2\int\varphi_x(\bx-\bx_0)^2\right]
$$
$$=\frac{2\pi^2\sigma^4}3\int\varphi(\bx-\bx_0)^2+\sigma^4\int\varphi_x(\bx-\bx_0)^2.$$


\begin{theorem}\label{Maintheoremdisparity} {\bf (Main disparity formula and exact noise error estimate)} Consider an optimal disparity $\mu(\bx_0)$ obtained as any absolute minimizer of $e_{\bx_0}(\mu)$ (defined  by \eqref{minimizationproblem}).
Then
\begin{equation}\label{formuledonnantm1}\mu(\bx_0)=
\frac{\int_{\varphi_{\bx_0}} [u_x(\bx+\varepsilon(\bx))]^2\varepsilon(\bx) d\bx}
{\int_{\varphi_{\bx_0}} [u_x(\bx+\varepsilon(\bx))]^2 d\bx}
 +\mathcal E_{\bx_0} +\mathcal F_{\bx_0} + \mathcal O_{\bx_0}\end{equation}
 where
 $$\mathcal E_{\bx_0}=\frac{\int_{\varphi_{\bx_0}} (u_x(\bx+\varepsilon(\bx)) \big(n_1(\bx)-n_2(\bx+\overline{\varepsilon})\big) d\bx}
 {\int_{\varphi_{\bx_0}} [u_x(\bx+\varepsilon(\bx))]^2 d\bx}$$ is the dominant noise term,

 $$\mathcal F_{\bx_0}=\frac{ \int_{\varphi_{\bx_0}}\big(n_1(\bx)-n_2(\bx+\overline{\varepsilon})\big)
(n_2)_x(\bx+\overline{\varepsilon}) d\bx}
 {\int_{\varphi_{\bx_0}} [u_x(\bx+\varepsilon(\bx))]^2 d\bx}$$
  and
 $\mathcal O_{\bx_0}$ is made of smaller terms. In addition the variances of the main error terms due to noise satisfy
 \begin{equation}\label{varianceofthenoiseterm1}
 \Var (\mathcal E_{\bx_0})=   2\sigma^2  \dfrac{ \int\left[\varphi(\bx-\bx_0)u_x(\bx+\varepsilon)\right]^2_Nd\bx}{\Big( \int \varphi(\bx-\bx_0)u_x(\bx+\varepsilon)^2d\bx\Big)^2};
 \end{equation}

  \begin{equation}\label{varianceofthenoiseterm2}
 \Var (\mathcal F_{\bx_0})\leq    \dfrac{ \frac{2\pi^2}{3} \sigma^4 \int \varphi (\bx-\bx_0)^2 d\bx  + \sigma^4 \int \varphi_x (\bx-\bx_0)^2 d\bx}{\Big( \int \varphi(\bx-\bx_0)u_x(\bx+\varepsilon)^2d\bx\Big)^2}.
 \end{equation}
 Finally, $$\mathcal O_{\bx_0}=\frac{\mathcal O_1+\mathcal O_2+\mathcal O_{\mathcal A}+
 \mathcal O_{\mathcal B}}{\int_{\varphi_{\bx_0}} [u_x(\bx+\varepsilon(\bx))]^2 d\bx},$$
 and $${\bf E}\mathcal O_{\bx_0} = O(\max_{\bx\in B_{x_0}}|\varepsilon(\bx)-\mu|),$$
 $$\Var(\mathcal O_{\bx_0}) = O(\max_{\bx\in B_{x_0}}|\varepsilon(\bx)-\mu|^2).$$

\end{theorem}

\noindent{\bf Proof:} This result is an immediate consequence of \eqref{formuledonnantm} completed with the variance estimates in Lemma \ref{estimatevariance}. The estimates for the higher order terms $\mathcal O$ are a straightforward application of Cauchy-Schwartz inequality.

\begin{remark} Theorem \ref{Maintheoremdisparity} makes sense only when the optimal disparity $\mu(\bx_0)$ is {\rm consistent}, namely satisfies  for $\bx$ in the  support $B_{\bx_0}$ of $\varphi(\bx-\bx_0)$,
\begin{equation}\label{hypothese1}|\varepsilon(\bx)-\mu(\bx_0)| <<1.\end{equation}   Thus, one of the main steps of block matching must be to eliminate inconsistent matches.
\end{remark}
\begin{remark}  In all treated examples, it will be observed
  that  $\Var(\mathcal B)\ll \Var(\mathcal A)$, which by Lemma  \ref{estimatevariance} directly follows from
 \begin{equation}\label{hypothese2}\sigma^2\left[\frac{2\pi^2}3\int \varphi(\bx-\bx_0)^2+\int\varphi_x(\bx-\bx_0)^2\right]\ll 2\int[\varphi(\bx-\bx_0)u_x(\bx+\overline{\varepsilon})]_N^2.\end{equation}

\end{remark}



\section{Mathematical definition of fattening, and its solution}

The previous mathematical formulation  tells us that the obtained minimizer for the SSD problem satisfies
\begin{equation}\label{eq:simple}
\mu(\bx_0)=
\frac{\int_{\varphi_{\bx_0}} [u_x(\bx+\varepsilon(\bx))]^2\varepsilon(\bx) d\bx}
{\int_{\varphi_{\bx_0}} [u_x(\bx+\varepsilon(\bx))]^2 d\bx}
\end{equation}
up to the noise terms. In other terms, the obtained minimizer will be an center of mass of the disparities at each pixel in the correlation window, each being weighted by its squared image gradient.

This explains the fattening effect, which actually occurs at {\it every pixel}: Whenever a pixel or a cluster of pixels have a large gradient with respect to their neighboring ones,  the estimated disparity for these neighboring pixels will be obtained by combining mainly the disparities of these few very contrasted pixels. It can even happen that a single pixel dominates the estimated disparity for all of its neighboring ones. This effect is strong in all textures and also near image edges, where a line of pixels dominates the SSD of all their neighboring ones. This case causes the so called {\it foreground fattening} phenomenon by which buildings looks fatter than they really are. Yet, the fattening effect happens everywhere, because a gradient barycenter is never exactly the center of the correlation window. Even if this is not very noticeable when looking at the disparity image,  this effect becomes conspicuous when looking at the 3D reconstruction of the estimated depth, (Fig. \ref{fig:base3d}).

The above calculations show that there is only one way to avoid the fattening: It is to remove the disparity imbalance in the comparison window.  One can compensate the effect of the squared gradients in the above integral by directly modifying the values of the window function $\varphi$, making it adaptive. By taking $\varphi_{x_0}(\bx) = \frac{\rho_{\bx_0}(\bx)}{u_x(\bx+\varepsilon(\bx))^2}$ in equation \eqref{eq:simple} we obtain
\begin{equation}\label{eq:simple2}
\mu(\bx_0)=
\frac{\int_{\rho_{\bx_0}} \varepsilon(\bx) d\bx}
{\int_{\rho_{\bx_0}} d\bx},
\end{equation}
which is equivalent to $$\mu(\bx_0) = \int  \rho(\bx - \bx_0) \varepsilon(\bx) d\bx,$$ since the function $\rho$ is normalized to have the integral equal to one. In that way the disparity becomes a weighted average of all disparities in the correlation neighborhood, which is no  more weighted by the image gradient. Therefore, the computed disparity is the convolution of the ground truth disparity $\varepsilon$ with a kernel,  which can incidently be fixed at will.
The most natural choice for the window $\rho$ is an isotropic kernel, for example a Gaussian $G_a$.  If we select such a kernel, the computed disparity writes $G_a \ast \varepsilon$, which can  be interpolated and could even be deconvolved to some extent.  The choice of the size of the window depends primarily on the noise variance. If there were no noise at all the window could be a Dirac.  In presence of noise,
 the dominant disparity error term due to the noise given by  Theorem \ref{Maintheoremdisparity}   rewrites
 \begin{equation}
 \Var (\mathcal E_{\bx_0})=   2\sigma^2   \int  \frac{ \rho(\bx-\bx_0)^2 } {u_x(\bx+\varepsilon(x))^2} d\bx.
 \end{equation}
Thus, the size of the window must be large enough to ensure this value to be low enough to compensate for $\sigma^2$. Indeed, the integral of $\rho$ being 1, the broader the support of $\rho$ the smaller the integral will be, because of  the presence of the $\rho^2$ term. This implies that the integral behaves like $1/n$, where $n$ is the number of pixels in the window.  A good point of the above result is that  adaptive window can be  larger without causing a fattening effect.

The discrete implementation of such an algorithm faces the problem of computing the true derivatives $u_x(\bx+\varepsilon(\bx))$ from the two available images $u_1$ and $u_2$.  We can compute the derivative on the first image, obtaining
$$u_1'(\bx) ^2= (u'(\bx+\varepsilon(\bx)) (1 + \varepsilon'(\bx)) + n_1 ' (\bx) )^2 . $$
Since this is a stochastic term, the right choice must be indicated by  its mean
$$E u_1'(\bx) ^2 = u'(\bx+\varepsilon(\bx))^2 (1 + \varepsilon'(\bx))^2 + 2 \sigma^2.$$
This identity shows that, because of the noise term, we will be only able to compute the actual derivatives if and when  $\varepsilon'(\bx)$ is small.   We shall make this assumption, which means that the relief is smooth. In order to avoid too small gradients due mainly to noise, we shall use the following weighting function
$$\varphi_{\bx_0}(\bx) = \frac{\rho_{\bx_0}(\bx)}{\max(u_x(\bx+\varepsilon(\bx))^2, 6 \sigma^2)},$$
where $\sigma$ is the noise standard deviation.


\section{Comparative experiments}
In order to illustrate and compare the performance of the classical SSD strategy and the proposed adaptive algorithm,  several tests were performed on synthetic and real stereo pairs, and the proposed method was compared with the two most classic fattening correction strategies.

The first experiments were simulated  pairs with a smooth disparity function. The disparity $\varepsilon$ in Fig. \ref{fig:base_gt} was applied to the reference texture images $u$  of Fig. \ref{fig:base}. Each image was warped by $\varepsilon$ to obtain the image pair. Gaussian white noise was added to both images of the pair. Texture images were used to make sure that around each pixel there was enough information to permit its correct matching. The first ground truth disparity varies slowly and smoothly while the other two are more oscillatory.


\begin{figure}[t]
\centering
\epsfig{file=object1.eps, width = 1.75 cm}
\epsfig{file=object2.eps, width = 1.75 cm}
\epsfig{file=object3.eps, width = 1.75 cm}
\epsfig{file=object4.eps, width = 1.75 cm}
\epsfig{file=object5.eps, width = 1.75 cm}
\epsfig{file=object6.eps, width = 1.75 cm}
\caption{Reference image warped by a known disparity to obtain an image pair.} \label{fig:base}
\end{figure}

\begin{figure}
\centering
\epsfig{file=sin2.gt.eps, width = 2.5cm}
\epsfig{file=sin1.gt.eps, width = 2.5cm}
\epsfig{file=sin3.gt.eps, width = 2.5cm}
\caption{Ground truth disparities applied to images in Fig. \ref{fig:base}.} \label{fig:base_gt}
\end{figure}

\begin{figure}
\centering
\epsfig{file=sin2.cor.p7.eps, width = 2.5cm}
\epsfig{file=sin2.cor.p11.eps, width = 2.5cm} \hspace{1 cm}
\epsfig{file=sin2.corb.p7.eps, width = 2.5cm}
\epsfig{file=sin2.corb.p11.eps, width = 2.5cm}


\epsfig{file=sin1.cor.p7.eps, width = 2.5cm}
\epsfig{file=sin1.cor.p11.eps, width = 2.5cm} \hspace{1cm}
\epsfig{file=sin1.corb.p7.eps, width = 2.5cm}
\epsfig{file=sin1.corb.p11.eps, width = 2.5cm}

\epsfig{file=sin3.cor.p7.eps, width = 2.5cm}
\epsfig{file=sin3.cor.p11.eps, width = 2.5cm} \hspace{1cm}
\epsfig{file=sin3.corb.p7.eps, width = 2.5cm}
\epsfig{file=sin3.corb.p11.eps, width = 2.5cm}

\caption{Obtained disparities for the first image in Fig. \ref{fig:base} and the three ground truth disparities in Fig. \ref{fig:base_gt}. The left column shows the disparities obtained with a classical SSD algorithm with an isotropic weighting window of size 7$\times$7 and 11$\times$11. In the right column same experiments, but with the proposed algorithm.} \label{fig:base_est}
\end{figure}

 Fig.  \ref{fig:base_est} presents the disparity maps obtained by both strategies for the first image of the data base. In this case, a noise with standard deviation 1 has been added, yielding a signal to noise ratio of about one hundred. The results with  SSD and with the proposed strategy are shown with prolate functions supported by 7$\times$7 and 11$\times$11 pixels.  Observe that the disparity obtained with the proposed strategy is more similar to the ground truth than the classical SSD algorithm. This improvement is conspicuous when the 11$\times$11  prolate is used or when the disparity map is more oscillatory. This experimental fact is in agreement with the mathematical arguments and formulas developed in the previous section.  The obtained disparity for the classical SSD strategy depends on the true disparity on the 7$\times$7 or 11$\times$11 neighborhood and is weighted by the square of the gradient. Thus, with a larger window the probability of having large gradients on the window is increased and the favored disparity by these large gradient points can be more different than the one of the reference pixel.

In Fig. \ref{fig:base3d} are displayed the three-dimensional representations of the central row in Fig. \ref{fig:base_est} with a 7$\times$7 prolate function. One better evaluates with this representation the difference between the classical and the adaptive SSD.  The surface obtained by the adaptive SSD is smooth and very similar to the ground truth. However, the surface by the classical SSD strategy presents many irregularities due to its dependence on the image gradients.

\begin{figure}
\centering
\epsfig{file=sin2.cor.surf.eps, width = 7 cm}
\epsfig{file=sin2.corb.surf.eps, width = 7 cm}
\caption{Three dimensional representation of the estimated disparity from the middle row of Fig. \ref{fig:base_est}. Top: estimated disparity by SSD with a 7$\times$7 correlation window. Bottom: proposed adaptive SSD with the same 7$\times$7 correlation window. The fattening effect creates evident irregularities in the reconstructed surface.} \label{fig:base3d}
\end{figure}


Table \ref{tb:mse} shows the average Euclidean distance between the obtained disparity and the ground truth for the six images in Fig.  \ref{fig:base}. The error values are very similar when the prolate is small or when the disparity varies slowly, while they increase for the classical SSD algorithm when a larger prolate or an oscillating ground truth is applied.
Table \ref{tb:mse2} shows the error committed by comparing the true normals to the surface of the ground truth with the normals to the surfaces of the obtained disparities. Are shown the ratio of points of the surface for which the normal has an error of more than 10 degrees with respect to the original normal. The accuracy gain is quite important by using the adaptive strategy.
Notice that the distance of normals is the right measure to estimate how two renderings of the same object differ visually. Indeed, most 3D visualizations are done by a Lambertian model. The grey  level of the rendered image  is the scalar product of the surface normal with the solar direction. Thus the above error measure is the right one to estimate the visual gain.

It is observed in Table 2 that with a small correlation window the use of the adaptive strategy is more sensitive to noise. This is not easily explained by comparing the precision terms in Theorem 1, but it can be explained by simple probabilistic arguments. When computing the weighted Euclidean distance of two noisy patches, the influence of noise on the distance is proportional to the energy of the window weight distribution. This influence is minimal when using a flat window or similarly an isotropic kernel. When using the proposed adaptive kernel, the weight of large gradient points is reduced and the weight of non gradient points increased.  This makes the window weighting less uniform. This noise sensitivity is reduced by increasing the window size, as  shown in the same table.

\begin{table}
\centering
\begin{tabular}{|c|ccc|} \hline
7$\times$7 & $\sigma = 0.0$ & $\sigma = 1.0$ (SNR=100) & $\sigma = 2.0$ (SNR=50) \\ \hline
SSD & 0.118 & 0.121 & 0.138 \\
Proposed &0.108 & 0.113 & 0.139\\ \hline
\end{tabular}

\vspace{0.5 cm}

\begin{tabular}{|c|ccc|} \hline
11$\times$11& $\sigma = 0.0$ & $\sigma = 1.0$ (SNR=100) & $\sigma = 2.0$ (SNR=50) \\ \hline
SSD & 0.135 & 0.136 & 0.139\\
Proposed & 0.107 & 0.109& 0.116\\ \hline
\end{tabular}

\caption{Average error on the disparity computation on the six images of Fig.  \ref{fig:base} and the middle ground truth of Fig.  \ref{fig:base_gt}. For the proposed method the distance is computed to the convolved ground truth as predicted by the formulas. The first table is obtained by using a correlation window of 7$\times$7 pixels while the second table is obtained by using a correlation prolate of size 11$\times$11. We observe that the SSD error increases when using a larger window. By using a larger window the ground truth disparity varies more and the possibility of having a large gradient increases, therefore making SSD more sensitive to adhesion. The obtained errors are quite similar for both algorithms, showing that the use of an adaptive SSD does not diminish the precision of SSD.} \label{tb:mse}
\end{table}




\begin{table}
\centering
\begin{tabular}{|c|ccc|} \hline
7$\times$7 & $\sigma = 0.0$ & $\sigma = 1.0$ (SNR=100) & $\sigma = 2.0$ (SNR=50) \\ \hline
SSD & 0.35 & 0.54 &  1.27  \\
Proposed &  0.04 & 0.20 & 1.26 \\ \hline
\end{tabular}

\vspace{0.5 cm}

\begin{tabular}{|c|ccc|} \hline
11$\times$11& $\sigma = 0.0$ & $\sigma = 1.0$ (SNR=100) & $\sigma = 2.0$ (SNR=50) \\ \hline
SSD &  0.48 & 0.50  & 0.64 \\
Proposed & 0.01 & 0.01 & 0.11 \\ \hline
\end{tabular}

\caption{Average on the six images of Fig.  \ref{fig:base} and the middle ground truth of Fig.  \ref{fig:base_gt} of the percentage of points  with an angular difference of the surface normal to the ground truth normal larger than 10 degrees.
For the proposed method the distance is computed to the convolved ground truth as predicted by the formulas. The first table is obtained by using a correlation window of 7$\times$7 pixels while the second table is obtained by using a correlation prolate of size 11$\times$11. Observe that with a larger correlation window   a surface more similar  to the original one is obtained. This result is notable: the obtained percentage of points with a very different normal to the surface is much higher for the classical SSD than the proposed algorithm.} \label{tb:mse2}
\end{table}


\begin{figure}
\centering
\epsfig{file=back.gt.eps, width=2.5 cm}
\epsfig{file=back.left.eps, width=2.5 cm}
\epsfig{file=back.right.eps, width=2.5 cm}
\caption{Synthetic image pair. Left: the disparity ground truth, the background has uniform disparity while the building simulates the slope of a roof. Center and right: image pair. }\label{fig:back_gt}
\end{figure}

\begin{figure}
\centering
\epsfig{file=back.cor.p7.eps, width=2.5 cm}
\epsfig{file=back.cor.p11.eps, width=2.5 cm} \hspace{0.5 cm}
\epsfig{file=back.corb.p7.eps, width=2.5 cm}
\epsfig{file=back.corb.p11.eps, width=2.5 cm}

\epsfig{file=back.difp7.eps, width=2.5 cm}
\epsfig{file=back.difp11.eps, width=2.5 cm} \hspace{0.5 cm}
\epsfig{file=back.difbp7.eps, width=2.5 cm}
\epsfig{file=back.difbp11.eps, width=2.5 cm}

\caption{Obtained disparities for the synthetic image pair in Fig. \ref{fig:back_gt}. The top left columns display the disparities obtained with a classical SSD algorithm with an isotropic weighting window of size 7$\times$7 and 11$\times$11. The top right columns show the same experiments but with the proposed algorithm. Bottom: image difference between the estimated disparities and the ground truth. For the proposed strategy the displayed image difference stands between the estimated disparities and the convolved ground truth by the same prolate. }\label{fig:back_est}
\end{figure}

The next experiment was performed with a synthetic disparity map applied to a building image. The background has uniform disparity but the building has a sloped roof. Since the background has uniform disparity, we can only observe the fattening effect in and near the building. The ground truth disparity and the simulated image pair are shown in Fig.  \ref{fig:back_gt}.   Fig.  \ref{fig:back_est} shows the estimated disparities with the classical SSD algorithm and with the proposed adaptive SSD, using again  prolate windows of 7$\times$7 and 11$\times$11 pixels. The same figure shows the error image, namely the difference  between the estimated disparities and the ground truth. With the proposed strategy the obtained image difference stands between the estimated disparities and the convolved ground truth by the same prolate. This is consistent with the formulation in the previous section, where we showed that the adaptive SSD estimates a convolved disparity, independent of the gradient of the image.  For the SSD algorithm, we observe a prominent error near the boundaries of the building, while for the proposed strategy this error passed unnoticed.

The next  experiment displays a more complicated case with occlusion and shadows containing nearly no information. Fig. \ref{fig:prison_gt} shows the image pair and its ground truth. In Fig. \ref{fig:prison_est} are displayed  the estimated disparities  and  the error image difference  between the estimated disparities and the ground truth. For the proposed strategy the image difference stands again between the estimated disparities and the convolved ground truth by the same prolate. Observe that the  error is mainly concentrated near the edges of the building, where the foreground fattening effect is  severe. Although in the synthetic case of Fig \ref{fig:back_est} we were able to nearly eliminate the error near the edges with the proposed strategy, this is not the case for this pair.
The error committed by the  SSD algorithm is reduced but not eliminated. This is due to the occlusions which make $\varepsilon$  discontinuous, and to the fact that near most of the building boundaries the shadow has removed all possible information that could be used to correct the match. Surprisingly, the error is much smaller at non shadowed edges, even if occlusions and discontinuities of the disparity are still present.


\begin{figure}
\centering
\epsfig{file=prison.gt.eps, width=2.5 cm}
\epsfig{file=prison.left.eps, width=2.5 cm}
\epsfig{file=prison.right.eps, width=2.5 cm}
\caption{Synthetic image pair. Left: the disparity ground truth. The background has uniform disparity while the building has a sloped roof. Center and right: image pair.}\label{fig:prison_gt}
\end{figure}


\begin{figure}
\centering
\epsfig{file=prison.p7.eps, width=2.5 cm}
\epsfig{file=prison.p11.eps, width=2.5 cm} \hspace{1 cm}
\epsfig{file=prisonb.p7.eps, width=2.5 cm}
\epsfig{file=prisonb.p11.eps, width=2.5 cm}


\epsfig{file=prison.difp7.eps, width=2.5 cm}
\epsfig{file=prison.difp11.eps, width=2.5 cm} \hspace{1 cm}
\epsfig{file=prisonb.difp7.eps, width=2.5cm}
\epsfig{file=prisonb.difp11.eps, width=2.5 cm}

a) Classical SSD      \hspace{3 cm}   b) Adaptive SSD

\caption{Estimated disparities with the classical SSD (a) and  adaptive SSD strategy (b) for the synthetic image pair in Fig. \ref{fig:back_gt}. Top:  disparities obtained with a weighting window of size 7$\times$7 and 11$\times$11. Bottom: image difference between the estimated disparities and ground truth.}\label{fig:prison_est}
\end{figure}


\subsection{Comparison with foreground fattening elimination strategies}

As exposed in the introduction, many strategies have been proposed to remove the fattening effect and are beautifully reviewed and compared in \cite{Scharstein02}.
Our goal now is to compare the proposed strategy with two of the more performing algorithms. Yoon et al. \cite{Yoon06}  selects an adaptive window containing only pixels with a  grey level similar to the reference one, in the spirit of bilateral filters  \cite{tomasi1998bilateral}. The main idea is to keep in the correlation window only points belonging to the same object, which are likely to have a similar grey level.  The Fusiello et al. \cite{fusiello2000symmetric} classic {\it min-filter} chooses among all the windows containing the reference pixel the one which has a minimal distance with its corresponding one in the second image.

Fig. \ref{fig:newalgoPrison} compares the adaptive strategy with these two algorithms on the pair of Fig. \ref{fig:prison_gt}.  The three estimated disparities remove the dilatation of buildings due to the fattening effect. But the estimated disparity by RAFA is more blurred than the  other two, since the recovered disparity is by our theorem a convolution of the original one with the correlation window.


\begin{figure}
\centering
\epsfig{file=newalgosPrison/outSubPixelBest.char.eps,width=3.5 cm}
\epsfig{file=newalgosPrison/outSubPixelAdapt.char.eps,width=3.5 cm}
\epsfig{file=newalgosPrison/outSubPixelRafa.char.eps,width=3.5 cm}
\caption{Estimated disparity on stereo pair in Fig. \ref{fig:prison_gt}.  From left to right:   Fusiello et al \cite{fusiello2000symmetric} min-filter, Yoon et al \cite{Yoon06} bilateral strategy, and the proposed adaptive window strategy. The three estimated disparities remove the dilatation of buildings due to the fattening effect. The estimated disparity by RAFA is more blurred than the other two,  because the recovered disparity is a convolution of the original one with the correlation window.   }\label{fig:newalgoPrison}
\end{figure}

%
%
%\begin{figure}
%\centering
%\epsfig{file=newalgosPrison/dome.char.eps,width=3.5 cm}
%\epsfig{file=newalgosPrison/difSubBest.char.eps,width=3.5 cm}
%\epsfig{file=newalgosPrison/difSubAdapt.char.eps,width=3.5 cm}
%\epsfig{file=newalgosPrison/difSubRafa.char.eps,width=3.5 cm}
%\caption{ }\label{fig:}
%\end{figure}


In order to evaluate the subpixel precision of the three methods, we applied the algorithms to the first texture image in Fig. \ref{fig:base} and the second simulated ground truth disparities in Fig. \ref{fig:base_gt}.  Fig. \ref{fig:newalgoTexture} displays the estimated disparity for the three algorithms. It is observed that the disparity estimated  by the min-filter produces a shock effect which creates discontinuities of the estimated disparity.  These shocks are not present when using the adaptive window of Yoon et al.. However, many irregularities are present in the estimated disparity, which are similar to the ones obtained by the classical SSD. The estimated disparity by RAFA algorithm is more similar to the ground truth. This can also be observed by looking at the 3D representation of the estimated disparities by the three algorithms.


\begin{figure}
\centering
\epsfig{file=newalgosText/outSubPixelBest.char.eps,width=3.5 cm}
\epsfig{file=newalgosText/outSubPixelAdapt.char.eps,width=3.5 cm}
\epsfig{file=newalgosText/outSubPixelRafa.char.eps,width=3.5 cm}
\caption{Estimated disparity on the first texture image in Fig. \ref{fig:base}. Independent additive white noise with standard deviation 2 was added to both images before matching.  From left to right: the min-filter  Fusiello et al \cite{fusiello2000symmetric}, the bilateral  Yoon et al \cite{Yoon06}, and the proposed adaptive window strategy. The RMSE (in pixels) are respectively 0.18, 0.19, and 0.11.}\label{fig:newalgoTexture}
\end{figure}

\begin{figure}
\centering
\epsfig{file=newalgosText/outSubPixelBestSurf.eps,width=5.5 cm}
\epsfig{file=newalgosText/outSubPixelAdaptSurf.eps,width=5.5 cm}
\epsfig{file=newalgosText/outSubPixelRafaSurf.eps,width=5.5 cm}
\caption{Three dimensional representation of the estimated disparity from  Fig. \ref{fig:newalgoTexture}. As predicted by the theorem, the fattening effect is optimally removed by the adaptive window. }\label{fig:newalgoTexture3d}
\end{figure}



The conclusion of these comparisons is that classical fattening removal techniques work correctly when for {\it foreground fattening}  due to the presence of important disparity discontinuities.
 Nevertheless, as was pointed out, fattening occurs everywhere, even  in the absence of strong depth discontinuities. In that case classical techniques do not attain the optimal precision obtained with the adaptive window technique.  This might be due to the fact that most existing databases for stereo benchmarks furnish only a  pixel precise ground truth. Therefore the benchmarks do not permit to detect or to evaluate the loss of precision due to the general surface fattening. To demonstrate this fact,  Fig. \ref{fig:middlebury} shows a detail of a stereo pair of the Middlebury dataset \cite{scharstein2002middlebury}. Are shown the estimated disparities by the classical SSD and by the adaptive RAFA.  Both disparities present many irregularities probably present in the object (a rough stone). However, the smoothness of the furnished ground truth did not allow for a numerical comparison of the two estimated disparities.

\begin{figure}
\centering
\epsfig{file=middlebury/rocksl.eps,width=4.5 cm}
\epsfig{file=middlebury/rocksr.eps,width=4.5 cm}
\caption{Stereo pair obtained from Middlebury dataset \cite{scharstein2002middlebury}. }\label{fig:middlebury}
\end{figure}


\begin{figure}
\centering
\epsfig{file=middlebury/crop.eps,width=3 cm}
\epsfig{file=middlebury/sub.char.eps,width=3 cm}

\epsfig{file=middlebury/gsub.char.eps,width=3 cm}
\epsfig{file=middlebury/disp.eps,width=3 cm}

\caption{Obtained disparities by classical SSD and RAFA on a piece of Fig. \ref{fig:middlebury}.  From top to bottom and left to right: detail  in the reference image, classical SSD disparity, adaptive RAFA disparity and ground truth furnished in the same dataset  \cite{scharstein2002middlebury}. The ground truth furnished is too smooth and quantized. Fattening effects cannot be compared on it.}\label{fig:middlebury2}
\end{figure}



\section{Conclusion} This paper has shown that in block matching methods the fattening phenomenon occurs everywhere.  A mathematical analysis has proved that  fattening could  be completely avoided in the regions with smoothly varying disparity,  by introducing adaptive weights in the SSD block matching. Experimental evidence on simulated data has been provided to confirm that fattening is indeed avoided with the adaptive window. Yet the adaptive window does not correct the strong foreground fattening, particularly annoying near large building edges in aerial imaging. However, the adaptive window promises to be a valuable and simple correction the fixed windows used widely in block matching method. Future work will consider how to insert this correction in a complete stereo reconstruction chain.




\bibliographystyle{plain}
\bibliography{biblio_stereo}

\end{document}

