Fusion of Range and Stereo Data for High-Resolution Scene-Modeling

This paper addresses the problem of range-stereo fusion, for the construction of high-resolution depth maps. In particular, we combine low-resolution depth data with high-resolution stereo data, in a maximum a posteriori (MAP) formulation. Unlike existing schemes that build on MRF optimizers, we infer the disparity map from a series of local energy minimization problems that are solved hierarchically, by growing sparse initial disparities obtained from the depth data. The accuracy of the method is not compromised, owing to three properties of the data-term in the energy function. Firstly, it incorporates a new correlation function that is capable of providing refined correlations and disparities, via subpixel correction. Secondly, the correlation scores rely on an adaptive cost aggregation step, based on the depth data. Thirdly, the stereo and depth likelihoods are adaptively fused, based on the scene texture and camera geometry. These properties lead to a more selective growing process which, unlike previous seed-growing methods, avoids the tendency to propagate incorrect disparities. The proposed method gives rise to an intrinsically efficient algorithm, which runs at 3FPS on 2.0MP images on a standard desktop computer. The strong performance of the new method is established both by quantitative comparisons with state-of-the-art methods, and by qualitative comparisons using real depth-stereo data-sets.


I. INTRODUCTION
Many computer vision methodologies, including dense 3D reconstruction [1], [2], gesture recognition [3], [4], and object detection [5] have benefited from recently-developed depth sensors. These sensors rely on active-light principles, including modulated-light and pulsed-light cameras, commonly denoted time-of-flight (TOF) [6], [7], or projected-pattern triangulation cameras [8]. Regardless of the working principle, however, these sensors provide low-resolution (LR) or mid-resolution depth maps that are inadequate for a number of applications such as 3DTV and film production. For example, many tasks in the film production industry can greatly benefit from a highresolution (HR) and high-quality depth map [9].
While HR depth maps can be obtained from multiple-view matching and reconstruction using standard color cameras, it is well known that stereo matching is problematic when the scene contains weakly textured areas, repetitive patterns, or occlusions; these situations are very common in both indoor and outdoor environments. Active-light sensors do not suffer from these limitations, although their own depth data are quite noisy in the presence of scattering, non-Lambertian materials, and slanted surfaces. The complementary nature of HR stereo and LR depth sensors leads to the design of mixed camera systems [10]- [18], which seem to be the most promising approach, at present, for high-quality 3D depth maps.
In this context, this paper addresses the problem of HR 3D reconstruction from the combination of a photometric camera pair and an active-light camera, provided that the multiplecamera setup is calibrated [19]. The combination of a stereo matching algorithm and of an active-light sensor raises the central question of devising a matching algorithm with the following features: (i) it considerably increases the resolution of the depth data, e.g., by a factor of ten, (ii) it eliminates depth-sensor errors wherever possible, (iii) it overcomes the limitations of stereo algorithms in textureless areas, and (iv) it is able to compete with a depth sensor in terms of speed. Hence, the availability of an efficient and robust stereo algorithm that takes advantage of LR depth sensors and that provides dense and accurate HR depth maps, possibly with subpixel accuracy, is particularly desirable.
To this end, we propose a 3D reconstruction method that merges depth-sensor measurements with photo-consistency stereo matching. We address the problem from the perspective of seed-growing, starting from a small number of control points whose disparities are then propagated to yield a dense disparity map. We show that this can be cast into maximum a posteriori (MAP) formulation (Sec. III), which leads, in turn, to a series of local optimization problems that are solved hierarchically by a novel region-growing process (Sec. V). While the proposed method may not reach the global optimum, it allows us to devise an intrinsically efficient methodology that bridges the gap between global optimizers based on Markov random fields (MRF) and locally-optimal winner-takeall (WTA) strategies (Sec V-E).
Efficient stereo-only or stereo-depth fusion methods rely on control points, by exploiting either feature correspondences in stereo [20], [21] or depth data in fusion [14], [15], and they assume that these points are of very good quality. A key contribution of this paper is that this requirement is relaxed in order to devise a method tolerant to bad control points. We propose to truly combine LR depth-sensor data with HR rich photometric information, whenever and wherever possible, showing that fusion is helpful, even in the early stage of depth initialization (Sec. IV-A & IV-B). The data term of the proposed MAP formulation benefits from a new cross- Fig. 1: The pipeline of the proposed depth-stereo fusion method. The low-resolution (LR) depth data are projected onto the color data and refined to yield a high-resolution (HR) sparse disparity map. Starting from these disparity seeds, an upsampling process provides an initial HR dense disparity map. Both the HR seeds and the initial dense disparity map are then used by the region-growing depth-stereo fusion to produce the final HR depth map. A prominent feature of our method is that fusion takes place at several data processing stages. correlation function, which provides real-valued disparities via subpixel corrections computed in closed-form (Sec.V-B1), and takes advantage of a depth-guided cross-product aggregation (Sec. V-B2). Moreover, the data-term adaptively merges the stereo-and depth-consistency terms guided by the scene texture and the camera geometry (Sec. V-C). These advantages lead to a more selective growing-of-disparities process, thus preventing the algorithm from propagating erroneous depthsensor data -a phenomenon that is often associated with propagation techniques -as experimentally verified (Sec. VI). It is important to note that the proposed method can be used 'as is' with any depth sensor, as it requires neither sensordependent confidence maps nor a sensor-dependent model. The proposed fusion pipeline is illustrated in Fig. 1.
Supplementary materials, in particular image datasets and Matlab code are available online. 1

II. RELATED WORK
We review pure stereo matching methods with emphasis on local algorithms, owing to their computational suitability for HR images. We also review upsampling methods and depth-stereo fusion methods. More detailed surveys of stereo matching algorithms and depth-stereo fusion methods can be found in [22] and [9], respectively.

