One-view occlusion detection for stereo matching with a fully connected CRF model

TL;DR

This paper introduces a novel one-view occlusion detection method within a fully connected CRF model for stereo matching, improving accuracy and efficiency by reducing the need for dual-view processing.

Contribution

It proposes the OVOD approach for occlusion detection that simplifies stereo matching by performing a single energy minimization, enhancing accuracy over standard methods.

Findings

Significantly improves stereo matching accuracy on Middlebury dataset.

Reduces computational cost by avoiding second view energy calculations.

Achieves state-of-the-art results in median, average, and mean squared error metrics.

Abstract

In this paper, we extend the standard belief propagation (BP) sequential technique proposed in the tree-reweighted sequential method to the fully connected CRF models with the geodesic distance affinity. The proposed method has been applied to the stereo matching problem. Also a new approach to the BP marginal solution is proposed that we call one-view occlusion detection (OVOD). In contrast to the standard winner takes all (WTA) estimation, the proposed OVOD solution allows to find occluded regions in the disparity map and simultaneously improve the matching result. As a result we can perform only one energy minimization process and avoid the cost calculation for the second view and the left-right check procedure. We show that the OVOD approach considerably improves results for cost augmentation and energy minimization techniques in comparison with the standard one-view affinity space implementation. We apply our method to the Middlebury data set and reach state-of-the-art especially for median, average and mean squared error metrics.

Tables

Table 1. TABLE I: Middlebury benchmark V3 for different error metrics. Here we compare our method with five other methods. Snapshot taken from website on December 15, 2017. Only average over all images of the full test data set is provided. The top part of this table with the all pixels mask and the bottom the non-occluded mask.

Error metrics:	bad 0.5	bad 1	bad 2	bad 4	avrg	rms	A50	A90	A95	A99
OVOD (our)	${45.8}_5$	${24.6}_5$	${15.8}_8$	${10.8}_7$	${5.19}_2$	${22.3}_2$	${0.50}_5$	${5.60}_5$	${\text{27.1}}_1$	${\text{108}}_1$
LocalExp [33]	${\text{44.2}}_1$	${\text{21.0}}_1$	${\text{11.7}}_1$	${8.83}_3$	${\text{5.13}}_1$	${\text{21.1}}_1$	${\text{0.47}}_1$	${5.27}_4$	${31.6}_6$	${109}_2$
3DMST [20]	${45.7}_4$	${22.0}_2$	${12.5}_3$	${8.81}_2$	${5.40}_3$	${23.5}_3$	${0.49}_4$	${4.23}_2$	${29.5}_3$	${114}_3$
PMSC [21]	${45.4}_3$	${22.8}_4$	${13.6}_4$	${9.71}_4$	${5.65}_4$	${23.9}_4$	${0.49}_3$	${5.20}_3$	${31.4}_5$	${115}_4$
MeshStereoExt [44]	${47.0}_9$	${24.8}_6$	${15.7}_7$	${12.0}_{10}$	${11.3}_{26}$	${39.0}_{36}$	${0.50}_6$	${26.4}_{20}$	${88.4}_{35}$	${193}_{40}$
NTDE [13]	${46.9}_7$	${24.9}_7$	${15.3}_6$	${10.4}_5$	${6.20}_7$	${26.1}_7$	${0.51}_7$	${5.82}_6$	${34.4}_7$	${128}_7$
OVOD (our)	${40.0}_5$	${17.1}_9$	${8.87}_{13}$	${5.03}_{12}$	${2.23}_2$	${\text{12.9}}_1$	${0.44}_3$	${2.08}_{11}$	${4.80}_6$	${\text{52.6}}_1$
LocalExp [33]	${\text{38.7}}_1$	${\text{13.9}}_1$	${\text{5.43}}_1$	${\text{3.69}}_1$	${2.24}_3$	${13.4}_4$	${\text{0.43}}_1$	${1.55}_3$	${4.81}_7$	${55.3}_3$
3DMST [20]	${39.9}_4$	${14.5}_2$	${5.92}_2$	${3.72}_2$	${\text{2.17}}_1$	${13.3}_3$	${0.44}_3$	${1.53}_2$	${\text{3.47}}_1$	${55.0}_2$
PMSC [21]	${39.1}_2$	${14.8}_3$	${6.71}_4$	${4.44}_6$	${2.26}_4$	${\text{12.9}}_1$	${\text{0.43}}_1$	${1.90}_8$	${5.27}_8$	${56.0}_4$
MeshStereoExt [44]	${40.1}_7$	${15.6}_5$	${7.08}_6$	${4.36}_5$	${2.49}_5$	${15.4}_7$	${0.44}_3$	${1.71}_5$	${4.05}_4$	${67.3}_7$
NTDE [13]	${40.4}_9$	${16.2}_8$	${7.44}_9$	${4.55}_8$	${2.58}_6$	${15.7}_{10}$	${0.44}_3$	${1.77}_6$	${4.74}_5$	${65.4}_6$

Equations

$$\mathrm{E}\left(s\right)=\sum _{p\in \mathcal{V}}{c}_p\left(s_p\right)+\sum _{\left(p,q\right)\in \mathcal{E}}{B}_{p,q}\left(s_p,s_q\right),$$

$$B_{p,q}\left(s_p,s_q\right)=b_q{w}_{p,q}\phi \left(s_p,s_q\right),$$

$$b_q={\left({\sum }_{p\in \mathcal{V}}{w}_{p,q}\right)}^{-1}$$

$$s_q=\underset{s\in S}{\mathrm{argmin}}{\bar{c}}_q\left(s\right)$$

