Patches, Planes and Probabilities:
A Non-local Prior for Volumetric 3D Reconstruction
Ali Osman Ulusoy
Michael J. Black
Andreas Geiger
Max Planck Institute for Intelligent Systems, Tübingen, Germany{osman.ulusoy,michael.black,andreas.geiger}@tue.mpg.de

Abstract
In this paper, we propose a non-local structured priorfor volumetric multi-view 3D reconstruction. Towards thisgoal, we present a novel Markov random field model basedon ray potentials in which assumptions about large 3D surface patches such as planarity or Manhattan world constraints can be efficiently encoded as probabilistic priors. We further derive an inference algorithm that reasons
jointly about voxels, pixels and image segments, and estimates marginal distributions of appearance, occupancy,
depth, normals and planarity. Key to tractable inferenceis a novel hybrid representation that spans both voxel andpixel space and that integrates non-local information from
2D image segmentations in a principled way. We compare
our non-local prior to commonly employed local smoothness assumptions and a variety of state-of-the-art volumetric reconstruction baselines on challenging outdoor scenes
with textureless and reflective surfaces. Our experiments
indicate that regularizing over larger distances has the potential to resolve ambiguities where local regularizers fail.

(a) Input Image(c) [50]+Pairwise Smoothness(b) Occupancy Probability [50](d) Our Result

Figure 1: Motivation: (a) The grass surface contains little texture, leading to reconstruction ambiguity. (b) Voxel
occupancy probabilities reveal this ambiguity where lighter
colors encode higher uncertainty [50]. (c) Pairwise smoothness priors cannot resolve the ambiguity and lead to a biased
3D model. (d) Our planarity prior regularizes over large distances and helps reconstruct the correct surface.

1. Introduction

caused by the textureless region. Ulusoy et al. [50] formulate 3D reconstruction as inference in a Markov random
field defined over the 3D voxel grid. Image evidence (input pixels) is modeled using ray potentials that accurately
incorporate visibility and free-space constraints.
While their method exposes reconstruction ambiguity,
it is not able to resolve this ambiguity to recover the correct surface because their model does not incorporate any
prior information; it models only image evidence. Luckily,
the 3D world we live in is not completely random but exhibits geometric structure. Previous works impose smoothness constraints via pairwise potentials that encourage adjacent voxels to take on the same occupancy state [18, 31] orcondition surface orientations on semantic information [40].

Dense 3D reconstruction from multiple RGB images isa long-standing problem in computer vision with numerous
practical applications. Unfortunately, it is also a highly illposed problem. Ambiguities arise in textureless areas orwhen photo-consistency assumptions are violated, e.g., atreflecting surfaces. For instance, consider the grass regionin Fig. 1a. The surface contains little texture, thus multiple
reconstructions satisfy the input images equally well.
Most previous work on multi-view stereo does not address such ambiguities and outputs a 3D model with no uncertainty information. In contrast, probabilistic approaches
model and expose the uncertainty in the reconstruction [1,5, 7, 8, 37, 50, 55]. Fig. 1b shows the result of a recent probabilistic method [50] that is able to expose the ambiguity
1

While these priors reduces surface noise to some extent,
they impose regularization only locally and are therefore
not sufficient to resolve large ambiguous regions as shownin Fig. 1c.
In this paper, we propose a novel prior formulation forvolumetric 3D reconstruction that encourages piece-wiseplanarity. We are inspired by the planar nature of many
elements in man-made environments, i.e., 3D range images of generic scenes can be approximated by piecewisesmooth regions with discontinuities at object boundaries[23]. Fig. 1d shows that our prior is able to disambiguate
large textureless regions to recover the correct surface.
Implementing a non-local prior in 3D is challenging.
Even in 2D, high-order spatial priors are expensive to represent and optimize [27, 33, 39]. Representing the planarity
prior directly in voxel space is complicated by the large variety of planes each single subvolume may contain, resulting
in numerous high-order cliques in the MRF.
Inspired by the success of non-local segmentation [25,29] and stereo matching [19, 32, 53] techniques, we encourage planarity within coherent image segments in all viewpoints. Our MRF reasons jointly about the occupancy andintensity of each voxel, the depth values observed at eachpixel, and the planarity and plane parameters in each image
segment. This hybrid 2D/3D representation with auxiliary
variables allows the inference algorithm to propagate viewbased planarity assumptions into 3D voxel space in a principled way and implicitly defines smoothness constraints oververy large neighborhoods in 3D.
The proposed MRF model is flexible in how plane priors can be integrated. In this work, we investigate a Manhattan world prior that encourages planes to align with thethree dominant orthogonal directions. Existing works onplanar multi-view stereo or Manhattan world representations treat these as hard constraints [13, 14, 32]. In contrast, we take a probabilistic approach where deviations
from the model are allowed as necessary (e.g., the sphereshaped dome in Fig. 1). Besides specifying the model,
we develop a message-passing algorithm for inferring approximate marginal distributions at every voxel, pixel andsegment. Our experiments demonstrate that the proposed
method improves upon state-of-the-art volumetric reconstruction techniques, in particular for challenging outdoor
scenes with large ambiguous (e.g., textureless) areas. Our
code and supplementary material are available at http://ps.
is.tue.mpg.de/research projects/volumetric-reconstruction.

2. Related WorkWe first review the most relevant work on probabilistic
volumetric reconstruction and then discuss approaches that
exploit primitives for scene modeling. For a more complete
review, we refer the reader to [15, 44].