A. Pure Stereo Matching
Stereo matching methods can be broadly classified into global and local [22]. Global algorithms [23] typically adopt an MRF formulation and solve a single optimization problem based on a MAP criterion. Despite their superiority over local methods, global algorithms are extremely time-consuming, and hence unattractive for fast fusing of depth-sensor data with high-resolution images. Local algorithms solve per-pixel optimization problems and the state-of-the-art methods build on adaptive cost aggregation [24]- [27]. Most methods, however, must visit the entire cost volume to find an optimal disparity value at each pixel. This volume grows rapidly with respect to the input, as the width of the image typically multiplies the number of disparities. Therefore, although they are able 1 https://team.inria.fr/perception/research/dsfusion/ to provide LR disparity maps in real-time, they remain slow and subject to memory issues in HR stereo. Note that global algorithms need several approximations to obtain LR disparity maps in real-time [28].
More interestingly, algorithms that rely on control-point correspondences [29] are drastically more efficient since they avoid visiting the whole cost volume. Region-growing approaches start from reliable but sparse correspondences (seeds) and propagate them in textured areas [20]. [30] suggests a similar propagation scheme where orientation-consistent disparities are propagated to neighbors at the cost of finding a plane equation per pixel. [21] proposed a generative model, where the prior disparity comes from a 2D triangulated mesh whose vertices (control points) are obtained from matches between low-level features. As with [20], textureless areas remain intractable and the final map is reliable only when the matches are dense and uniformly distributed over the images. Notice that the proposed method relies on the idea of control points that are transferred from a depth sensor, thus avoiding the limitation owing to untextured areas.
To obtain continuous disparities that are required in many scenarios, e.g., 3D reconstruction [31], local stereo algorithms typically employ two strategies: (i) fitting a curve around the correlation peak [14], [22] or (ii) integrating an intensity interpolation kernel into the (dis)similarity function [32] whose optimization leads to subpixel correction. The latter is also the case in the fusion framework of [16], [17] that inherently takes advantage of inter-pixel depth estimations.

B. Depth Upsampling and Depth-Stereo Fusion
Any prior depth information, even at low resolution, is likely to help dense disparity estimation. Apart from a naive interpolation, the bilateral filter [33] can post-process an interpolated map using color HR images for guidance [31]. Alternatively, a joint bilateral filter applies spatial and range kernels to the LR depth map and HR color image respectively, so that upsampling is a by-product of filtering [34]. Upsampling methods, however, are limited to clean depth LR data and cannot reconstruct accurate HR maps when LR data are delivered by depth sensors.
The above limitation of the upsampling methods gives rise to fusion approaches that can merge depth and stereo data in either early or late stages, once the sensors are calibrated. Late fusion suggests merging two depth maps, one obtained with stereo and one from the range sensor, possibly upsampled [18], [35]. The majority of fusion methods, however, merge the stereo and range-sensor data at an earlier processing level. [10] estimates a LR TOF-based disparity map, which initializes a coarse-to-fine stereo scheme [36]. A semi-global dynamic programming approach is followed in [14] with the TOF-based disparities, wherever available, being considered as error-free matches. In a global framework, [11] produces mid-resolution depth maps by merging TOF and stereo data within a graphcut stereo-matching paradigm; each energy term exploits both modalities. Likewise, MRF-based formulations have been proposed [12], [13], [37]. In [37], ground control points reflect an extra regularization term in the MRF energy function. The work described in [12] uses an MRF scheme to merge depth distributions of each sensor alone, but the goal is a LR depth map. [13] extends [12] by means of weighted fusion. A balanced fusion (based on several confidence maps) within a total variation model is also proposed in [17]. A similar variational model that infers the HR map in a coarse-to-fine manner is adopted by [16]. Note that, unlike most fusion methods and similar to the proposed one, [16] is not tuned to a specific depth sensor.
More closely to the present work, [15] fuses the data within the seed-growing method of [20]. In particular, TOFbased disparities constitute seeds while a triangulation-based interpolated TOF map regularizes the seed-growing process. When the TOF data are noisy, however, this approach tends to produce incorrect disparities, and to propagate false positives during the growing process. The proposed method differs considerably from [15] in terms of depth initialization, cost function, and fusion strategy. The proposed initial map is robust to depth discontinuities while it also guides the cost aggregation inside a window. Moreover, our likelihood term integrates functions that are capable of providing sub-pixel corrections [32]. This turns out to be very beneficial, not only for the continuous nature of the final map, but also for the growing process itself, thus propagating high-confidence disparities only. Note that the subpixel disparity correction is obtained from a closed-form solution -an interesting feature for efficiency. Finally, our algorithm benefits from an adaptive fusion scheme that better balances the contribution of each modality (depth or color), and that results in fewer unmatched pixels.

III. PROBLEM FORMULATION
The main mathematical notations that are used throughout the paper are summarized in Table I. As discussed, the direct upsampling of LR depth data suffers from limitations, in particular when the upsampling factor is high. 2 Therefore, our goal is to build D by jointly taking advantage of both sensing modalities. Given a proper calibration, (e.g., [19], [38]), the  with D, D 0 and D↓ respectively, with dp ∈ D, d 0 p ∈ D 0 and d p↓ ∈D↓ tp: Subpixel disparity correction of p µp = (dp; tp): Disparity-correction pair referred here to as metadisparity with |tp| < 1 M = {µp}: Set of meta-disparities S = {sp}: Set of observed stereo pixel intensities Ip: Intensity of pixel p u(x, y): a vectorized (zero-mean) form of an intensity window centered at the 2D position (x, y) Np: Neighborhood of pixel p med: 2-d median operator I R , D R : Intensity and disparity maps defined on a sub-region R E S (dp), E D (dp): Stereo-based and depth-based energy of p for given disparity g(ξ; γ) = e −|ξ|/γ : Exponential mapping of ξ with scale γ mapping of the LR depth image onto the rectified HR color images will typically yield a sparse disparity map D↓, which is almost evenly distributed across the HR grid. Note that in [20], [21], the initial matches between the stereo images do not correspond to a uniform sparse version of D, as they are unpredictably distributed, due to the reliability of properly detecting interest points in images. Instead, the sparse map obtained with a depth sensor can be used to guide a stereo algorithm, regardless of the presence or absence of scene texture.
We propose to model the estimation of D, and therefore the map D, as a maximum a posteriori (MAP) problem, based on the available depth and stereo observations. However, instead of immediately using D↓, we first estimate a dense initial map D 0 and its associated set D 0 . Then, we obtain the final disparity map by solving the following optimization problem: where P (D|S, D 0 ) is the posterior distribution of disparities given the observations S and D 0 . The proposed depth initialization method is described in Sec. IV and the proposed solution to the MAP formulation (1) is described in Sec. V. Note that a reliable estimation of D 0 is quite important, since it guides several components of the fusion methodology.
IV. DEPTH INITIALIZATION A two-step approach is proposed in order to obtain the initial disparity map D 0 . First, we refine D↓ to deal with mapping errors. Second, we upsample the refined sparse map in a novel way using color information to obtain the initial dense map D 0 . As shown below, this leads to an initialization robust to depth discontinuities, which in turn helps the growing.