$${\bar{c}}_q^t\left(s\right)=c_q\left(s\right)+\sum _{p,q\in \mathcal{E}}{m}_{p\to q}^{t-1}\left(s\right),$$

$${\bar{c}}_q^t\left(s\right)=c_q\left(s\right)+\sum _{p|p,q\in \mathcal{E}}{m}_{p\to q}^{t-1}\left(s\right)$$

$$m_{p\to q}^t\left(i\right)=\underset{j\in S}{\min }\left({\beta }_{p,q}{\bar{c}}_p^t\left(j\right)-m_{p\to q}^{t-1}\left(j\right)+B_{p,q}\left(i,j\right)\right),$$

$$\begin{array}[]{l}m_{p\to q}^{t}\left(i\right)={b_{q}}{w_{p,q}}m_{p}^{t}\left(i\right),\\ m_{p}^{t}\left(i\right)=\mathop{\min}\limits_{j\in S}\left({\bar{c}_{p}^{t}\left(j\right)+{\varphi_{p,q}}\left({i,j}\right)}\right)\\ \end{array}$$

$${\bar{c}}_q^t\left(s\right)=c_q\left(s\right)+\frac{b_q}{1-b_q}\sum _{p\in \mathcal{V}}{w}_{p,q}{m}_p^{t-1}\left(s\right)-m_q^{t-1}\left(s\right)$$

$$\begin{array}[]{*{20}{c}}{{{\tilde{f}}_{q}}=\frac{1}{{{W_{q}}}}\sum\limits_{p\in{\cal V}}{{e^{-a{d_{p,q}}}}}{f_{p}}}\\ {{W_{q}}=\sum\limits_{p\in{\cal V}}{{e^{-a{d_{p,q}}}}}},\\ \end{array}$$

$$\begin{array}[]{l}{d_{q,p}}=\mathop{\min}\limits_{{P_{p,q}}}\;\sum\limits_{\varepsilon\in{P_{p,q}}}{{u_{\varepsilon}}},\\ {u_{\varepsilon=\left({q,p}\right)}}=\left\|{{I_{q}}-{I_{p}}}\right\|+\delta,\\ \end{array}$$

$$a=\frac{2}{{\sigma }_r^2},\delta =\frac{{\sigma }_r^2}{{\sigma }_s^2},$$

$$w_{p,q}=e^{\underset{P_{p,q}}{\min }\sum _{\epsilon \in P_{p,q}}{u}_{\epsilon }}=\underset{P_{q,p}}{\max }\prod _{\epsilon \in P_{p,q}}{e}^{-a{u}_{\epsilon }}$$

$${\bar{c}}_q^t\left(s\right)=c_q\left(s\right)+\frac{1}{W_q-1}\sum _{p\in \mathcal{V}}{w}_{p,q}{m}_p^t\left(s\right)-m_q^t\left(s\right)$$