Volumetric Reconstruction: Following early work [1, 7,
28, 43], Pollard and Mundy [37] propose a volumetric reconstruction method that updates the occupancy and colorof each voxel sequentially for each image. GPU implementations of this framework show impressive results [9, 49].
However, their framework lacks a global probabilistic formulation leading to evidence overcounting [38, 50]. To address this, a number of recent approaches have phrased 3Dvolumetric reconstruction as MRF inference, exploiting thespecial characteristics of high-order ray potentials to accurately model the image formation process [18, 31, 40, 50].
Reconstruction with Primitives: Several methods exploit
planar patches to represent piece-wise planar [3, 6, 10, 12,
34,41,45,46,56] or Manhattan world [13,14,42,47] scenes.
While the approaches produce impressive results, they enforce planarity as a hard constraint and thus only apply topiece-wise planar scenes. In contrast, our spatial prior canbe viewed as a soft constraint on planarity because it allows
deviations from planarity where it does not hold.
In a similar spirit, Häne et al. [22] propose a model forpiecewise planar depth map fusion. Their method takes as
input depth maps and a dictionary of patches and integratethese patches as soft constraints in a total variation framework that leads to improved results wrt. classical TV priors.
In contrast to the proposed hybrid pixel/voxel approach,their method is restricted to a 2.5D image representationand handles only very small patches (3 − 5 pixels) while weregularize over much larger regions (up to 10k pixels).
Gallup et al. [17] sample planes from initial depth maps
and exploit a semantic classifier to classify the image intoplanar and non-planar regions. Inference is then performed
via graph cuts, segmenting the image into regions explainedby planes and non-planar regions. Similar to [22], their
method uses a 2.5D image representation and requires depth
maps as input, while our approach integrates all constraints
into a single joint volumetric reconstruction and directly
takes RGB images as input.
Lafarge et al. [30] propose a method that simultaneously optimizes 3D primitives and a mesh using an objective that combines photo-consistency terms, mesh smoothness and priors on pairwise primitive arrangements. While
their method demonstrate impressive results, it is limited by
the topology and shape of the mesh initialization. Additionally, their method outputs a deterministic 3D model while
our approach yields a probabilistic 3D interpretation.
More recently, semantic and shape information has been
leveraged as prior knowledge for stereo matching [19] andmulti-view reconstruction [17, 20, 21, 35, 57], e.g., by constraining the set of plausible geometries [2, 11, 19] or bymodeling class specific normal distributions [20,21]. While
our focus in this work is on planarity and Manhattan worldpriors, semantic information can be easily integrated into
our framework.

3. Probabilistic Model
This section introduces our hybrid model for probabilistic volumetric 3D reconstruction with segment-based priors.
As input we assume a set of images and camera poses which
we obtain using structure-from-motion [51, 52]. As ourwork extends [50], we use their notation whenever possible.
To make this paper self-contained, we briefly repeat the image formation process described in [50] in Section 3.2. We
then specify our model in Section 3.3. Details about ourinference algorithm will be given in Section 4.

Ray Potentials & Occupancy Potentials

Patch Prior

3.1. Notation
The 3D space is decomposed into a grid of voxels. Eachvoxel is assigned a unique index from the index set X . We
associate each voxel i ∈ X with two random variables: abinary occupancy variable oi ∈ {0, 1} which signals if thevoxel is occupied (oi = 1) or free (oi = 0), and an appearance variable ai ∈ R describing the voxel intensity (or moregenerally, color).
Let R denote the set of viewing rays of all cameras. Notethat we model one viewing ray per pixel, thus R also corresponds to the total set of pixels. For a single ray r ∈ R,
let or = {or1 , . . . , orNr } and ar = {ar1 , . . . , arNr } denote theordered sets of occupancy and appearance variables associated with voxels intersecting ray r. The ordering is definedby the distance to the respective camera. We further associate each pixel/ray r ∈ R with an auxiliary depth variable
dr , discretized according to the depth of each voxel along
the ray r.
Each input image is segmented using the superpixelization algorithm of [54], yielding a set S that comprises all
segments from all input images. We associate each segments ∈ S with two random variables: a binary planarity variable ps ∈ {0, 1} indicating whether the segment is planar(ps = 1) or not (ps = 0), and a variable ns ∈ R3 specifying 3D plane parameters (i.e., xT ns = 1 if x ∈ R3 on
plane ns ) for this segment. We abbreviate the total set ofoccupancy and appearance variables in the voxel grid witho = {oi |i ∈ X } and a = {ai |i ∈ X }. We further summarize all depth, planarity and plane normal variables byd = {dr |r ∈ R}, p = {ps |s ∈ S} and n = {ns |s ∈ S},
respectively.

3.2. Image Formation
An image is formed by assigning each pixel the appearance of the first occupied voxel along ray r [50]:
Ir =

Nr
X
i=1

ori

Y

(1 − orj ) ari + 

(1)

j<i

where  ∼ N (0, σ) is a noise term, Ir denotes the intensity(or color) at the pixel corresponding to ray r. Note that the

Figure 2: Factor Graph. Our graphical model in plate notation. R comprises the set of all pixels/rays in all images
and S comprises all segments. Rs is the set of pixels inside
segment s. Ray depth variables dr connect voxels with segp
ments in our hybrid 3D/2D representation. Note that ψsrconnects only to those dr ’s for wich r ∈ Rs .
Q
term ori j<i (1 − orj ) evaluates to 1 for the first occupiedvoxel along the ray and to 0 for all other voxels.

3.3. Markov Random Field
We formulate volumetric 3D reconstruction as inference
in a Markov random field. We specify the joint distribution
over o, a, d, p and n asp(o, a, d, p, n) =Y
1 Y o
ϕi (oi )ψra (or , ar ) ψrd (or , dr )
Z
i∈X
r∈R
Y
Y
p
d
×
ϕs (ps ) ϕnψsr
(dr , ps , ns )s (ns )
s∈S

(2)

r∈Rs