A. Sparse-Depth Refinement
We assume a camera setup with a depth camera mounted in between the two color cameras; other depth-stereo setups are (a) (b) (c) Fig. 2: (Best viewed on screen) (a) The mapping of depth data onto the left image causes artifacts in the presence of depth discontinuities. A cascade of (b) geometry-consistency and (c) color-consistency filters refines the sparse disparity map. Depth values are color-coded from red (close) to blue (far). Fig. 3: (Best viewed on screen) Depth upsampling results using (a) triangulation-based interpolation [15] after cutting big triangles, (b) joint bilateral filter [34] and (c) our method. The depth values are color-coded from red (close) to blue (far), while black areas correspond to non-available values. The white edges in close-ups show the color edges of the image. discussed in [9]. Regardless of the depth sensor technology and type, the projection of the depth map onto the left and right images implies a parallax effect, and hence occlusions. Moreover, this causes gaps as well as areas with overlapping depth data close to depth boundaries [15]. In the case of TOF cameras, these areas are further contaminated from jumpedge errors [39], or 'flying pixels' [40], while a structuredlight camera, e.g., Kinect, leaves more gaps due to the offset between the position of the light projector and the position of the infrared sensor. Fig. 2(a) illustrates the artifacts that we briefly discussed: flying pixels and depth-data overlap in the top and bottom closeups, respectively. In order to eliminate these artifacts we apply a geometry-consistency cascade of two filters: the first one removes isolated pixels (mostly flying pixels) and the second one keeps the foremost pixel inside a window to compensate for the above mentioned overlap. An example of applying this filtering is shown in Fig. 2(b). In practice this does not fully refine the sparse depth map. We still observe mismatches near depth discontinuities, because of depth bias and calibration errors. Therefore, a second filter that imposes color consistency is applied, as described below.
We consider a window centered at p↓ and split into four equally sized sub-windows W i , i = {1 . . . 4}, such that their intersection is only the pixel p↓ (see Fig. 4). The output, d p↓ , of the filter is: with The output d p↓ is the median disparity of the adjacent subwindow whose median intensity is closest to that at p↓. For color images, the term |I p↓ −med(I Wi )| can be replaced by the average deviation from the median, over the color channels. This filter leads to a further refinement near depth discontinuities which are pathological areas for stereo algorithms. The result of this kind of filtering is shown in Fig. 2(c). Note that both refinement filters apply to sparse locations only so that their complexity is negligible. Fig. 4: The window split for the color-consistency filter. The pixel p is linked with the closest (shaded) sub-window in terms of the color consistency (links represent color distances from I p to med(I Wi )).

B. Sparse-Depth Upsampling
While a naive upsampling of the sparse disparity map could be performed, e.g., [15], thereby producing an initial map, strong depth discontinuities are likely to contaminate such an interpolation. Alternatively, cross-bilateral filtering [41] or joint bilateral upsampling [34] may upsample the map using the HR color image as a guide. The filter support of these methods jointly counts on spatial and range (color) kernels. However, both methods apply a linear smoothing once the filter support per pixel is computed. Instead, we propose a non-linear upsampling strategy that builds on the median filter.
Specifically, the depth (or disparity) at a dense pixel location p is initialized by where N c p is a constrained neighborhood of p, that is N c p ⊂ N p , which contains only sparse depth measurements whose color is consistent with I p : with g(ξ; γ) being an exponential mapping (see Table I).
Unlike common bilateral filters, our upsampling process makes a more definitive selection of pixels, thus preserving the depth discontinuities of the scene, while filtering out some of the noise in the depth data. Once the N c p is defined, one can optionally consider a spatial kernel and compute a weighted average instead, in order to better deal with slanted surfaces. Fig. 3 compares our initialization with the upsampling results obtained by [15] and [34]; in the detailed views, the intensity edges are also shown. The proposed method provides more discriminative depth boundaries that coincide with intensity edges. Missing values may be observed in highly textured areas, since N c p may be an empty set. In this example, the radius of N p is 20, γ c = 10, and c = 0.2. The same radius is used for the method of [34] while the scales for spatial and color kernel are 10 and 20 pixels, respectively. The geometryconsistency filter reasonably applies in all cases while our method benefits from our color-consistency refinement as well.
Since the median operator is chosen to account for outliers within a window, the mean operator can be invoked instead when the depth variance almost vanishes (homogeneous areas), thus drastically reducing the computational burden of the upsampling process. Since the vast majority of pixels belong to such areas, the complexity of our filter approaches that of joint upsampling filter [34]. Note that the latter has been extended in [18] by integrating color segmentation results.

V. DEPTH-STEREO FUSION
Let d p ∈ D have N possible discrete states; the goal is to estimate the disparity (state) of each HR image location through the MAP formulation (1). Once D 0 has been initialized, one can assume that S and D 0 are conditionally independent, so that the posterior distribution of (1) can be decomposed as As mentioned, global solutions are prohibitively expensive for high-resolution disparity maps. Therefore, we focus on approximate solutions that allow for the decomposition of the global optimization problem into many local (per-pixel) optimization problems. The proposed method is based on the seeded region-growing framework [20], [42], where the known message of a location (the parent) is propagated to its neighbor (the child). This implies that the estimation of d p is also conditioned by a parent known disparity d pa p , hence dealing with the principal where d pa p is the parent of d p and the probability has a uniform distribution. In other words, the disparity range for each node is a function of the assigned disparity of his parent. We consider a narrow support area, i.e., a constant low value for r, e.g., 1 or 2; node-dependent parametrization of r is left to future work.