$${o_{{q^{L}}}}=\left\{\begin{array}[]{l}\begin{array}[]{*{20}{c}}{occ}&{{\text{if }}{s^{L}}\left({x_{q}^{L},{y_{q}^{L}}}\right)\neq{s^{R}}\left({x_{q}^{L}-s_{q}^{L},{y_{q}^{L}}}\right)}\\ \end{array}\\ \begin{array}[]{*{20}{c}}{nocc}&{{\text{elsewhere}}}\\ \end{array}\\ \end{array}\right.$$

$$x_q^R=x_q^L-s$$

$$s_q^R=\underset{s\in S}{\mathrm{argmin}}\bar{c}\left(x_q^L-s,y_q,s\right)$$

$$s_q=\min \left(s_q^L,s_q^R\right)$$

$$c_p=\min \left(1.5c_p^{CNN},\sqrt{c_p^{CNN}}\right)$$

$$\phi \left(s_p,s_q\right)=\lambda \min \left(\left|s_p-s_q\right|,{\Delta }_{\max }\right),$$

$${w^{RC}}\left({z,p}\right)=\left\{\begin{array}[]{l}\begin{array}[]{*{20}{c}}0&{{\rm{if}}\left\|{zM-p}\right\|>1}\\ \end{array}\\ \begin{array}[]{*{20}{c}}{\frac{1}{2}\left({1+\cos\left({\frac{{\pi\left\|{zM-p}\right\|}}{{M}}}\right)}\right)}&{{\rm{else}}}\\ \end{array},\\ \end{array}\right.$$

$$\begin{array}[]{l}{s^{in}}\left({x,y}\right)=\left\{\begin{array}[]{l}\begin{array}[]{*{20}{c}}{{s^{dwsc}}\left({\xi,\zeta}\right)}&\begin{array}[]{l}x=M\xi\\ y=M\zeta\\ \end{array}\\ \end{array}\\ \begin{array}[]{*{20}{c}}0&\text{elsewhere}\\ \end{array}\\ \end{array}\right.\\ ms{k^{in}}\left({x,y}\right)=\left\{\begin{array}[]{l}\begin{array}[]{*{20}{c}}1&\begin{array}[]{l}x=M\xi\\ y=M\zeta\\ \end{array}\\ \end{array}\\ \begin{array}[]{*{20}{c}}0&\text{elsewhere}\\ \end{array}\\ \end{array}\right.\\ \end{array}$$

$$\begin{array}[]{l}{s^{out}}=\frac{1}{{{W^{out}}}}{F^{GD}}\left\{{{s^{in}}|{I^{L}},{\sigma_{ru}},{\sigma_{su}}}\right\}\\ {W^{out}}={F^{GD}}\left\{{ms{k^{in}}|{I^{L}},{\sigma_{ru}},{\sigma_{su}}}\right\},\\ \end{array}$$

$$h_q\left(i\right)=\sum _{p\in N_{\rho }}{e}^{-\frac{{|p-q|}^2}{2{\rho }^2}-\frac{{|I_p-I_q|}^2}{2{\sigma }_m^2}}{\delta }_K\left(s_p,i\right)s_p,$$

$$\underset{i\in S}{s_q=\arg \mathrm{med}}\left(h_q\left(i\right)\right)$$

$${\tilde{s}}_q=\sum _{p\in N_{\rho }}{e}^{-\frac{{|{\pi }_p-s_q|}^2}{2{\sigma }_{\pi }^2}-\frac{{|\nabla s_p-\nabla s_q|}^2}{2{\sigma }_g^2}}{\pi }_p\left(\nabla s_p,s_p,q\right),$$

Full text

I Introduction

Stereo matching is considered to be one of the fundamental problems of computer vision  [4, 30]. It is important for a large variety of computer vision applications, such a intermediate view generation, 3D scene reconstruction, autonomous driving systems and robotics.

Methods that solve the stereo matching problem can be divided into two categories: cost filtering methods and energy minimization methods. The earliest cost filter approaches use fixed or adaptive windows as the matching neighborhood. Later they have been extended into algorithms which operate in the bilateral color-image space [20, 22, 23, 34, 38, 41]. An efficient alternative to bilateral cost filtering is proposed by Rhemann et al. [28], which outperforms the earliest bilateral filter in speed and quality. This algorithm is based on the guided filter approach [10].

A common way to solve the stereo matching problem with energy minimization is by performing inference in Bayesian networks. The network should converge to the maximum a posteriori (MAP) probability or to the energy minimum of the conditional random field (CRF) functional defined over the image graph, where pixels or image patches (superpixels) are vertices and pairwise relationship between pixels are usually encoded with graph edges [17, 11]. Early work considers pairwise potentials of the CRF energy functional in restricted neighborhoods (locally connected models)  [3, 16, 31, 12, 32]. This constraint makes the probabilistic graphical model less effective and flexible.

The matching cost formation is an important component of any stereo matching algorithm. Traditionally, costs can be divided in two groups. The first group is based on per-pixel matching dissimilarity measures [30]. The second group is based on non-parametric transforms such as rank and census [42] or normalized cross correlation [19]. Recent progress in the field of convolutional neural networks (CNNs) has yielded more robust matching costs for stereo [43]. Consequently, the traditional cost computation has been replaced by the CNN based cost in most recent works. Also in this work cost computation is based on CNNs [33].

I-A Contributions

In this paper we propose two main contributions: a) we derive the energy minimization solution for the fully connected CRF model in the framework of BP [32] instead of mean field (MF), as used in [18]. We adapt the sequential technique proposed in [15] to the fully connected CRF model, and show that the sequential algorithm for the fully connected model accelerates convergence when compared to the non-sequential version of the algorithm; b) we propose the one-view occlusion detection (OVOD) approach that provides a new marginal solution for BP. In contrast to the standard WTA estimation the proposed OVOD solution allows to find occluded regions in the disparity map and simultaneously improves the matching result in these regions. As a consequence, we are only required to perform one energy minimization process, we can omit the cost calculation step for the second view and avoid the left-right check procedure. To demonstrate the general ability of our algorithm we apply our method to the Middlebury data set and reach state-of-the-art level for median, average and mean squared error metrics.

Note that the method proposed in [15] considers only a locally-connected CRF model. In [24] we proposed a recursive filter that is based on the geodesic affinity. However, this version does not apply to the energy minimization problem. Here we use that approach as a tool that allows us to extend the locally-connected CRF technique in [15] to the fully-connected CRF model. Furthermore, the proposed extension of [15] is theoretically proven and can be used not only in stereo but also in any other problem, which is based on energy minimization with a fully-connected CRF model.

I-B Related work

Generally cost aggregation and global optimization are considered as different techniques. Nevertheless, in two state-of-the-art works authors have proposed methods to combine local preprocessing within a global optimization framework. They combine segmentation based cost aggregation with BP energy minimization [14, 35]. Firstly, an initial image is divided into relatively small segments also called superpixels. Then an aggregated cost is used for the unary potential in the energy minimization problem. There are two problems that one faces implementing this approach: firstly, segmentation is not a trivial problem especially in context of stereo matching; secondly, the assumption that pixels inside the same segment belong to similar disparity values can possibly produce estimation errors, which are impossible to correct at the stage of disparity map post-processing.

An extension of locally connected CRF models to fully-connected CRF models was proposed in [25]. However in that method only one iteration is used and this iteration is equivalent to the bilateral cost filtering. An iterative solution of fully connected CRF models was proposed in [2], but there the energy optimization was done by a Quasi-Newton (L-BFGS) method.

In the segmentation paper of Krahenbuhl and Koltun [18] fully connected CRF models have been proposed. They solve the minimization problem by MF approximations. However, comparison with belief propagation methods show that BP typically performs much better than MF [37]. Thus, we propose to use the BP approach developed in the paper by Kolmogorov [15], where a sequential calculation scheme for BP convergence was developed. This method significantly reduces the number of iterations needed to reach the same level of energy compared to its non-sequential analogue proposed in [36].

The algorithms [18, 2] use the color-space bilateral distance to compute the pairwise potential in the CRF energy functional. In the presented work the same binary potential is constructed with the geodesic distance affinity [45, 5]. The motivation of this choice is twofold. Firstly, fast realization of the bilateral filter is necessary to solve the CRF energy minimization problem for the fully connected model in a reasonable computational time. However, the fast bilateral filter algorithm [1] that is used in [18, 2] is in general slower than the recursive filtering based on the geodesic distance [8, 39], and for several kernel intrinsic parameters is considerably slower. Secondly and more important, all known fast algorithms of the bilateral filter [27, 1, 9, 26] do not allow sequential update of the belief propagation messages on the image graph tree. The last observation is also relevant to the algorithms based on the guided filter approach [10]. Formally, affinity weights in the guided filter are calculated recursively, however such a recursion could not be adapted to the sequential process using the standard scheme in [10].

To realize the fast geodesic distance based filter, which is a part of our CRF energy minimization algorithm we follow the idea in [5] and especially in [24], where recursive summation trees are proposed. Note that the algorithm in [24] was developed for denoising filters. We extend this technique and adapt the recursive trees to the sequential calculation scheme for BP message passing.

In [38] a non-local cost aggregation method for stereo matching was proposed, and the proposed algorithm uses the geodesic affinity for the cost filtering approach. In contrast, we use the geodesic distance not for the cost filtering, but for a message update procedure with a fully-connected CRF model.

Now we will briefly discuss some of the most recent advances in the field of accurate stereo matching. Two recently proposed methods, the Local Expansion Moves [33] and PatchMatch-based superpixel cut [21], solve the energy minimization labeling problem with a locally connected graph cut technique. The methods construct the energy smoothness term that fits the solution to oriented planes, while traditional methods consider only front-parallel planes. Such a smoothness term results in increased accuracy, but increases the dimensionality of the label space from 1D to 3D, consequently increasing the general computational complexity of the methods. We discuss this further in the experimental section of our paper.

In the paper [20], the authors propose a cost-aggregation method that embeds a minimum spanning tree based support region filtering in the general procedure. Furthermore, it combines this approach with a PatchMatch 3D label search, thus performing the search process with an adaptive patch size. An efficient stereo matching algorithm, which applies adaptive smoothness constraints using texture and edge information, is proposed in [13]. The approach includes two base procedures: the first one segments the disparity map domain into three classes with different smoothness terms, and the second procedure obtains disparity values by minimizing a local cost volume energy. Finally, a 3D mesh technique is proposed in [44]. This novel global stereo model was initially designed for view interpolation. The method estimates a disparity map and a 3D triangular mesh in such a way that both problems are solved simultaneously.

II Energy minimization with fully connected CRF model

In this section we consider energy minimization with fully connected CRF models. We extend the locally connected model proposed in [15] to the fully connected one. Also we join the two step energy minimization approach [25] (which consists of the bilateral cost filtering and the locally connected energy minimization algorithm [15]) into one iterative energy minimization process of a fully connected model with the geodesic distance affinity.

The principal difference between fully connected models and locally connected models is the number of pixels that are included in the considered pixel neighborhood. In fully connected models all image pixels are potentially connected to each other. However the strength of each connection depends on a chosen affinity metric between pixels. This connectivity strength is formalized by the reciprocal image affinity weights. Consequently, these weights become a part of binary potentials in the smoothness term of the energy functional. First, let us consider the energy function of the following form:

[TABLE]

where (and further in this paper) $p,q,k,l\in{\cal V}$ corresponds to pixels or vertices and set $\varepsilon=\left({p,q}\right)\in{\cal E}$ to edges of an image graph ${\cal G}=\left\{{{\cal V},{\cal E}}\right\}$. Where ${s_{p}}$ denotes the label of pixel $p$ that belongs to some discrete set of disparities $s_{p}\in S$; ${c_{p}}\left(\cdot\right)$ defines a unary potential which corresponds to the conventional penalty cost or the logarithm of the likelihood probability. ${B_{p,q}}\left({\cdot,\cdot}\right)$ is a binary potential which defines edge interaction between pixels $\left({p,q}\right)$. In the case of the fully connected CRF model each combination of pixels $\left({p,q}\right)$ of the image graph ${\cal G}$ belongs to the edge set ${\cal E}$. And the binary potential for the model can be expressed as:

[TABLE]

where $\varphi\left({\cdot,\cdot}\right)$ is a pairwise penalty that enforces smoothness of the estimated label function, and ${w_{p,q}}$ is a pairwise weight that defines some reciprocal influence between every node. The variable ${b_{q}}$ is a normalization factor set to

[TABLE]

The problem of energy minimization in Eq. (1) can be solved in the framework of BP. The aim is to find the marginal function ${\bar{c}_{q}}\left(\cdot\right)$ and then estimate labels in each pixel $q$ by the simple WTA formula, which finds the maximum likelihood:

[TABLE]

The core operation of BP is passing a message ${m_{p\to q}}$ from node $p$ to node $q$ for the directed edge $\left({p\to q}\right)$, where $\left({p,q}\right)\in{{\cal E}}$. The marginal function is then expressed as follows:

[TABLE]

or the following

[TABLE]

Here $t$ is the number of the iteration, and it is supposed that the related iterative message passing procedure is given [32]:

[TABLE]

where indexes ${i,j}\in S$ belong to the discrete disparity domain $S$, and ${\beta_{p,q}}$ is a normalization factor that is set to ${{\beta_{p,q}}={b_{q}}{w_{p,q}}}$. The term $m_{p\to q}^{t-1}\left(j\right)$ can be neglected for the fully connected model, because in general this value vanishes relative to the sum of all other messages incoming to pixel $q$. Thus we can rewrite Eq. (7) as follows:

[TABLE]

The marginal of Eq. (5) can also be rewritten:

[TABLE]

One can see that the sum term ${\sum\nolimits_{p\in{\cal V}}{{w_{p,q}}m_{p}^{t-1}}}$ in Eq. (9) coincides with the discrete image filtering definition for a large class of filter kernels ${{w_{p,q}}}$, which include traditional Gaussian, bilateral or geodesic distance based kernels. In this paper we propose to use the fast geodesic distance based filter described in the next Section, because its recursive graph-tree based calculation allows sequential updating of messages using Eq. (8). Thus the solution of the energy minimization problem with the CRF fully connected model can be obtained by an iterative process using Eq. (9), Eq. (8), and Eq. (4).

III Geodesic distance based recursive calculation tree

As pointed out in Section II we propose to use the fast geodesic distance, because its recursive calculation allows sequential updating of messages. Consequently this leads to faster and more accurate convergence of the energy minimization process in contrast to the bilateral or the guided filter. In addition, the topology of the standard bilateral filter with Gaussian kernel does not fit well to the stereo matching problem. This problem arises because disparity values of two different but near pixels can considerably differ while belonging to the same color of a cover image. In this case, the bilateral filter blurs both output values. This ability of the bilateral filter to communicate over color edges could be useful for example for color based segmentation problem, but can produce estimation artifacts for stereo. Consequently, we found that edge preserving filters based on the geodesic distance measure are more suitable to the stereo matching problem.

The geodesic distance based filter is usually chosen in the following form, which makes the filter recursive

[TABLE]

where$f_{p}$ and ${{{\tilde{f}}_{q}}}$ are the input and the output of the filter respectively, the weight $e^{-a{d_{p,q}}}$ defines a geodesic distance based affinity between any two image pixels.

The variable ${d_{p,q}}$ in Eq. (10) is the geodesic distance between image pixels $\left({p,q}\right)$ which for an image ${I_{p}}$ can be defined on the discrete grid graph as

[TABLE]

where ${P_{q,p}}$ is a path between two graph vertices $\left({p,q}\right)$ and $\delta$ is the space term. Remember that in graph theory, a path in a graph is a sequence of edges which connect a sequence of vertices.

Note that the filter intrinsic parameters $a$ and $\delta$ in Eq. (11) approximately correspond to the parameters of the classic bilateral filter with the Gaussian kernel as follows

[TABLE]

where $\sigma_{r}^{2}$ and $\sigma_{s}^{2}$ are the range and the space variance respectively. Also from Eq. (11) one can derive the following useful relation

[TABLE]

In [24] it is shown that the filter sum in Eq. (9) can be calculated efficiently and this calculation scheme allows for sequential message updates.

The full tree, which corresponds to this fast algorithm [24], is composed by four quadrant-domains (or branches of the tree). One of four branches is illustrated in Fig. 1. Here $W_{q}$ represents the sum of weights $w_{p,q}$ for the first quadrant, when ${{\tilde{f}}_{q}}$ gives filter results also for this quadrant. Note that the optimal tree and all weights in Eq. (13) can be calculated for all other procedures.

III-A Message sequential update

The difference between sequential and iterative BP message updates is defined in  [15]. According to this definition one can rewrite the marginal Eq. (9) as

[TABLE]

Taking into account Eq. (5) one can say that a message coming from pixel $p$ depends on all incoming messages of the same iteration. However, because the sequential process assumes that the message calculation comes recursively on the tree branch (quadrant-domain) from leaves (all previous pixels of the branch) to root $p$, only one branch update matters and the full iteration is composed of all four branch passes.

To provide more intuition on the sequential idea of BP message updating, let us consider the one dimensional case (e.g. a string in an image) of BP convergence. It is known that the energy converges after finite number of iterations and that this number is proportional or equal to the number of nodes of which the graph of this one dimensional data consists [7]. On the other hand, a dynamic programming (DP) approach can achieve the solution at once, which is equivalent to one iteration. Note that a recursive transition function from vertex to vertex of DP in this case is exactly equivalent to the BP message coming from the same vertices. Thus the proposed sequential approach applies this dynamic programming principle to the belief propagation technique.

IV One-view occlusion detection concept

Before we explain the idea of one-view occlusion detection (OVOD), we briefly describe the base approach that detects occlusion in stereo and that is used in almost all state-of-the-art algorithms. It is the so called left-right consistency check [6]. This technique involves two independent cost volume calculations. Our method allows to use only one cost volume and therefore only one energy minimization procedure is necessary.

The left-right consistency check approach supposes that firstly two independent cost volumes $c^{L}(q,s)$ and $c^{R}(q,s)$ have to be calculated, where pixels ${q^{L}}$ and ${q^{R}}$ can be represented as two corresponding points with coordinates in the left ${q^{L}}={\left[{x_{q}^{L},y_{q}^{L}}\right]^{T}}$ and the right ${q^{R}}={\left[{x_{q}^{R},y_{q}^{R}}\right]^{T}}$ image planes respectively. Note that here ${s^{L}}$ and ${s^{R}}$ are disparity coordinates of the left and the right view cost volumes respectively. Then both cost volumes are usually processed by energy minimization or cost filtering algorithms to obtain marginal $\bar{c}^{L}(q,s)$ and $\bar{c}^{R}(q,s)$ . Consequently, two independent solutions $s_{q}^{L}$ and $s_{q}^{R}$ can be obtained by using Eq. (4).

Finally, consistency check formula is used to detect occluded pixels:

[TABLE]

where ${o_{{q^{L}}}}$ is the consistency check output occlusion map. If a value ${o_{{q^{L}}}}=occ$ it means that an occluded pixel is detected, and the previously estimated value $s_{q}^{L}$ is uncertain and should be replaced by any gap-filling post-processing procedure.

The idea of OVOD is to avoid the double-cost volume procedure and to omit the gap-filling procedure. It is based on the observation that any visible 3D world point ${q^{W}}$ can be mapped into both cost volume coordinate systems by ${q^{W}}\Rightarrow\left[{x_{q}^{L},y_{q}^{L},s}\right]$ and ${q^{W}}\Rightarrow\left[{x_{q}^{R},y_{q}^{R},s}\right]$ with clear relationship between both mapping: $x_{q}^{L}-x_{q}^{R}=s$ and $y_{q}^{L}=y_{q}^{R}$. A cost volume value ${e^{-c\left({{q^{W}}}\right)}}$ defines a probability of the world point ${q^{W}}$ to belong to a point on the surface of the object. The point ${q^{W}}$ can be mapped into both left and right cost volumes with the same cost value $c\left({{q^{W}}}\right)$. Thus, one can assume that it is possible to use the left-view cost volume as the right-view cost volume after a simple shearing transformation that changes only the X-coordinate:

[TABLE]

Then the left-view marginal volume $\bar{c}^{L}(q,s)$ also can be considered as the right-view marginal volume $\bar{c}^{R}(q,s)$ but with the transformed X-coordinate given by Eq. (16). In this case the right-view solution $s_{q}^{R}$ is obtained by modified Eq. (4):

[TABLE]

As is illustrated in Fig. 2 a marginal solution $s_{q}^{L}$ for a fixed pixel ${q^{L}}$ is found at the minimum in the vertical direction (L) of the unified marginal volume. Another solution is found at the minimum in the slanted direction (R) defined in Eq. (17) and illustrated in Fig. 2. To give an intuition: one can compare the (L) and the (R) directions with light rays coming from a cost volume mapped world point ${p^{W}}$ (see Fig. 2) to the left and to the right lenses of a stereo system respectively. If the value of ${o_{{q^{L}}}}=nocc$ it means that both WTA minimization directions (L) in Eq. (4) and (R) in Eq. (17) should intersect in the solution point $s_{q}^{L}=s_{q}^{R}$ . Therefore, when $s_{q}^{L}=s_{q}^{R}$ this means that an occluded pixel is detected, or the world point ${q^{W}}$ does not belong to the reflected surface of any object. In stereo matching the occluded pixel should be treated. We use the following simple occlusion handling to avoid a gap-filling procedure:

[TABLE]

The idea behind this choice is simple and illustrated in Fig. 2: if pixel $x_{q}^{L}$ is occluded then the neighboring pixel to the right in Fig. 2 is also occluded or its disparity value is larger than $s_{q}^{L}$. Note that we cannot use larger values to fill the gap due to the occlusion definition: disparity of an occluded pixel should always be less than the disparity of pixels that occlude. Thus, to replace the uncertain value $s_{q}^{L}$ we have to choose the disparity value from the left neighborhood pixels and Eq. (18) satisfies the proposed assumption (see Fig. 2). However, it is possible to use other algorithms to fill the occluded regions, because in our algorithm the occlusion map ${o_{{q^{L}}}}$ is calculated as a by-product of the main procedure.

V Stereo matching algorithm

In this section we describe the application of our theory for the specific stereo matching problem.

V-A Energy potentials

Firstly, one has to define the main functional in Eq. (1) for the stereo matching specific problem. The unary potential Eq. (1) or the cost variable ${c_{p}}\left({{s_{p}}}\right)$ in our experiments is directly incorporated from the state-of-the-art CNN-based matching cost function by [33]. This cost is computed by the approach proposed in [43]. The CNN-cost calculation function $c_{p}^{CNN}\left({p,p-s}\right)$ computes a matching cost between left and right image patches of size $11\times 11$ centered at $p$ and $p-s$, respectively, using the MC-CNN neural network [43]. The cost values relevant to the Middlebury data set [29] (half-size) are directly taken from the web site of the project [33] (see relevant reference in the Middlebury evaluation table). Note that in our experiments we use only the left-to-right cost. We use truncated version of this cost $c_{p}^{CNN}$:

[TABLE]

The binary potential ${B_{p,q}}\left({\cdot,\cdot}\right)$ which is redefined via $\varphi\left({{s_{p}},{s_{q}}}\right)$ in Eq. (2) is a truncated linear function:

[TABLE]

where $\lambda$ in our experiments is set to 0.12 and the base value of $\Delta_{\max}$ is 8. For an additional experiment to compare our approach with similar method [25] we use the cost calculating algorithm taken from that paper.

V-B Image-cost downscaling

To make our experiments tractable on a CPU computer we perform a downscale procedure for the cost and stereo images, which are used to form the CRF model affinity space. In general we use different scale factor for the different size of the input image. The base values of the scale factors in our experiments are: $M=4$. There are exceptions for the stereo images with the disparity range $S>300$: $M=5$ and for images with full-size less than $1000\times 500$: $M=3$. Downscaling is performed with the raised cosine kernel, which weight function is:

[TABLE]

where $z$ is a pixel on the downscaled image grid, and $p$ is a pixel on the original grid.

For the cost volume we also perform downscaling in the label-disparity domain with the factor $M$, but in this case we take the min-value of the cost inside the scale interval and remember the real disparity value of this minimum to recover disparity values with initial cost accuracy.

V-C Fully connected model parameters

To solve the functional in Eq. (1) and thus solve the main part of our stereo matching algorithm one needs to perform the message sequential update process defined in Eq. (14). To do that we need to define intrinsic metric weights $w_{p,q}$ of our fully connected model. These weights $w_{p,q}$ are calculated with Eqs. (11-12) and in our experiments $\sigma_{r}=30$ and $\sigma_{s}=8$ are the range and the space standard deviation of the model affinity space respectively.

V-D Disparity map upscaling

Because the output results ${s^{out}}\left({x,y}\right)$ should be the same size as an input image it is necessary to upscale the disparity map ${s^{dwsc}}\left({\xi,\zeta}\right)$ of the energy minimization on the downscaled grid $\left({\xi,\zeta}\right)$ of the main process. For this purpose we propose to use the geodesic distance filter with further notation ${F^{GD}}\left\{\cdot\right\}$ that described in Eq. (10). Firstly, let us define the space functions $s^{in}$ and $msk^{in}$ that are input for the the filter ${F^{GD}}\left\{\cdot\right\}$:

[TABLE]

Note that here the variables $\xi$ and $\zeta$ belong to the integer values.

Then we apply filter to the input functions $s^{in}$ and $msk^{in}$ to obtain the output results ${s^{out}}\left({x,y}\right)$:

[TABLE]

where $I^{L}$ is the left stereo image that defines the geodesic distance affinity of the filter in Eq. (10), and ${\sigma_{ru}}$ and ${\sigma_{su}}$ are intrinsic parameters of the filter relevant to Eq. (12). In our experiments ${\sigma_{ru}=0.85M}$ and ${\sigma_{su}=0.85M}$.

V-E Disparity map post-processing

To obtain state-of-the-art results on stereo estimation several post-processing steps are required. We briefly summarize them here.

Weighted median filter   Note that this filtering could be done before upscaling. Weighted median filter usually is a robust extension of the bilateral filter and widely used in stereo matching post-processing, for example in [40]. The main idea is to accumulate a weighted histogram ${h_{q}}\left(i\in S\right)$ in every pixel based on a previously estimated disparity map:

[TABLE]

where ${{s_{p}}}$ is a value of the disparity map obtained in the previous step, $\delta_{K}$ is the Kronecker delta function, $N_{\rho}$ is a squared neighborhood window of size $\left(2\rho+1\right)\times\left(2\rho+1\right)$, and ${\sigma_{m}}$ and ${\rho}$ are two intrinsic parameters of the weighted median filter. We take the value of ${\rho=3}$ and ${\sigma_{m}=10}$ for not occluded region, and ${\rho=8}$ and ${\sigma_{m}=10}$ for occluded. Then the new value of the desired disparity map is estimated as follows:

[TABLE]

Plane fitting filter   The filter exploits the assumption that the desired disparity function $s_{q}$ belongs locally to the geometric plane with constant gradient ${\nabla{s_{q}}}$. The idea of this filter is partially taken from [33], and we define it as:

[TABLE]

where ${{\tilde{s}}_{q}}$ is the output of the filter and ${\pi_{p}}=\pi_{p}\left({\nabla{s_{p}},{s_{p}},q}\right)$ is a simple predictor function, which can predict value of $s_{q}$ with a known value of the gradient $\nabla{s_{p}}$ and the value ${s_{p}}$ at the pixel $q$. The intrinsic parameter ${\sigma_{\pi}}$ defines the affinity between the predicted value of ${\pi_{p}}$ and the given value $s_{q}$ and it is equal to 5; another intrinsic parameter of the filter $\sigma_{g}=0.2$ enforces the weight $\exp\left({-{{\left|{\nabla{s_{p}}-\nabla{s_{q}}}\right|}^{2}}/2\sigma_{g}^{2}}\right)$ to be zero if the gradients of the desired function $s$ in pixels $q$ and $p$ are considerably different, and thus belong to the different planes; ${{N_{\rho}}}$ is the local neighborhood window of size $15\times 15$.

Disparity sub-pixel enhancement   For this post-processing step we use the approach realized in [40]. The method interpolates the initial cost near the base disparity map with a quadratic polynomial, and then finds the local minimum of the disparity in each pixel.

Small area outliers suppression   This post-processing step approach was realized in [25]. The idea is to detect outliers by segmenting the output disparity map and then finding regions with relatively small area. Consequently, the detected regions are filled by neighborhood values.

Note that the choice of post-processing methods and full pipeline intrinsic parameters of our algorithm is motivated by the output accuracy results based on the Middlebury v3 training data set processing. Our method post-processing steps are applied in in the same order as it is described in this Section. Our readers can find more technical details inside our C++ project for Microsoft Visual Studio 2017 that is referred in the Middlebury evaluation table.

VI Experiments

The experiments have been designed to verify the main contributions of the paper. They are divided into four parts where:

• •

we compare the proposed approach to state-of-the-art;

• •

we evaluate the advantage of the proposed sequential approach over non-sequential approach;

• •

we analyze the advantage of using the OVOD model for energy minimization;

• •

we analyze the results of the proposed fully connected CRF model based on geodesic distance affinity implementation applied to the stereo matching problem.

In our experiments we use the Middlebury data set [29] (half-size).

VI-A Comparison to state-of-the-art methods

Here we compare our method against the state-of-the-art methods on the Middlebury benchmark dataset. We consider the top five methods from the Middlebury benchmark table: LocalExp [33], 3DMST [20], PMSC [21], MeshStereoExt [44], NTDE [13]. Our algorithm achieves similar results as state-of-the-art for average (avrg), root-mean-squared (rms), A90, and A95 error metrics. Results are provided in Fig. 3. In this root-mean-squared error metric our method achieves second rank for the all pixels mask and first-second rank for the non-occluded pixels mask. Also results of 7 images with different masks for this metric achieve first rank in comparison with the second method LocalExp [33] method that reaches only 3 first positions.

We also make Table I to show the position of our algorithm among other top methods. Here the average errors over all the Middlebury data set images for all metrics are presented. From this table we can see that for bad pixels metrics the results of our algorithm is not so good as for other error metrics, but still in the top ten methods for the all pixels mask. We think, that this is the effect of our scaling-up-down procedure. Our main energy minimization process is severely restricted by the factor $M^{3}$ of the cost volume. To apply our method on the original image grid we have to increase memory usage and calculation time up to 125. Thus there is a trade-off between accuracy and computational resources. We also want to point out that the average best algorithm, LocalExp [33], needs to compute the additional right-to-left cost and the corresponding disparity map. In contrast our OVOD approach uses only the left-to-right matching cost and we do not need to estimate the right-to-left disparity map. On the other hand we obtain top position for the A95 and A99 error metrics. These metrics show that our method is very accurate for the majority of the pixels; the average error over the 95% (and 99%) of pixels with least error is very low.

To complete our comparison, we prepare Fig. 4 for the run-time of the methods. One can see that our method is on the first position while all other methods are in inverse order in comparison with the results that are shown in Fig. 3, where the accuracy is considered. Thus, we can see that in general methods that are more accurate need more time to complete the same estimation task. Considering the relative fastness of our method we can say that there are two reasons behind this fact: the OVOD approach uses only one cost volume, thereby reducing its run-time by almost a factor two. Furthermore, our method uses a cost-volume downscaling procedure that considerably decreases computational complexity. For example, if our algorithm works with the original image scale, the largest cost volume consists of 500 mega voxels. It means that we need approximately 50 GB RAM on double floating-point precision. Here we count that every voxel corresponds to at least 8 message values. With our base scale factor that equals to 4 we can reduce the memory usage and the algorithm run-time by approximately 64 times.

We prepare Fig. 5 for visual comparison. In this figure we show the resultant disparity maps of Adirondack, Motorcycle, PianoL, Recycle, Teddy obtained by the top stereo methods.

VI-B Comparison of sequential and non-sequential approaches

For the next experiment on stereo matching we use the single Piano stereo image, because the average energy values over the full data set are not very representative due to different size, the number of disparity labels of an image and thus different energy levels for each stereo image. We have verified that the convergence curves for the other images are almost identical to the Piano stereo image case. As cost we use the cost computed by [25]. In Fig. 6 the result of this experiment is provided. Here the energy value is calculated with Eq.  (1) using intermediate disparity solutions on each iteration. Upper bound values are calculated by methodology proposed in [15]. Pixel accuracy improvement is presented in Fig. 6(b). We can see that the sequential approach yields faster convergence and the gap between energy values and the lower bound estimation is less than for the non-sequential approach.

Note that the convergence of our energy minimization algorithm depends on the used cost. The CNN based cost [33] which we use in our final experiments is more robust than the traditional cost (used in the previous experiment), and in this case the minimization process converges considerably faster than, for example, with the cost taken from [25]. In our experiments with the Middlebury training data set for most stereo images results do not improve after three iterations. In this section we perform our second experiment with the full training set. The sequential realization of our algorithm was applied with three iterations and non-sequential with six. The result of this experiment is given in Fig. 7. We again can see that the sequential approach yields faster convergence and the accuracy of the sequential approach is higher even in the case, when the number of iterations is half as many as in the non-sequential case.

VI-C OVOD versus WTA

This experiment demonstrates the advantages of the proposed OVOD approach. We found that this approach can improve the matching results even for the case when the cost processing is not implemented; the cost is considered as the BP marginals. This case is illustrated in Fig. 8. We can see that the OVOD solution is considerably more accurate than the solution obtained with the WTA operator, especially for metrics that count the occluded regions (all pixels). The drop of the error varies from 15% on Teddy to 65% on Vintage. Fig. 9 demonstrates the advantages of the proposed OVOD approach in the case when the BP process in our method pipeline is replaced by the cost filtering based on the geodesic distance filter (an approach similar to [38]). Again we can see that the proposed approach improves the accuracy of stereo matching.

VI-D Fully connected CRF models

In this experiment we compare our method with two similar approach. The first approach, the TSGO algorithm [25], uses the two step energy minimization technique and performs a left to right consistency check. To obtain a fair comparison we take the same cost as was used in [25]. Fig. 10 demonstrates that our approach leads to improved accuracy. The average rank on the Middlebury of the TSGO is 47 for this data set, while the ranking result of the proposed method with the TSGO cost would be 22.

The second approach we compare to is the cost filtering with geodesic distance affinity. This kind of filtering was used in [38]. To obtain a fair comparison we take the same CNN cost as was used in our main pipeline. Fig. 11 demonstrates that our approach leads to a significant drop of over 30% on average.

VII Conclusions

Our work proposes two main contributions. First is a new theoretical contribution in the field of MAP: we derive a fully connected model for the BP approach. We demonstrate that this model can be applied to the stereo matching problem and can improve the matching results in comparison with locally connected models. Our second contribution is the OVOD approach that allows to perform only one energy minimization process and avoids the cost calculation for the second view and the left-right check procedure. Application of our algorithm on the Middlebury data set reach state-of-the-art results and for several error metrics even first rank position. We think that this results can be improved by applying the fine tune part of the local expansion moves algorithm [33] as a post-processing step.

Acknowledgements

This work has been supported by the Spanish project TIN2015-65464-R, TIN2016-79717-R (MINECO/FEDER) and the COST Action IC1307 iV&L Net (European Network on Integrating Vision and Language), supported by COST (European Cooperation in Science and Technology). We acknowledge the CERCA Programme of Generalitat de Catalunya. We also acknowledge the generous GPU support from Nvidia.