where Z denotes the partition function, Rs contains all
rays/pixels associated with segment s, and ϕ and ψ denote
unary and high-order potentials, respectively. The corresponding factor graph is illustrated in Fig. 2 (top).
Voxel Occupancy Prior: We model our prior belief about
the state of the occupancy variables using a Bernoulli distribution
ϕoi (oi ) = γ oi (1 − γ)1−oi(3)where γ is the prior probability that voxel i is occupied.
Appearance Ray Potential: Our ray potentials model theimage generation process as specified by Eq. 1, i.e., theyencourage the appearance of the first occupied voxel along
ray r to agree with the image observation Ir at pixel r:ψra (or , ar ) =

Nr
X
i=1

ori

Y
j<i

(1 − orj ) νr (ari )(4)Here, νr (a) denotes the probability of observing intensity aat ray r which we model as νr (a) = N (a|Ir , σ).
Depth Ray Potential: The depth ray potential models theconstraint that the depth dr at pixel r must equal the depthof the first occupied voxel along the ray,(PNr r Qoi j<i (1 − orj ) dri
1 if dr = i=1dψr (or , dr ) =0 otherwise(5)
where dri denotes the depth of voxel i along ray r. Note thatthe potential ψrd (or , dr ) is 1 if and only if dr is equal to thedepth of the first occupied voxel along the ray. Otherwise,
the potential evaluates to 0, indicating an invalid state.
Planarity Prior: We now describe our planarity prior,which favors piece-wise planar depth maps in each of the
input views. As illustrated in Fig. 2 (right), the prior comprises the three components below.
Planarity Potential: We encourage segment planarity via
ϕps (ps ) = exp(λs |Rs | ps )(6)

where ps ∈ {0, 1} is the planarity indicator variable of segment s. Note that each segment comprises many pixels.
Thus, we weight this planarity potential by the segment area
|Rs |. The weight λs determines the strength of the prior,i.e., how much planarity should be enforced.
Plane Normal Potential: Our prior belief about the orientation of planar patches is modeled as a mixture distribution
on the plane normal ns
ϕns (ns ) =

K
X
k=1


wk M

ns
µk , κk
kns k


(7)

where M (· | µ, κ) denotes the von Mises-Fisher distribution with parameters µ (mean direction) and κ (concentration parameter). While any kind of surface orientation information can be incorporated into this potential (e.g., semantic information), in Section 6, we will investigate a rather
generic Manhattan world prior (K = 3, κ > 0).
Plane Depth Potential: To ensure that the plane parameters of a planar patch agree with the depth variables of thecorresponding image pixels, we define(exp (−λp η(dr − Dr (ns ))) if ps = 1pψsr (dr , ps , ns ) =1otherw.
(8)
where λp is a consistency weight, η(·) denotes a penalty
function and Dr (ns ) returns the depth of plane ns along rayr. The intuition behind this potential is simple: if ps = 0,then all possible depth values dr are equally likely for eachr ∈ Rs . In contrast, when ps = 1, we favor depth values

dr that are close to the plane ns . To account for segmentation errors, which can lead to outliers within a segment, we
model this factor using a robust Lorentzian penalty η(·).

4. Inference
We are interested in estimating the marginal distributionsof o, a, d, p and n in the proposed MRF. Given thesemarginals, we can easily calculate several quantities of interest, e.g., the probability of voxel occupancy and intensity, the probability of a segment s being planar and the expected plane normal, as well as the probability distributionof depth along each ray. Importantly, these depth distributions enable the computation of depth maps that are optimalin terms of Bayes decision theory [50].
Unfortunately, inference in our graphical model is challenging due to the large number of variables (o, a, d),the high-order potentials for modeling visibility constraints(ψra , ψrd ) and the mixed discrete (o, d, p) and continuous(a, n) state spaces of the variables. Furthermore, our factor
graph in Fig. 2 contains a large number of loops due to intersecting viewing rays R. Thus, exact inference is intractable.
In this section, we show how an approximation to the desired marginals can be obtained using message passing. In
particular, we derive an algorithm based on sum-productparticle belief propagation [24] in factor graphs [26].

4.1. Message Passing
Let µf →x denote the message sent from factor f to variable x, and let µx→f denote the corresponding variable-tofactor message. The messages from the unary factors to
the variables, as well as the variable-to-factor messages are
straightforward and we omit them here to save space. In
the following, we present the message equations for the appearance ray potential ψra , the depth ray potential ψrd andp. The supplementary document
the plane depth potential ψsrcontains detailed derivations of all the equations. Below, we
assume that all incoming messages to a factor are normalized such that they sum/integrate to 1.
Plane Depth Messages: The continuous message fromp
the plane depth potential ψsr
to the plane parameters ns isgiven by
XX
p
p
µψsr
ψsr
(dr , ps , ns )µ(ps )µ(dr ) (9)→ns (ns ) =ps= µ(ps = 1)

X

dr

p
ψsr
(dr , ps = 1, ns )µ(dr ) + µ(ps = 0)drwhere we have abbreviated the incoming messages usingµ(ps ) = µps →ψrd (ps ) and µ(dr ) = µdr →ψrd (dr ). Notepthat the message µψsr→ns (ns ) becomes uniform if there isstrong evidence of non-planarity, i.e., µ(ps = 0) = 1, fromthe other pixels in the segment, e.g., in case of a highlycurved surface. Otherwise, i.e. µ(ps = 1) = 1, the message evaluates high for planes ns that agree with the depth
distribution µ(dr ).
The message to the binary planarity variable ps reads as
Z X
p
p
ψsr
(dr , ps = 1, ns )µ(ns )µ(dr )µψsr →ps (ps = 1) =
ns d
r
p
µψsr
→ps (ps = 0) = 1(10)with µ(ns ) = µps →ψrd (ns ) and µ(dr ) = µdr →ψrd (dr ).
The planarity message is high if the depths of likely planes(where µ(ns ) is high) coincide with likely depths (highvalue of µ(dr )). Otherwise, the depth distribution at the
pixel cannot be explained with the incoming plane distribution, and therefore, the planarity message is low.
Finally, the message to the depth variable dr is given by
XZ
p
p
ψsr
(dr , ps , ns )µ(ps )µ(ns ) (11)µψsr →dr (dr ) =

j<i

p
ψsr
(dr , ps = 1, ns )µ(ns ) + µ(ps = 0)

= µ(ps = 1)ns

If there is strong evidence of non-planarity from the otherpixels in the segment, i.e. µ(ps = 0) = 1, then the message to the depth variable becomes uniform, i.e. the depth
variables are not affected. Otherwise, the message is highfor values dr that match the depth of likely planes, i.e.,dr ≈ Dr (ns ) where µ(ns ) is high.
Appearance Ray Messages: In a naı̈ve application ofbelief propagation, computing the factor-to-variable messages requires exponential time in the number of variables
involved in the potential. Since each viewing ray intersects
hundreds of voxels, the appearance ray potentials typicallyinvolve hundreds of variables, making this computation intractable. However, the special structure of the ray potentials allows for reducing the computation time from exponential to linear in the number of variables [31,50]. Exploiting this property, Ulusoy et al. [50] derived the sum-productmessage equations for this factor, which we include below
for completeness:XYµψra →ori (ori = 1) =µ(orj = 1)µ(ork = 0) ρrj
j<i

+

Y

k<j

µ(ork = 0) ρri

(12)

k<i

µψra →ori (ori = 0) =Xµ(orj = 1)

j<i

+

X
j>i

Depth Ray Messages: The message from ψrd to dr is readily given byYµψrd →dr (dr = dri ) ∝ µ(ori )µ(orj = 0)
(14)

ns

ps

Z

message has an intuitive interpretation: it increases the occupancy probability of voxels that are photo-consistent andvisible. The probability of voxels between the camera andthe likely surface location are decreased. For occluded voxels the message is uniform.
Ulusoy et al. also showed that the messages to the continuous appearance variables, i.e. µψra →ai , can be computed
analytically and that they can be compactly represented
as a constant plus a weighted Gaussian distribution. The
variable-to-factor messages µai →ψra cannot be computed
analytically. We follow [50] and approximate them using
Mixture-of-Gaussians (MoG) distributions. We refer thereader to [50] (Section 4.1) for all necessary details.

µ(orj = 1)Yµ(ork = 0) ρrj

where the incoming occupancy messages are again abbreviated by µ(ori ) = µori →ψrd (ori ).
By following a similar argument to the derivation of theappearance ray messages (see supplementary document),
the message to the occupancy variable ori is given byXYµψrd →ori (ori = 1) =µ(orj = 1)µ(ork = 0) µ(drj )
j<i

+

Y

k<j

µ(ork = 0) µ(dri )

(15)

k<i

µψrd →ori (ori = 0) =Xµ(orj = 1)

j>i

+

X
j<i

µ(orj = 1)Yµ(ork = 0) µ(drj )

k<j
k6=i

Y

µ(ork = 0) µ(drj ) (16)k<jwhere the incoming depth messages are abbreviated asµ(dri ) = µdr →ψrd (dr = dri ). Note that the messages
Eq. 15+16 are very similar to those of the appearance ray
factor in Eq. 12+13. The difference is that the photoconsistency measure ρ in Eq. 12+13 is replaced with theincoming depth message µ(d) which carries informationfrom our planarity prior. The messages in Eq. 15+16 intuitively increase the occupancy probability of voxels at likely
depths, i.e., where µ(drj ) is high. The probability of voxels between the camera and the likely depth are decreased
whereas the message to the occluded voxels are uniform.

k<jYµ(ork = 0) ρrj

(13)

k<j
k6=i

The incoming occupancy messages are abbreviated byµ(ori ) = µori →ψra (ori ) and ρrj is a photo-consistency measure at the jth voxel along ray r, see [50] for details. This

4.2. Particle Belief Propagation
Unfortunately, the continuous plane parameter variablesn complicate the message computations. In particular, Eq. 9
is of continuous form and therefore difficult to represent.
Further, the integrals that arise in Eq. 10+11 cannot be calculated in closed form.

1
K

K X
X
k=1 dr(17)(k)µ(ns )pψsr(dr , 1, ns )µ(dr )(k)
Ws (ns )Num pixels (%)

Num pixels (%)(a) C APITOL (all)
100
80
60
40
20
00

(k)

Given the set of plane particles {ns }K
k=1 , we build akernel density estimate in order to obtain an approximation
to the proposal distribution Ws (n)
K

(18)

k=1

where σkde denotes the kernel bandwidth, which we empirically set to a fixed value for all our experiments.

2

SP+pairwise
SP
LC+MoG
LC
PM
Ours

4
6
8
Error threshold (m)10(b) C APITOL (grass)

100
80
SP+pairwise
60
SP
LC+MoG
40
LC
PM
20
Ours
00
2
4
6
8
10
Error threshold (m)(c) BARUS &H OLLEY (all)where Ws (·) denotes an appropriate proposal distribution.
The approximation of the depth variable message in Eq. 11
is similar and can be found in the supplementary material.
Note that the quality of these approximations depend on
how well the set of particles explores the continuous space.
For instance, consider a perfectly planar patch. If none ofthe particles come close to the correct plane, the planarity
message in Eq. 10 will be underestimated, hence incorrectly
lowering the planarity belief of the segment.
To avoid this situation, in this work, we take advantage
of a data-driven strategy [48] that generates samples froman initial 3D reconstruction of the scene. More specifically,
we first run inference for our model without the planarity
patch prior and compute the most likely depth at each pixel.
We then repeatedly draw three pixels from the segment, prioritizing pixels with low depth certainty (negative entropy)
and fit a plane to them. Promising planes, i.e., planes that
explain most of the depth values inside the segment, areadded to the particle set. In addition to this particle set, we
also generate particles by conditioning on the plane prior
implied by the normal potentials. We draw plane normals
from the prior in Eq. 7 and then for each sample, optimize
the depth of the plane to best match the depth estimates
within the segment. Compared to a naı̈ve approach suchas uniformly sampling the space of plane parameters, ourdata-driven strategy avoids unlikely particles. Further details of our particle sampling strategy are presented in the
supplementary document.

1 X
Ws (n) =N (n; µ = n(k)s , Σ = σkde I)
K

100
80
SP+pairwise
60
SP
LC+MoG
40
LC
PM
20
Ours
00
2
4
6
8
10
Error threshold (m)Num pixels (%)

pµψsr→ps (ps = 1) ≈

100
80
60
40
20
00

Num pixels (%)

To tackle these challenges, we exploit particle belief
propagation [24], which has been adopted by many works
with great success [4, 19, 36, 53]. The main idea is to discretize the continuous space with a finite set of particles anduse this discretization to approximate the integral equations
with a Monte Carlo estimate.
(k)
For each segment, we draw K particles, {ns }K
k=1 , froma proposal distribution Ws (n). Using these particles, weapproximate the integral in Eq. 10 via importance sampling2

SP+pairwise
SP
LC+MoG
LC
PM
Ours

4
6
8
Error threshold (m)10(d) BARUS &H OLLEY (roof)

Figure 3: Reconstruction accuracy plots for the C APITOL
and BARUS &H OLLEY datasets.

5. Implementation
This section provides details of our implementation. We
initialize all ray messages uniformly and first pass messages
without the planarity prior, iterating over each image at least
once. This yields an initial 3D model. The beliefs in the occupancy variables are then propagated to the per-pixel depth
variables via Eq. 11. The median of these depth distributions yields a depth map which we use to segment the images using a depth-aware superpixelization algorithm [54].
We further use these depth maps to generate plane par(k)ticles {ns }K
k=1 for each segment as discussed in Section 4.2, using K = 64 throughout all our experiments.
Next, we compute the messages from the plane depth facp
to the plane parameters (Eq. 9) and to the planarity
tor ψsrvariables (Eq. 10). Finally, the information aggregated fromthe plane depth factor and all other pixels in the image segment is passed back to the depth variables using Eq. 11. The
depth variables in turn influence the occupancies along each
image ray via Eq. 15+16.
We interleave this process with the message updates forthe appearance ray factors and iterate over each image untilconvergence. All necessary algorithmic details can be foundin the supplementary document. Our implementation uses
GPU amenable octree data structures and GPU parallelization for message computations. Our current implementation
takes roughly 20 seconds to process a 1 Megapixel image
and 30 million voxels.

6. Experimental Evaluation
We evaluate our algorithm on three aerial datasets with
LIDAR ground truth1 provided by Restrepo et al. [38].
The datasets exhibit several challenges, including large fea1 The original datasets are distributed with sparse LIDAR point clouds.
Ulusoy et al. triangulated these point clouds to obtain a dense ground truth
mesh [50]. We use their meshes for a fair comparison to the baselines.

(a) Image(b) PM(c) LC(d) LC+MoG(e) SP(f) SP+pairwise(g) Ours(h) Our Depth

Figure 4: Qualitative Results. This figure shows qualitative results on the C APITOL (top), BARUS &H OLLEY (middle) andD OWNTOWN (bottom) dataset. (a) Region of interest. (b-g) Visualization of errors. Cooler colors correspond to lower error.
(h) Depth map predicted by our model. We encourage the reader to zoom in for details.
tureless regions, reflective and specular surfaces, severe
occlusions and transient objects. The input images areroughly one megapixel in size and each dataset contains∼ 200 images. The datasets are referred to as C APITOL,
BARUS &H OLLEY and D OWNTOWN, and example image
crops are presented in Fig. 4 (first column).
Baselines: We compare our results to several state-of-theart baselines. First, we compare to the sum-product algorithm of Ulusoy et al. [50], whose formulation is equivalent to our model without the patch prior and which wecall “SP” in the following. As their method does not encode any spatial regularization, we introduce an additional
baseline “SP+pairwise” with pairwise smoothness potentials that encourage adjacent voxels to take the same occupancy label [31] . We optimize the parameters of thispotential using the CAPITOL dataset. Further details canbe found in the supplementary document. Next, we compare to the approach of Liu and Cooper (“LC”) whose
max-product formulation utilizes ray potentials [31] as our
method, but suffers from a systematic bias in ambiguous regions as shown in [50]. We also include comparisons to an
improved version of Liu and Cooper’s approach as proposedby [50] with a more robust voxel color model, which we refer to as “LC+MoG”. Finally, we include a comparison to
Pollard and Mundy’s approach (“PM”) [37], which lacks aglobal probabilistic formulation but nevertheless achieves
very good results on several datasets [9, 49], including the
Middlebury benchmark [44].
Evaluation: We evaluate reconstruction accuracy by comparing depth maps predicted by each method to groundtruth. More specifically, we compute depth maps for all
input viewpoints and report the sum of per-pixel absolute
errors over all pixels in all views. We create ground truth
depth maps by projecting the LIDAR mesh onto the view.
For algorithms producing deterministic MAP outputs, i.e.,“LC” and “LC+MoG”, we consider the first occupied voxel
along the ray as a depth prediction. For algorithms which
compute occupancy marginals, we compute Bayes optimal
depth estimates at each pixel under the `1 loss as describedin [50] and the supplementary document.
We set the parameters of the superpixelization algorithm [54] to product roughly 500 segments. We found thissegmentation granularity to yield a reasonable tradoff between over- and undersegmentation of the images. We empirically chose a single set of parameters that we use for all
three datasets: λs = 5, λd = 1, κ = 20.
All three datasets contain many surfaces that are orthogonal to each other as can be seen in Fig. 4 (first column).
This structure motivates a Manhattan world prior on theplane orientations in Eq. 7, which favors planes oriented
along the X, Y or Z directions. However, this prior requires
orienting the scene such that the dominant directions in thescene coincide with the X, Y and Z direction. Towards
this goal, we first compute the ground plane normal (i.e.,
the Z direction) as the RANSAC fit to a sparse point cloud
generated by running our model without the planarity prior
and accumulating the 3D point clouds from each depth map.
The only remaining unknown is the rotation around the Z
axis, which we compute as the entropy minimizer of thepoint cloud projections onto canonical orthogonal X and Y
planes (see [16] for details).

Num pixels (%)

100
80
60
40
20
00

2

SP+pairwise
SP
LC+MoG
LC
PM
Ours

4
6
8
Error threshold (m)10(a) D OWNTOWN (all)(b) Image(c) SP(d) SP+pairwise(e) Ours

Figure 5: (5a) Accuracy plots for DOWNTOWN, (5b) Regionof interest, (5c-5e) Visualization of errors.
Discussion of Results: Fig. 3 and Fig. 5a present cumulative accuracy plots for the three datasets and all methods.
Overall our algorithm outperforms previous methods for the
C APITOL and D OWNTOWN datasets and achieves similar
performance for BARUS &H OLLEY.
We visualize the error maps for all datasets and methods
in Fig. 4. For the C APITOL and BARUS &H OLLEY datasets(first two rows), the majority of errors are localized on thelarge featureless regions, e.g., the grass region for C API TOL and the black rooftop in BARUS &H OLLEY . The results show that “PM” as well as algorithms that estimate a
MAP solution, i.e., “LC” and “LC+MoG”, yield a systematic bias in textureless regions, leading to large errors. The
sum-product (“SP”) result is able to reveal the ambiguity inthe region (see Fig. 1b) and obtains better results. Nonetheless, “SP” cannot resolve the ambiguity without additional
prior information. The model with the pairwise smoothness potentials (“SP+pairwise”) yields denser and smootherresults, leading to lower errors on most structures, e.g. thebuilding facades in C APITOL and BARUS &H OLLEY. However, results on the grass region in C APITOL and the rooftop
in BARUS &H OLLEY indicate that the pairwise model is unable to resolve ambiguities on the large textureless regions.
In contrast, our algorithm achieves significantly more accurate results in these regions as seen in Fig. 4g. To investigate
this more closely, we evaluate errors only on the large featureless surfaces. The results are displayed in Fig. 3b+3d
and show the clear improvement of our algorithm. Fig. 1d
depicts a rendering of our reconstruction where the ambiguity in the model is resolved successfully.
Besides textureless regions, our planarity prior helps im-prove reconstruction of reflective surfaces as well. The lastrow of Fig. 4a displays a building with a reflective surface
from DOWNTOWN. All baselines yield large errors on thebuilding surface since it violates the Lambertian surface assumption. Our planarity prior helps correct these errors toachieve a denser and more accurate result. Fig. 5 displays
results for another reflective building. All algorithms produce large errors for the side surface of the building which
has a mirror like reflectivity. However, the front side of thebuilding contains a mixture of reflective and non-reflective
materials. While all algorithms produce some correct depth
values in this region, the results are in general very noisy
and contain large holes. The model with pairwise smoothness potentials (“SP+pairwise”) fails to bridge these large
gaps in the reconstruction. In contrast, our algorithm regularizes over the entire facade and is able to reconstruct thisfacade surface correctly.
The results confirm that even though our prior favors
piecewise planar reconstructions, non-planar structures arepreserved, e.g. the building dome in Fig. 4 (first row) shows
very little error. Our algorithm is robust to non-planar structures due to two facts: First, our model allows for turning
off the planarity prior wherever necessary. Second, we also
allow for outliers within a segment, owing to our robust
penalty function η(·) in Eq. 8. The latter property is important to deal with imprecise segmentation boundaries.
For BARUS &H OLLEY, our planarity priors introduces
errors around some of the shrubs and trees, which lowers theoverall accuracy (Fig. 3c). However, for some of these regions, the LIDAR ground truth is not accurate either due tothe tree tops which have been extruded to the ground level.
In the supplementary document, we present an additional
evaluation excluding such regions, as well as further examples and failure cases.

7. Conclusion
We have presented a novel non-local prior for probabilistic volumetric 3D reconstruction that encourages planarity
within image segments and regularizes over large voxel
neighborhoods. Our experiments show that the proposed
prior is able to resolve reconstruction ambiguities of textureless and partially reflective surfaces and achieves stateof-the-art results in reconstruction accuracy for highly challenging aerial datasets. In our future work, we plan to incorporate semantic information in our model. Furthermore, ourprior formulation also allows for geometries beyond planar
segments. We believe that integrating more complex primitives such as spheres, cylinders, deformable shapes or evenwholistic 3D scene models will be promising extensions of
the presented model.

References[1] M. Agrawal and L. S. Davis. A probabilistic framework forsurface reconstruction from multiple images. In Proc. IEEE
Conf. on Computer Vision and Pattern Recognition (CVPR),
2001. 1, 2
[2] S. Bao, M. Chandraker, Y. Lin, and S. Savarese. Dense object
reconstruction with semantic priors. In Proc. IEEE Conf. on
Computer Vision and Pattern Recognition (CVPR), 2013. 2[3] A. Bartoli. A random sampling strategy for piecewise planar scene segmentation. Computer Vision and Image Understanding (CVIU), 105(1):42–59, 2007. 2[4] F. Besse, C. Rother, A. Fitzgibbon, and J. Kautz. PMBP:
PatchMatch Belief Propagation for correspondence field estimation. International Journal of Computer Vision (IJCV),
110(1):2–13, 2014. 6[5] R. Bhotika, D. J. Fleet, and K. N. Kutulakos. A probabilistic
theory of occupancy and emptiness. In Proc. of the European
Conf. on Computer Vision (ECCV), 2002. 1[6] A. Bodis-Szomoru, H. Riemenschneider, and L. Van Gool.
Fast, approximate piecewise-planar modeling based onsparse structure-from-motion and superpixels. In Proc. IEEE
Conf. on Computer Vision and Pattern Recognition (CVPR),
2014. 2[7] J. D. Bonet and P. Viola. Poxels: Probabilistic voxelized
volume reconstruction. In Proc. of the IEEE International
Conf. on Computer Vision (ICCV), 1999. 1, 2
[8] A. Broadhurst, T. W. Drummond, and R. Cipolla. A probabilistic framework for space carving. ICCV, 2001. 1[9] F. Calakli, A. O. Ulusoy, M. I. Restrepo, G. Taubin, and J. L.
Mundy. High resolution surface reconstruction from multiview aerial imagery. In 3DIMPVT, 2012. 2, 7
[10] A.-L. Chauve, P. Labatut, and J.-P. Pons. Robust piecewiseplanar 3d reconstruction and completion from large-scale unstructured point data. In Proc. IEEE Conf. on Computer Vision and Pattern Recognition (CVPR), 2010. 2[11] A. Dame, V. Prisacariu, C. Ren, and I. Reid. Dense reconstruction using 3D object shape priors. In Proc. IEEE Conf.
on Computer Vision and Pattern Recognition (CVPR), 2013.
2[12] F. Fraundorfer, K. Schindler, and H. Bischof. Piecewise planar scene reconstruction from sparse correspondences. Image and Vision Computing (IVC), 24(4):395–406, 2006. 2[13] Y. Furukawa, B. Curless, S. M. Seitz, and R. Szeliski.
Manhattan-world stereo. In Proc. IEEE Conf. on Computer
Vision and Pattern Recognition (CVPR), 2009. 2[14] Y. Furukawa, B. Curless, S. M. Seitz, and R. Szeliski. Reconstructing building interiors from images. In Proc. of theIEEE International Conf. on Computer Vision (ICCV), 2009.
2[15] Y. Furukawa and C. Hernandez. Multi-view stereo: A tutorial. Foundations and Trends in Computer Graphics andVision, 9(1-2):1–148, 2013. 2[16] D. Gallup, J.-M. Frahm, P. Mordohai, Q. Yang, and M. Pollefeys. Real-time plane-sweeping stereo with multiple sweeping directions. In Proc. IEEE Conf. on Computer Vision and
Pattern Recognition (CVPR), 2007. 7[17] D. Gallup, J.-M. Frahm, and M. Pollefeys. Piecewise planar and non-planar stereo for urban scene reconstruction. In
Proc. IEEE Conf. on Computer Vision and Pattern Recognition (CVPR), 2010. 2[18] P. Gargallo, P. Sturm, and S. Pujades. An occupancy–depth
generative model of multi-view images. In Proc. of the Asian
Conf. on Computer Vision (ACCV), pages 373–383, 2007. 1,
2[19] F. Güney and A. Geiger. Displets: Resolving stereo ambiguities using object knowledge. In Proc. IEEE Conf. on Computer Vision and Pattern Recognition (CVPR), 2015. 2, 6
[20] C. Haene, N. Savinov, and M. Pollefeys. Class specific 3d
object shape priors using surface normals. In Proc. IEEE
Conf. on Computer Vision and Pattern Recognition (CVPR),
2014. 2[21] C. Haene, C. Zach, A. Cohen, R. Angst, and M. Pollefeys.
Joint 3D scene reconstruction and class segmentation. In
Proc. IEEE Conf. on Computer Vision and Pattern Recognition (CVPR), 2013. 2[22] C. Haene, C. Zach, B. Zeisl, and M. Pollefeys. A patchprior for dense 3d reconstruction in man-made environments.
In Proc. of the International Conf. on 3D Digital Imaging,
Modeling, Data Processing, Visualization and Transmission(THREEDIMPVT), 2012. 2[23] J. Huang, A. B. Lee, and D. Mumford. Statistics of range
images. In Proc. IEEE Conf. on Computer Vision and Pattern
Recognition (CVPR), 2000. 2[24] A. Ihler and D. McAllester. Particle belief propagation. Conference on Artificial Intelligence and Statistics (AISTATS),2009. 4, 6
[25] P. Kohli, L. Ladicky, and P. H. S. Torr. Robust higher order potentials for enforcing label consistency. International
Journal of Computer Vision (IJCV), 82(3):302–324, 2009. 2[26] D. Koller and N. Friedman. Probabilistic Graphical Models:Principles and Techniques. MIT Press, 2009. 4[27] N. Komodakis and N. Paragios. Beyond pairwise energies:
Efficient optimization for higher-order MRFs. In Proc. IEEE
Conf. on Computer Vision and Pattern Recognition (CVPR),
2009. 2[28] K. N. Kutulakos and S. M. Seitz. A theory of shape by spacecarving. International Journal of Computer Vision (IJCV),
38(3):199–218, 2000. 2[29] L. Ladicky, C. Russell, P. Kohli, and P. H. S. Torr. Associative hierarchical random fields. IEEE Trans. on Pattern Analysis and Machine Intelligence (PAMI), 36(6):1056–
1077, 2014. 2[30] F. Lafarge, R. Keriven, M. Bredif, and H.-H. Vu. A hybrid
multiview stereo algorithm for modeling urban scenes. IEEE
Trans. on Pattern Analysis and Machine Intelligence (PAMI),
35(1):5–17, 2013. 2[31] S. Liu and D. Cooper. Statistical inverse ray tracing forimage-based 3d modeling. IEEE Trans. on Pattern Analysis
and Machine Intelligence (PAMI), 36(10):2074–2088, 2014.
1, 2, 5, 7[32] W. Luo and H. Maitre. Using surface model to correct and fit
disparity data in stereo vision. In Proc. of the International
Conf. on Pattern Recognition (ICPR), 1990. 2

[33] P. Marquez-Neila, P. Kohli, C. Rother, and L. Baumela. Nonparametric higher-order random fields for image segmentation. In Proc. of the European Conf. on Computer Vision(ECCV), 2014. 2[34] B. Micusik and J. Kosecka. Multi-view superpixel stereo
in urban environments. International Journal of Computer
Vision (IJCV), 89(1):106–119, 2010. 2[35] A. Owens, J. Xiao, A. Torralba, and W. T. Freeman. Shapeanchors for data-driven multi-view reconstruction. In Proc.
of the IEEE International Conf. on Computer Vision (ICCV),
2013. 2[36] J. Pacheco, S. Zuffi, M. J. Black, and E. Sudderth. Preserving modes and messages via diverse particle selection.
In Proc. of the International Conf. on Machine learning(ICML), 2014. 6[37] T. Pollard and J. L. Mundy. Change detection in a 3-d world.
In Proc. IEEE Conf. on Computer Vision and Pattern Recognition (CVPR), 2007. 1, 2, 7
[38] M. I. Restrepo. Characterization of Probabilistic Volumetric
Models for 3-d Computer Vision. PhD thesis, Brown University, 2013. 2, 6
[39] S. Roth and M. J. Black. Fields of experts. International
Journal of Computer Vision (IJCV), 82(2):205–229, 2009. 2[40] N. Savinov, L. Ladicky, C. Häne, and M. Pollefeys. Discrete
optimization of ray potentials for semantic 3d reconstruction.
In Proc. IEEE Conf. on Computer Vision and Pattern Recognition (CVPR), 2015. 1, 2[41] A. Saxena, M. Sun, and A. Y. Ng. Make3D: learning 3D
scene structure from a single still image. IEEE Trans. on
Pattern Analysis and Machine Intelligence (PAMI), 31:824–
840, 2009. 2[42] M. Schönbein and A. Geiger. Omnidirectional 3d reconstruction in augmented manhattan worlds. In Proc. IEEE International Conf. on Intelligent Robots and Systems (IROS), 2014.
2[43] S. Seitz and C. Dyer. Photorealistic scene reconstruction by
voxel coloring. Proc. IEEE Conf. on Computer Vision and
Pattern Recognition (CVPR), 1997. 2[44] S. M. Seitz, B. Curless, J. Diebel, D. Scharstein, andR. Szeliski. A comparison and evaluation of multi-view
stereo reconstruction algorithms. In Proc. IEEE Conf. on
Computer Vision and Pattern Recognition (CVPR), 2006. 2,
7[45] S. N. Sinha, D. Steedly, and R. Szeliski. Piecewise planarstereo for image-based rendering. In Proc. of the IEEE International Conf. on Computer Vision (ICCV), 2009. 2[46] N. Srinivasan and F. Dellaert. A rao-blackwellized mcmc
algorithm for recovering piecewise planar 3d models frommultiple view rgbd images. In Proc. IEEE International
Conf. on Image Processing (ICIP), 2014. 2[47] J. Straub, G. Rosman, O. Freifeld, J. J. Leonard, and J. W.
Fisher. A mixture of manhattan frames : Beyond the manhattan world. In Proc. IEEE Conf. on Computer Vision and
Pattern Recognition (CVPR), 2014. 2[48] Z. Tu, S. C. Zhu, and H. Shum. Image segmentation by data
driven markov chain monte carlo. In Proc. of the IEEE International Conf. on Computer Vision (ICCV), 2001. 6[49] A. O. Ulusoy, O. Biris, and J. L. Mundy. Dynamic probabilistic volumetric models. In Proc. of the IEEE International Conf. on Computer Vision (ICCV), 2013. 2, 7[50] A. O. Ulusoy, A. Geiger, and M. J. Black. Towards probabilistic volumetric reconstruction using ray potentials. In
Proc. of the International Conf. on 3D Vision (3DV), 2015.
1, 2, 3, 4, 5, 6, 7
[51] C. Wu. Towards linear-time incremental structure from motion. Proc. of the International Conf. on 3D Vision (3DV),
2013. 3[52] C. Wu, S. Agarwal, B. Curless, and S. Seitz. Multicore bundle adjustment. In Proc. IEEE Conf. on Computer Vision and
Pattern Recognition (CVPR), 2011. 3[53] K. Yamaguchi, T. Hazan, D. McAllester, and R. Urtasun.
Continuous markov random fields for robust stereo estimation. In Proc. of the European Conf. on Computer Vision(ECCV), 2012. 2, 6[54] K. Yamaguchi, D. McAllester, and R. Urtasun. Robust
monocular epipolar flow estimation. In Proc. IEEE Conf.
on Computer Vision and Pattern Recognition (CVPR), 2013.
3, 6, 7[55] A. Yao and A. Calway. Dense 3-d structure from image sequences using probabilistic depth carving. In Proc. of the
British Machine Vision Conf. (BMVC), 2003. 1[56] L. Zebedin, J. Bauer, K. Karner, and H. Bischof. Fusion offeature-and area-based information for urban buildings modeling from aerial imagery. In Proc. of the European Conf. on
Computer Vision (ECCV), 2008. 2[57] C. Zhou, F. Güney, Y. Wang, and A. Geiger. Exploiting object similarity in 3d reconstruction. In Proc. of the IEEE
International Conf. on Computer Vision (ICCV), 2015. 2