A. MAP as Energy Minimization
As is customary, likelihoods are chosen from the exponential family which leads to an energy minimization framework. Our model assumes an energy-dependent distribution (Boltzmann) for P (s p |d p ) and a Laplacian one for P (d 0 p |d p ): Based on (7-10), Fig. 6 shows an example with the distribution of d p being constrained by single or joint observations. Because of the exponential terms, the pixel-wise maximization of the posterior distribution reduces to the minimization of the local energy where is the (stereo) data-term which is defined below. The term E D (d p ) guides the inference in textureless areas, while it penalizes mismatches due to depth discontinuities when the latter are well preserved in D 0 . Notice that a smoothness constraint is implicitly enforced because of the prior term P (d p |d pa p ) owing to the low value of r. In other words, the support area of E(d p ) is truncated as shown in Fig. 6, and the pixel-wise optimization problem becomes The order of visiting pixels and solving (12) obeys a mostconfident first-solved rule, as discussed in Sec V-D. Below, once the data-term is defined, we modify (11) to adaptively combine the data-term with the regularizer.

B. The Data Term
The proposed data-term benefits from two properties: its ability to provide energy for subpixel disparities and to adaptively consider pixel-wise similarities within a window. Apart from the advantage of a locally continuous map [31], it is more important to note that these properties lead to a more selective growing in areas with texture (owing to the former) and varying depth (owing to the latter) because of more accurate and reliable local energies.
1) Subpixel disparity: To obtain energy at interpixel locations, we make use of meta-disparities (d p ; t p ) (see Table I). The estimation of t p builds on [32], which has been shown to be superior to the parabola-fitting around the peak correlation: the pixel p = (x, y) of the left image is matched with the subpixel p = (x + d p + t p , y) of the right image by defining a bounded correlation function C dp (t p ) (|C dp (t p )| < 1) which is maximized with respect to t p for a given d p . This allows us to define our data-term as where Consequently, the total energy (11) is parameterized by t * p as well. Note that if t * p can be estimated via an analytic solution, the extra computational cost of subpixel disparity estimation is negligible. It becomes clear that, since C dp (t * p ) ≥ C dp (0) in (13), more accurate energy values are estimated and more reliable messages are propagated.
Since invariance to photometric distortions is important in stereo matching, we propose to adopt the normalized correlation coefficient, and one of its variants [43], for C dp (t p ). While the former has been already extended to deal with subpixel disparities [32], the latter has never been extended before. Both use zero-mean vectorized forms of the windows around p and p , let u L (x, y) and u R (x + d + t, y), with the latter being written via a first-order Taylor approximation as where ∆ is a difference operator along the x−axis. This is also the case in variational framework [16], [17] where subpixel accuracy is obtained by the early interpolation of the intensity, rather than the late interpolation of correlation around the peak [22].
Enhanced Correlation Coefficient (ECC): The ECC function [32] results from the integration of the above linear interpolation kernel into Pearson coefficient: If the denominator of (15) is non-degenerate, then C P d (t) is a quasi-concave function of t and its maximization results in a closed-form solution [32].
Enhanced Moravec Correlation Coefficient (EMCC): The Moravec coefficient [43] replaces the denominator of (15) with the mean of the variances. This also allows us to introduce a left-right symmetry, thereby estimating a left disparity −t/2 and a right disparity t/2 instead of t (such a modification with ECC leads to a complex optimization problem). The enhanced Moravec correlation coefficient (EMCC) is defined by Note that one can easily show that the integration of an interpolation kernel into the cost function of [15] is equivalent with the EMCC scheme. Although (16) is a rational function of t, the next proposition guarantees that the maximizer has an analytic form. We refer the reader to the appendix for the proof and the exact maximizer.
Proposition I: A rational function of two second-degree polynomials, as in (16), attains at most one global maximum if the denominator is non-degenerate; its maximizer is given by a closed-form solution.
Note that, the estimation of t could be unreliable for strictly homogeneous areas. The value of the window variance or entropy is a good criterion to assess the reliability of subpixel correction, and to enable it accordingly.
2) Adaptive similarity aggregation: The best-performing local stereo algorithms benefit from an adaptive cost aggregation strategy [24]. This strategy is based on the assumption that depth discontinuities are most likely to coincide with color discontinuities, so that each pixel within a window contributes differently to the (dis)similarity cost based on its spatial and color distance from the central pixel. However, only a few color edges correspond to depth edges and the above assumption should be followed only in the absence of any prior information about the depth. Since in our scenario the prior depth information is available, the spatial and color consistency can be replaced by a depth consistency term.
To be specific, we adopt here the exponential g(·) (see Table I) to compute pixel-wise weights w q : with q ∈ N p . The weights apply element-wise to u L , u R , ∆u R and ∆u L in (15) and (16). In other words, we compute the subpixel correction and the optimum local energy after down-weighting pixels in the window that belong to another surface compared to the one of the central pixel. It becomes clear now that not only the term E D (d p ) but also the stereo term E S (d p ) in (11) benefit from an upsampling method that is robust to depth discontinuities. Note that, even if the initial depth map is biased, it is sufficient enough to guide the aggregation step within the window.

C. Adaptive Fusion
While a constant fusion may be reasonable for specific types of scenes (e.g. highly-textured scenes), an adaptive balance of the terms in (11) is usually preferred. This implies that the less we count on the data term E S (d p ), the more we should count on the regularization term E D (d p ) during the inference. 4 This suggests a convex combination when the scene point of p is viewed by all cameras.
A summary of methods that perform weighted fusion is discussed in [17]. However, most of them consider TOF-based weights for the regularizer (e.g., [13]) which contradicts our goal of a sensor independent fusion. Moreover, directly using the confidence map of a TOF image is not a good strategy [40]. Therefore, we only rely on stereo data to obtain the mixing coefficients. It is well known in stereo or optical flow that the matching of a point is reliable when its associated image patch contains sufficient texture [44], [45]. Since a good indicator for the texture presence is the image entropy [44], an entropy filter provides us with am adequate reliability factor e p for each window centered at p.
Let us now consider the left image as reference and compute the initial left-to-right disparities. Likewise, we can build a right-to-left disparity map based on the right image, and a cross-checking of these maps can provide an estimation of the major occlusions due to strong depth discontinuities, with respect to the reference image. We refer to these areas as stereo-occlusions and we denote them as Ω SO . Recall now that some points in the left image are not seen by the depth camera, and that this gives rise to gaps in the initial disparity map which can be easily detected. We refer to these areas as depth-occlusions and we denote them as Ω DO . It becomes obvious that the evaluation of E S (d p ) and E D (d p ), in Ω SO and Ω DO respectively, should be avoided. Hence, we propose the following adaptive fusion with the pair (η S p , η D p ) being defined as where Ω defines the whole image area and e p is the normalized output of the entropy filter. We intentionally add the last 4 Here we omit the disparity correction t. case in (19) which shows that the fusion in Ω SO Ω DO is meaningless and a post-filling method should be followed.
It is important to note that the so-called ordering constraint in stereo is valid when large foreground objects appear in the scene, while it is violated when very thin objects are close to the camera (see Fig. 7). The former implies Ω DO Ω SO = Ω DO and the latter implies Ω SO Ω DO = ∅. While the area Ω SO \ Ω DO can always be predetected, the area Ω DO \ Ω SO is safely predetected only when the ordering constraint is not valid. 5 This is because Ω SO is detected from the crosschecking of disparity maps that already suffer from depthocclusions, since they are computed from the depth-to-stereo mapping. Ideally, if the complexity is not an issue, a stereoocclusion detection scheme (e.g. [46]) based on stereo image pair could be enabled beforehand. Therefore, we prefer to not grow the area of depth-occlusions. Optionally, these areas can be filled in a post-processing step. Fig. 7(c) shows the areas for the example of Fig 2. Points that have been removed during the refinement step (Sec. IV-A) are also marked as depthocclusions, while outliers in the initial maps can produce false positives for Ω SO , e.g. the red points on the right arm of the person. Our region-growing method, however, is able to compensate for such errors.

D. The region-growing algorithm
As has been explained, our method solves pixel-based optimization problems in a region-growing manner, based on the seeds contained in D↓. Since the initial disparities may be noisy and biased, they do not reflect true matches as opposed to [15]. This means that we exploit d p↓ to restrict the disparity range of its neighbor, but once a disparity value is assigned to the latter, p↓ is reset to a pixel p with unknown disparity. For our convenience, the set D↓ is augmented to an initial set M of meta-disparities (d p ; 0) with the same cardinality. We also denote with N (µ p ) = {µ j p } 4 j=1 the image-based neighbors of µ p , that is, µ j p = (dp; tp) withp being any of the four immediate neighbors of p. Note that µ j p does not necessarily belong to M.
Algorithm 1 describes the growing process. The algorithm starts by sorting the elements of M based on their energy 5 It would be possible to detect all areas if one knows a priori that the constraint is everywhere valid or invalid. Set f v (µ p ) = true 8: for all µ j p ∈ N (µ p ) with f a (µ j p ) = false do 9: µ j * p = arg min E(µ j p ) (Eq. 20) 10: if E(µ j * p ) < T then 11: Set f a (µ j * p ) = true 12: Push µ j * p in M w.r.t. sorting 13: end if 14: end for 15: until card(M) is fixed 16: Compute the dense disparity map D from M 17: return D.
value, while it initializes to false the visit-and assignmentflag of all candidate meta-disparities. Next, it considers the lowest-energy µ p with false visit-flag and switches this flag to true. Then, it assigns values at each µ j p ∈ N (µ p ) with false assignment-flag based on the following minimization scheme where E(dp; tp) = η Ŝ p E S (dp; tp) + η D p E D (dp) (21) and d * p is the optimum disparity of the parent node ofp. Note that (21) extends (11) by adding the subpixel disparity parameter into the stereo term, and making the fusion pixeldependent. In other words, what we do for each neighborp with false assignment-flag is the following. For each candidate integer value dp ∈ [d * p − r, d * p + r] the optimum t * p that minimizes the term E S (dp; tp) is obtained based on (16) or (15) and the local energy E(d p ; t * p ) is computed from (21). Among the 2r + 1 values, the disparity minimizer d * p is finally chosen and the corresponding subpixel correction is assigned, as it is shown in (20). Recall that r has a low value in contrast to conventional stereo algorithms where it equals the whole disparity range. If E(d * p ; t * p ) < T , then the meta-disparity (d * p ; t * p ) is pushed into M with respect to the sorting, and its assignment-flag becomes true. The above process is repeated with unvisited meta-disparities until the cardinality of M remains fixed. The visiting order depends on local energies, since the lowest-energy meta-disparity is always picked from the stack. The final set M corresponds to the final dense disparity map while a post-filling method can deal with missing disparities; their number heavily depends on threshold T . It is important to note that the algorithm cannot get stuck in a loop, because it propagates disparities in a tree structure, and a true visit-flag can never be reset to false.

E. Comparison with Other Inference Formulations
State-of-the-art stereo or fusion methods adopt an MRF model that provides a straightforward way to integrate multiple sensor data. If we recall equation (6), the term P (D) can be written as a Gibbs distribution whose energy is a sum of potential functions over maximal cliques [47], hence a sum of pairwise potentials when a 4-pixel neighborhood is considered. This in turn offers a smoothness term in the global energy equation that leads to piecewise smooth disparity map. There exist exact solutions for such models, under very specific conditions [23]. Their computational complexity, however, becomes more and more severe as the number of states increases, thus becoming prohibitively expensive in the case of HRstereo. Therefore, one has to focus on feasible approximate inferences that decompose the global optimization problem into a series of local optimization problems.
MRF stereo generalizes the winner-take-all (WTA) approach [22] which can be seen as the simplest inference. Fig. 8(a) and (b) show an 1D-grid example of the MRF and WTA networks. In essence, links between disparity nodes in WTA disappear and each node is connected only to input data (observations). This implies a uniform prior distribution P (D) and the maximization of P (D|S, D 0 ) reduces to a set of independent pixel-wise maximization problems of type P (d p |s p , d 0 p ). As with the proposed method, we tacitly assume, that s p are not intensities but they represent the stereo data of local windows centered at p and its candidate correspondence. This is necessary in WTA since the likelihood cannot count on single pixels only, while it can optionally used in global MRF solutions as well. It is worth noticing that our algorithm switches to a WTA solution when r = N , since the visiting order of nodes and their connectivity becomes meaningless and the uniform distribution P (d pa p |d p ) is not truncated anymore. Another MRF approximation that finds pixel-wise optimizers without breaking the connectivity between nodes is the iterated conditional modes (ICM) method [48]. After a proper initialization of all disparity nodes, ICM-like schemes visit one node (pixel) at a time and assign the disparity value that maximally contributes to the global posterior distribution. On the contrary, our approach requires the initialization of few nodes only, at least one in principle. Moreover, the graph in our case is a set of directed trees, i.e., a forest, so that each non-root node has only one parent. It is important to note, however, that it is the output of our algorithm that defines the final network since nodes may be re-ordered, links may be cancelled and arrows may be reversed during the inference.
To be more specific, let us consider the graph of Fig. 8(c) and let us assume the clamping of d 2 and d 5 to the observation d 0 2 and d 0 5 . Starting from d 2 , we can 'propagate' d 0 2 to its neighbor d 3 (and d 1 ), i.e., to look for the optimum value of d 3 but being strictly conditioned by the value d 0 2 . Note that ICM would look for the best assignment of d 3 value by taking into account both initial values d 0 2 and d 0 4 . Moreover, ICM would search, in principle, among all (here N ) states for the optimum d 3 's assignment, while our scheme looks only around d 0 2 , namely d 3 takes values in {d 0 2 − r, ..., d 0 2 + r}, with r being a small integer. Once d 3 's assignment is done, d 0 2 can be passed to d 1 in a similar way. Next, another node is visited, here one of d 3 and d 5 , and its disparity is propagated to its neighbors. As a result, the principal graph of Fig. 5 is iteratively considered. Note that ICM assigns a new value to d 3 anyway, while our scheme invokes a criterion that validates the assignment. If d 3 's assignment is not valid, our algorithm will possibly assign a new disparity value to d 3 only after d 4 's assignment. Fig 8(d) and (e) show possible final graphs obtained by different realizations of our algorithm. In the example of Fig 8(e), the initial disparity of d 5 was cancelled by the validation process, thus all nodes were filled due to d 2 .
It is now clear that, as opposed to WTA solution, the final inference we obtain depends on the visiting order of nodes. The visiting order of the ICM scheme is either fixed in advance (e.g. raster-scanning), or random. Inspired by [47] and [20], however, we instead adopt a highest confidence first (HCF) scheme that suggests visiting the nodes based on their local evidence (energy). This means that we keep the nodes sorted with respect to their energy, and we visit each time the leastenergy node that has not been visited yet. For instance, in Fig. 8(d), starting with d 2 implies that d 2 is more confident than d 5 . The above validation process relies on thresholding the local energy as explained in Sec. V-D. Table II summarizes some properties of MRF-based solutions (Graph-cut [49] and ICM [48]), WTA, and our approach.

VI. EXPERIMENTS
In this section, we evaluate our algorithm, and quantitatively compare it to the state-of-the-art, based on both simulated and real data-sets. We also test our algorithm and provide qualitative comparison on a new and challenging dataset.

A. Simulated Data
We use the Middlebury database, and focus on a challenging data-set which contains 1.5MP images (1300 × 1100) along with ground-truth maps (GTM) [50]. To simulate an LR disparity map from another viewpoint, we proceed as follows. Given the calibration parameters, we transform the groundtruth disparities into 3D points, as viewed from the midpoint of the baseline. We then apply a 3D rotation to the pointset and we downsample the points by a factor of 10. The rotation is such that the average disparity bias is more than 2 pixels. Finally, we back-translate the points into sparse biased disparities and we add colored noise, that is, a 2D midfrequency sinusoidal signal with peak-to-peak distance equal to 2σ, where σ denotes the noise deviation. Note that [15], [31], [34] only downsample the GTMs in their experimental setup. Algorithm performance is quantified in terms of the so-called bad matching pixels (BMP) percentage in the nonoccluded areas [22], i.e., (1/N o ) p (|D p − G p | > δ) where G is the GTM and N o is the number of the non-occluded pixels. While the threshold level δ = 1 is mostly used for mid-resolution images, HR stereo justifies the value δ = 2 as well [21]; a value δ < 1 is chosen when subpixel accuracy is to be evaluated.
The same parameter settings are used for our method, in all of the experiments. The radius of the upsampling filter is 20 and the values γ c and e c are 10 and 0.2 respectively. Because of the propagation strategy, we choose a relatively small window, i.e., 9 × 9. The local energy in (11) and the weights in (17) are obtained with λ = 0.01 and γ d = 5, respectively. We enable subpixel correction when the normalized entropy in the (left) window is above 0.4. As mentioned, we use a fixed (and strict) search range around the disparity parent, that is r = 1, which leads to the most efficient solution. As for the threshold, we set T = 0.5. Note that the energy validation threshold implies a tradeoff between accuracy and density. We recommend setting a middle threshold value and post-filling sparse missing disparities, e.g., with the upsampling filter, rather than using a high threshold that incorporates erroneous disparities in a fully dense map. The density obtained with this strategy is about 90% in HR images. As with all algorithms, large remaining gaps are filled with a streak-based filling [21]. We refer to our methods as Fusion-ECC (F-ECC) and Fusion-EMCC (F-EMCC).
Before comparing with the state-of-the-art, we show the performance gain in terms of the new modules that are integrated compared to EPC method [15], thus quantifying the contribution of the new data-term and the adaptive fusion. Note that [15] (EPC) follows a seed-growing approach by using a quadratic model for both stereo and depth consistency terms respectively. To better evaluate the contribution of the enhanced correlation coefficients presented in the data-terms, we use LR images (450 × 375 on an average) whose GTM contains subpixel disparities. The noise deviation in the sparse map is σ = 2. We intentionally do not fill any remaining large gaps, in order to assess the net contribution of each module, and we compute the error for the filled area only (85% density). Table III shows the BMP error averaged over eight images. We also evaluate our approach when none of the modules are enabled, i.e. pure correlation is used along with a fixed fusion of the terms; we refer to this method as simple fusion. All of the variants start from the same initial depth map, obtained by our upsampling method. As can be seen, both the data-term and the adaptive fusion process contribute to a better reconstruction, compared to simple fusion. The use of the proposed energy data-term leads to more accurate  [23], [49] results, while the adaptive fusion eliminates large errors (see the error with δ = 2). Even the simple fusion has a lower BMP than EPC, owing to the different stereo-and depth-consistency terms (we use a linear model for the latter). The two proposed criteria behave similarly, with the F-ECC being slightly more accurate, since it achieves higher correlation values (see also Table IV). Note that the results are systematically worse if we use the initialization of [15].
To compare with the state-of-the-art, we implemented the upsampling methods of [34] and [31] (two-view version), as well as the MRF-based fusion [12] by using the MRFstereo toolbox of [23], referred here to as F-MRF. While [12] uses belief propagation, we experimentally found that Graph-Cuts [49] perform better. Specifically, we tried ten different parameter settings and we found that the best performing algorithm is the expansion mode with Birchfield-Tomasi cost [51] truncated at 7, linear disparity differences truncated at 5 and quadratic cost for the smoothness-energy; the weights for the stereo-, depth-and smoothness terms were 1, 1.2 and 10. We refer here to [12] instead of [13] for MRF-based fusion since the latter relies on a TOF-based reliability fusion which cannot be implemented here. Fig. 9 plots the BMP curves of the upsampling and fusion algorithms as a function of the noise deviation for the challenging low-texture HR image Lampshade1. We just add noise here in the down-sampled GTMs. Except for [34], all schemes start from the same HR map, obtained by a naive interpolation, while its BMP curve is plotted as well. As can be seen, F-MRF and EPC are more affected by the noisy prior disparity, in contrast to the proposed algorithm, which is less sensitive to initialization. It is clear that the pure up-sampling methods provide acceptable results only when the initial LR disparity map is very accurate.
We now proceed with a detailed comparison including efficient and well-known stereo algorithms as well. Specifically, we include four recently proposed methods, [21] (ELAS), [20] (CGS), [25] (FastAgg), [27] (NonLocalAgg) and two MRFbased stereo algorithms, Graph Cuts (GC) [52] and constantspace belief propagation (CSBP) [28]. The top-performing lo- cal algorithms, FastAgg and NonLocalAgg, build cost volumes that depend on both the image size and the disparity range. This leads to a huge memory footprint (∼ 3GB) in the case of HR images, and the authors' implementations could not be run as is. In order to be able to run FastAgg, the cost volume has been split into slices and cached on disk. The NonLocalAgg was run using the maximum allowed resolution while the disparity map produced by the algorithm was finally upsampled. Note that both methods invoke a left-right consistency checking, combine color and gradient information and enable refinement steps. Authors' implementations for ELAS, GCS, GraphCut, CSBP (local-minima version+bilateral postprocessing) were used in the comparisons, with the default settings suggested by the authors. We also implemented the ICM algorithm for the MRF-based fusion using the same parameters with GC. For a fair comparison, all fusion schemes merge stereo data with the same initial disparity map, which is obtained here by our upsampling process. Due to the simulated experimental setup, however, the error of the initial map obtained from this process is very close to that of Kopf et al.'s method [34] (the average difference of their BMP error is below 0.5) and is thus omitted. Table IV provides the BMP error for eight HR images with error threshold δ = 1, while the corresponding table for δ = 2 is given in the appendix. Bold and underlined numbers mark the lowest and second lowest errors per column. Weakly textured scenes (Lampshade1, Monopoly) seem to be problematic for conventional stereo algorithms, while cylindrical surfaces (Bowling2, Baby1) present another challenge. It is not surprising, however, that stereo methods outperform fusion methods in highly textured images (e.g. Rocks2), or in images with many thin objects (e.g. Art), since the sparse noisy initial map negatively affects the fusion. ELAS and FastAgg behave better than other stereo algorithms. Similar results would be expected from the NonLocalAgg method, if we were able to run it at full resolution.
Stereoscopic and depth data are better fused in general by the proposed criteria than the other fusion methods. EPC  and up-sampling methods verify the sensitiveness to their initialization, with the former being more effective due to the fusion process. As far as the MRF solution is concerned, the benefit due to the depth data is verified from the results, i.e. the F-MRF (GC) scheme behaves better than pure-stereo GC. Moreover, GC in fusion provides better results than the simple ICM algorithm. Unlike the proposed criteria, F-MRF deals better with thin objects, since it does not aggregate costs in a window, hence the lower error in Art and Reindeer. Recall that [12] refines the TOF depth map based on stereo data, without increasing the resolution. By putting aside the high complexity, it seems that MRF-based solutions need to be reformulated for HR stereo-depth fusion, e.g. high-order connectivity might be more helpful, semi-global solutions could be investigated, and conditional random fields [50] might need to be extended to the fusion framework.
Between the two proposed criteria in our fusion scheme, ECC and EMCC, it is the image content that makes one outperform the other. While F-ECC may be slightly better on average, F-EMCC deals better with images of very low texture, e.g. Monopoly, which supports Moravec's argument for introducing MCC [43]. Moreover, F-EMCC is more affected by the filling, as it provides less dense maps than F-ECC, provided that the threshold is the same. Fig. 10 shows the average performance of all fusion competitors as a function of the error threshold δ, while the stereo baseline of FastAgg is added for reference. In essence, this figure reflects the distribution of errors. Evidently, the contribution of the LR depth prior in the fusion schemes is verified, as opposed to the stereo baseline whose performance is bounded. The proposed schemes are more accurate compared to the fusion baselines. However, F-MRF provides lower errors when the tolerance is not that strict (δ > 2) owing to its global smoothness constraint. Note that the performance of [34] approaches the performance of the fusion schemes as δ is increasing. The other upsampling method of [31] seems to produce large errors.
Computational efficiency is an important feature of any depth-stereo fusion method. Table V shows the execution times of the algorithms for the simulated data, as well as some of their implementation details. A combined Matlab-C version of our fusion algorithm requires 2.0s per image triplet while it takes 0.3s (for 2MP images) when GPU hardware is used for some initializations, and the SSE instruction set accounts for the online computation of correlation values between windows. We also developed a GPU-based implementation of a triangulation-based upsampling (interpolation) that takes less than 50ms. This allows one to envisage real-time execution of the proposed depth-stereo fusion framework, despite the high resolution. Note that it is the Matlab implementation that makes FastAgg and upsampling schemes slow, while the time for the NonLocalAgg method is based on matching at half resolution. We also point out that CSBP attains a solution in a reasonable time, despite its MRF-based formulation. 1.0 GC [52] > 10 3 CSBP [28] 15 FastAgg [25] > 10 3 NonLocalAgg [27] 5.0 Yang et al. [31] > 10 2 Kopf et al [34] 88 EPC [15] 1.5 F-MRF (GC) [12] > 10 3 F-MRF (ICM) [12] > 10 2 Our upsampling 95 F-ECC, F-EMCC 2.0 Upsampling of [15] < 0.1 * F-ECC, F-EMCC 0.3 * * time needed for the 2MP images of our real data-set (see Fig.13)

B. Real Data
Dal Mutto et al. [18] provide real TOF-stereo data along with ground truth disparities, shown in Fig. 12. To be consistent with [18], we upsample the depth in a similar way, using a bilateral filter, which benefits from color segmentation. This procedure is used to initialize the depth in all fusion schemes. Table VI shows the mean square error (MSE) of the disparity estimation for several algorithms. The contribution of the stereo data in fusion is unquestionable in scenes A and B, where the depth varies locally. All of the fusion methods obtain a more accurate map than the initial one. However, scene C contains only planar objects, and the upsampling methods provide good results. Although the proposed scheme does not perform best in all examples, it always improves the initial estimate, which demonstrates the advantage of adaptive fusion (very similar numbers are obtained with F-EMCC). As expected, the use of stereo data only (e.g., [25]) performs well only with the textured scene B. Our upsamping filter is more accurate than [34] and less accurate than the filter proposed by [18]. Note that our implementation of [34] achieves better results than those reported in [18].
We also assess the performance of the fusion methods on the HCIbox data-set [17]. The scene shows the interior of a box that contains some objects (Fig. 11). Note that there is no texture, apart from some horizontal lines on the stairs and the ramp, hence stereo methods tend to fail. We follow the experimental setup of [17], thus evaluating the depth estimation based on some statistics of the absolute error, after excluding inter-reflection areas (see Fig. 11). We do not include the results of [17], since the authors provided a different inter-reflection mask with larger support area than the one used in [17]. 6 Table VII shows the error statistics of the algorithms. All of the fusion methods start from the same initial map, obtained by our upsampling method. The proposed fusion method achieves the lowest mean and median error (similar results are obtained with F-EMCC). The variance of F-MRF is increased (a local bias was observed owing to the global smoothness), while its median remains low. Because of the depth discontinuities, [34] yields a less accurate result compared to our upsampling.
We also captured our own challenging TOF-stereo data-set using a synchronized camera setup, developed in collaboration with 4D View Solutions 7 . Two HR (1624 × 1224) color cameras and a MESA Imaging SR4000 TOF (range) camera The three scenes (cropped) of the dataset used in [18].
(176 × 144) are mounted on a rail. The stereo baseline is 50cm approximately while the TOF camera is mounted in the middle. The algorithm of [19] allows us to transform the depth-sensor measurements into a very sparse stereo disparity map (1% sparsity), with an average error of 1.0 pixel. The sparse map is refined and upsampled as discussed in Sec. IV (see also Fig.3). For a fully real-time upsampling, we use our GPU-based implementation of the triangulation-based interpolation [15]. Note that our stereo rectification puts the principal points in the same position, rather than making the optical axes parallel; this maximizes the overlap between the images, given the relatively wide baseline. Our 'MIXCAM' data-set contains challenging cases, e.g., periodicities, weakly textured areas, thin objects, depth discontinuities, and so on. The fusion algorithm merges the stereo and depth data and the outcome is a dense HR disparity map. We reuse our upsampling filter in a post-processing step to fill missing disparities. A streak-based method fills any remaining gaps in all algorithms. The results are shown in Fig. 13. White areas denote unmatched pixels, while black areas mark the detected TOF-occlusions. The left column shows the left image, with the TOF image shown in the bottom-right corner at the true scale. Next columns show the results of ELAS, FastAgg, F-MRF (GC) and F-ECC algorithms; the last column show the disparity maps of F-ECC after post-processing.
ELAS fills local areas, surrounded by textured points, through an interpolation scheme. We intentionally show the results of ELAS before the streak-based filling; as opposed to FastAgg, where missing disparities after the left-right consistency check are filled. Clearly, a pure stereo algorithm cannot deal with large untextured areas, and the post-filling is unreliable. F-MRF provides fully dense results. Note that we run F-MRF with half-resolution images (812 × 612), due  to its tremendous memory requirements. Moreover, we set a fixed value along the disparity range for the data-term of all TOF-occlusion points, so that the global inference becomes independent of this area. F-MRF provides artifacts in stereo occlusions, that are next to TOF-occlusions when the scene contains large foreground objects. As with [13], the results of F-MRF scheme verify the lack of an adaptive fusion of the depth-and stereo-consistency data terms, as opposed to our methods. However, F-MRF seems to deal better with very thin objects (e.g. the branch of the plant), as already discussed above. Note that the biased range measurements of very slanted surfaces (e.g. the table-top) negatively affect the fusion schemes, in particular when the table surface lacks texture (e.g. first example). The proposed scheme provides very good results on average, especially after the post-processing step, which fills the gaps and refines the disparities. We obtain very similar results with the F-EMCC method, while EPC provides results visually close to ours, but with more gaps. The bilateral upsampling of [34] provides visually good results, but with blurred depth discontinuities (see also Fig. 3).
VII. CONCLUSIONS We have presented a high-resolution stereo matching algorithm that is guided by low-resolution depth data, thus helping the algorithm to compensate for its difficulty in estimating disparities over weakly textured areas. We cast the problem into a MAP formulation whose inference is obtained through a series of local optimization problems, solved hierarchically in a seed-growing manner. The latter characteristic yields an intrinsically efficient solution that allows for near realtime matching of 2.0MP images. The data-term of the energy function benefits from a correlation function that is capable of providing scores at subpixel disparities, from an adaptive cost aggregation step inside the window based on the depth data, and from an adaptive fusion of stereo-and depth-consistency terms based on the scene texture and the camera geometry. These properties lead to a more selective growing process that prevents the algorithm from propagating incorrect disparities. As a result, a low-complexity method builds an accurate highresolution disparity map. A quantitative comparison against pure stereo and stereo-depth fusion algorithms, as well as a qualitative assessment on real data, has validated the strong performance of the proposed method. Future research will include the optimum visiting order for seeds in the growing framework, as well as an adaptive window size, based on the local surface orientation.