arXiv:1710.01749v1 [cs.CV] 4 Oct 2017

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

1

Semantic 3D Reconstruction with Finite
Element Bases
Audrey Richard1,†

1

Photogrammetry & Remote Sensing
ETH Zurich, Switzerland

2

Institute of Computer Graphics & Vision
TU Graz, Austria

3

Austrian Institute of Technology

audrey.richard@geod.baug.ethz.ch

Christoph Vogel2,†
christoph.vogel@icg.tugraz.at

Maroš Bláha1
maros.blaha@geod.baug.ethz.ch

Thomas Pock2,3
thomas.pock@icg.tugraz.at

Konrad Schindler1

konrad.schindler@geod.baug.ethz.ch

1

†

shared first authorship

Abstract
We propose a novel framework for the discretisation of multi-label problems on arbitrary, continuous domains. Our work bridges the gap between general FEM discretisations, and labeling problems that arise in a variety of computer vision tasks, including
for instance those derived from the generalised Potts model. Starting from the popular
formulation of labeling as a convex relaxation by functional lifting, we show that FEM
discretisation is valid for the most general case, where the regulariser is anisotropic and
non-metric. While our findings are generic and applicable to different vision problems,
we demonstrate their practical implementation in the context of semantic 3D reconstruction, where such regularisers have proved particularly beneficial. The proposed FEM
approach leads to a smaller memory footprint as well as faster computation, and it constitutes a very simple way to enable variable, adaptive resolution within the same model.

Introduction

A number of computer vision tasks, such as segmentation, multiview reconstruction, stitching and inpainting, can be formulated as multi-label problems on continuous domains, by
functional lifting [10, 13, 28, 34, 36]. A recent example is semantic 3D reconstruction (e.g.
[4, 18]), which solves the following problem: Given a set of images of a scene, reconstruct
both its 3D shape and a segmentation into semantic object classes. The task is particularly
challenging, because the evidence is irregularly distributed in the 3D domain; but it also
possesses a rich, anisotropic prior structure that can be exploited. Jointly reasoning about
shape and class allows one to take into account class-specific shape priors (e.g., building
walls should be smooth and vertical, and vice versa smooth, vertical surfaces are likely to
be building walls), leading to improved reconstruction results. So far, models for the mentioned multi-label problems, and in particular for semantic 3D reconstruction, have been
limited to axis-aligned discretisations. Unless the scenes are aligned with the coordinate
c 2017. The copyright of this document resides with its authors.
It may be distributed unchanged freely in print or electronic forms.

2

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

Figure 1: Semantic 3D model, estimated from aerial views with our FEM method.
axes, this leads to an unnecessarily large number of elements. Moreover, since the evidence
is (inevitably) distributed unevenly in 3D, it also causes biased reconstructions. Thus, it is
desirable to adapt the discretisation to the scene content (as often done for purely geometric
surface reconstruction, e.g. [25]).
Our formulation makes it possible to employ a finer tesselation in regions that are likely
to contain a surface, exploiting the fact that both high spatial resolution and high numerical
precision are only required in those regions. Our discretisation scheme leads to a smaller
memory footprint and faster computation, and it constitues a very simple technique to allow for arbitrary adaptive resolution levels within the same problem. I.e., we can refine or
coarsen the discretisation as appropriate, to adapt to the scene to be reconstructed. While
our scheme is applicable to a whole family of finite element bases, we investigate two particularly interesting cases: Lagrange (P1) and Raviart-Thomas elements of first order. We
further show that the grid-based voxel discretisation is a special case of our P1 basis, such
that minimum energy solutions of “identical” discretisations (same vertex set) are equivalent.

2

Related Work

Since the seminal work [14] volumetric reconstruction from image data has evolved remarkably [12, 16, 21, 22, 23, 29, 43, 46]. Most methods use depth maps or 2.5D range scans
for evidence [51, 52], represent the scene via an indicator or signed distance function in the
volumetric domain, and extract the surface as its zero level set, e.g., [30, 42].
Joint estimation of geometry and semantic labels, which had earlier been attempted only
for single depth maps [26], has recently emerged as a powerful extension of volumetric 3D
reconstruction from multiple views [1, 4, 18, 24, 39, 44, 45]. A common trait of these
works is the integration of depth estimates and appearance-based labeling information from
multiple images, with class-specific regularisation via shape priors.
Multi-label problems are in general NP-hard, but under certain conditions on the pairwise interactions, the original non-convex problem can be converted into a convex one via
functional lifting and subsequent relaxation, e.g. [10]. This construction was further extended to anisotropic (direction-dependent) regularisers [41]. Moreover, [53] also relaxed
the requirement that the regulariser forms a metric on the label set, yet its construction can
only be applied after discretisation [28]. In this paper, we consider the relaxation in its most
general form [53], but are not restricted to it. The latter construction is also the basis to the
model of [18], whose energy model we adapt for our semantic 3D reconstruction method.
Their voxel-based formulation can be seen as a special case of our discretisation scheme.
For (non-semantic) surface reconstruction, several authors prefer a data-dependent discretisation, normally a Delaunay tetrahedralisation of a 3D point cloud [20, 25, 47]. The occupancy states of the tetrahedra are found by discrete (binary) labeling, and the final surface
is composed of the triangles that separate different labels. Loosely speaking, our proposed
methodology can be seen either as an extension of [18] to arbitrary simplex partitions of the
domain; or as an extension of [25] to semantic (multi-label) reconstruction.

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

3

We note that regular voxel partitioning of the volume leads to a high memory consumption and computation time. Yet, we are essentially reconstructing a 2D manifold in 3D space,
and this can be exploited to reduce run-time and memory footprint. [24] use an octree instead of equally sized voxels to adapt to the local density of the input data. [4] go one step
further and propose an adaptive octree, where the discretisation is refined on-the-fly, during
optimisation. In our framework the energy is independent of the discretisation, it can thus be
combined directly with such an adaptive procedure.
Also volumetric fusion via signed distance functions [32] benefits from irregular tesselations of 3D space, e.g., octrees [40] or hashmaps [33]. In contrast to our work, these
target real-time reconstruction and refrain from global optimisation, instead locally fusing
depth maps. Their input normally is a densely sampled, overcomplete RGB-D video-stream,
whereas we deal with noisy and incomplete inputs. To achieve high-quality reconstructions
in our setting, we incorporate semantic information, leading to a multi-label problem.
Our work is based on the finite element method (FEM), e.g. [8, 37]. Introduced by Ritz
[38] more than a century ago, and refined by Galerkin and Courant [11], FEM serves to
numerically solve variational problems, by partitioning the domain into finite, parametrised
elements. In computer vision FEM has been applied in the context of level-set methods [49]
and for Total Variation [2]. To our knowledge, we are the first to apply it to multi-labeling.

3

Method

The multi-labeling problem [10, 28, 41, 53] in the domain Ω ⊂ Rd is defined by finding m
labeling functions xi : Ω → {0, 1}, i = 1 . . . m as the solution of:
m

inf ∑

Z

xi i=1 Ω

m

ρ i (z)xi (z)dz + J(xi ),

s. t.

∑ xi (z) = 1 ∀z ∈ Ω,

(1)

i=1

where ρ models the data term for a specific label at location z ∈ Ω and J denotes a convex
regularisation functional that enforces the spatial consistency of the labels. One prominent
example is to chose J := k·k 2 , known as Total Variation, which penalises the perimeter of the
individual regions [10, 34]. Note that in the two-label case (Potts model), this relaxation is
exact after thresholding with any threshold from the open unit interval [10]. Although we are
ultimately interested in non-metric regularisation, we start with the continuous, anisotropic
model [41], and postpone the extension to the non-metric case to Sec. 3.5.

3.1

Convex Relaxation

The continuous model allows for an anisotropic regulariser in J: label transitions can be
penalised on the area of the shared surface, as well as on the surface normal direction. This
is achieved with problem-specific 1-homogeneous functions that emerge from convex sets,
so called Wulff-shapes. A relaxation of xi (z) ∈ {0, 1} to xi (z) ∈ [0, 1] then leads to a convex
energy, which can be written as the following saddle point problem, with primal functions x
and dual functions λ :
Z

m

i Ω

i=1

min max ∑ ρ i(z)xi (z)+hxi (z), ∇·λ i (z)idz, s.t. λ i (z)−λ j (z) ∈ W i j, ∑ xi (z)=1, xi (z)≥0. (2)
xi

λi

The constraints have to be fulfilled for all z ∈ Ω. In addition to the primal variables xi ,
we have introduced the dual vector-field λ i : Rd → Rd , whose pairwise differences are
constrained to lie in the convex sets (Wulff-shapes) W i j . By letting these shapes take an
anisotropic form, one can then encode scene structure, e.g. [18, 41]. For our problem we

4

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

demand Neumann conditions at the boundary of Ω, i.e. hλ i , νi = 0, ∀z ∈ ∂ Ω, because the
scene will continue beyond our domain (ν is the normal of the domain boundary ∂ Ω).

3.2

Finite Element Spaces

Here, we can only informally introduce the basic idea of FEM and explain its suitability for
problems of the form (2). We refer to textbooks [8, 15, 27] for a deeper and formal treatment.
One way to solve (2) is to approximate it at a finite number of regular grid points, using
finite differences. FEM instead searches for a solution in a finite-dimensional vector space;
this trial space is a subspace of the space in which the exact solution is defined. To that
end, one chooses a suitable basis for the trial space, with basis functions of finite support, as
well as an appropriate test function space. FEM methods then find approximate solutions to
variational problems by identifying the element from the trial space that is orthogonal to all
functions of the test function space. For our saddle-point problem, we can instead identify
the trial space with our primal function space and the test space with its dual counterpart.
Now, we can apply the same principles, and after discretisation our solution corresponds to
the continuous solution defined by the respective basis. As (2) is already a relaxation of the
original problem (1), we do not present an analysis of convergence at this point. Instead, the
reader is referred to [2] for an introduction to this somewhat involved topic.
In order to choose a space with good approximation properties and suitable basis functions, we tesselate our domain into simplices. More formally, we define M = {F,V, S} to be
a simplex mesh with vertices v ∈ V, v ∈ Rd , faces f ∈ F defined by d, and simplices s ∈ S,
defined by d+1 vertices that partition Ω: ∪k sk = Ω, sl ∩ sk = fl,k ∈ F – i.e. two adjacent simplices share only a single face. In this work, for a specific set of vertices V , we select M to
be the corresponding Delaunay tetrahedralisation of Ω and only consider explicit bases. In
particular, we focus on the Lagrangian (P1) basis, which we use in the following to derive
our framework; and on the Raviart-Thomas (RT) basis. Details for the latter are given in the
supplementary material. The main difference between them is that P1 leads to piecewise linear solutions, which must be thresholded, while RT leads to a constant labeling function per
simplex, similar to discrete MAP solutions on CRFs. We note that constant labeling can lead
to artefacts, such that the adaptiveness of the FEM model becomes even more important.
The idea of both derivations is similar: (i) select a basis for our primal (P1) or dual (RT)
variable set, (ii) find a suitable form via the divergence theorem and Fenchel duality, (iii)
extend to the non-metric case, following a principle we term ”label mass preservation”.

3.3

Lagrange Elements

The Lagrange Pk (M) basis functions describe a conforming polynomial basis of order k+1
on our simplex mesh M, i.e. its elements belong to the Hilbert space of differentiable function with finite Lebesgue measure on the domain Ω: Pk (M) ⊂ H 1 (Ω) := {p ∈ L2 (Ω), ∇p ∈
(L2 (Ω)d )}. We are interested in the Lagrange basis of first order, P1 (M):
P1 (M):= {p : Ω → R|p ∈C(Ω), p(x) := ∑ cTs x + ds , cs ∈ Rd , ds ∈ R, if x ∈ s and 0 else}. (3)
s∈S

We construct our linear basis with functions that are defined for each vertex v of a simplex s
and can be described in local form with barycentric coordinates:
p1s,v (x) := αv with x = ∑ αv v, ∑ αv = 1, αv ≥ 0 if x ∈ s and 0 else.
(4)
v∈s

v∈s

In each simplex, one can define a scalar field φs (x) ∈ R and compute a gradient in this basis
that will be constant per simplex s (cf. Fig. 2):

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

5

Figure 2: Left: Illustration of P1 basis function shape. Middle: Scalar field defined as a
convex combination of basis coefficients. Right: Gradient definition in a simplex (5).
φs (x) := ∑ φv p1s,v and ∇φs = ∑ φv Jv ,
v∈s

(5)

v∈s

with coefficients φv ∈ R. Jv ∈ Rd denotes a vector that is normal to the face fv opposite node
| fv |
, and points towards the simplex centre. | fv | denotes the area of the face fv ,
v, has length |s|d
and |s| the volume of the simplex s (cf . Fig. 2 and supplementary material).

3.4

Discretisation

To apply our Lagrange basis to (2) we first make use of the divergence theorem:
Z

hxi (z), ∇·λ i (z)idz =

Ω

Z

h∇xi (z), λ i (z)idz −

Z

Ω

∂Ω

|

hν(z), λ i (z)idz .
{z
}

(6)

=0

The latter summand vanishes by our choice of λ . Our approach for a discretisation in the
Lagrange basis is to choose the labeling function xi ∈ P1 (M). This implies that our dual
space consists of constant vector-fields per simplex: λsi ∈ Rd . To fulfill the constraint set in
(2) we have to verify that, per simplex, the λsi lie in the respective Wulff-shape. The simplex
constraints on the xi have to be modeled per vertex. According to (4), the labeling functions
are convex combinations of their values at the vertices and thus stay within the simplex.
We also have to convert the continuous data costs ρ i into a cost per vertex ρvi , which
can be achieved by convolving the continuous cost with the respective basis function: ρvi :=
R
i
Ω ∑s∈N (v) φs (x)ρ (x)dx. In practice, the integral can be computed by sampling ρ. Integrating the right hand side in (6) over the simplex s leads to a weighting with its volume |s| and
the energy (2) in the discrete setting becomes:
m

min max ∑ ρvi xvi +∑ |s|h∇xi , λsi i s.t. (λsi − λsj)∈W i j ∀i< j, s ∈ S, ∑ xvi =1, xvi ≥0 ∀v∈V. (7)
xi

λi

3.5

v,i

s,i

i=1

Non-metric extension

To start with, we note that a non-metric model does not exist in the continuous case [31]
and our extension works only after the discretisation into the FEM basis. Please refer to
the supplementary material for an in-depth discussion. Note that our label set of semantic
classes does not have a natural order (in contrast to, e.g., stereo depth or denoised brightness);
and also the direction-dependent regulariser is unordered and does not induce a metric cost.
To allow for non-metric regularisation we transform the constraint set (λ i − λ j ∈ W i j ), by
introducing auxiliary variables zi j and Lagrange multipliers yi j , and use Fenchel-Duality:
max min ∑ h(λsi − λsj ) − zi j , yisj i − δW i j (zisj ) = max min ∑ −h(λsi − λsj ), yisj i + ||yisj ||W i j . (8)
ij

λsi ,zs

ij

ys i< j

λsi

ij

ys i< j

The dual functions of the indicator functions for the convex sets W i j are 1-homogeneous,
of the form || · ||W i j := supw∈W i j wT ·. Recall that our label costs are not metric: ∀i < j < k :

6

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

Figure 3: Left: Without our non-metric extension, optimisation w.r.t. (7) can lower transition
costs by inserting another label (here 1 between 0 and 2). Right: A solution is to split the
gradients of the indicator functions and use direction-dependent variables xi j .
|yi j |W i j + |y jk |W jk ≥ |yik |W ik , does not hold. It was shown [53] that a regulariser of the form
(8) transforms any non-metric cost to the metric case. Figure 3 shows an example. Here,
an expensive transition between labels 0 and 2 will be replaced by two cheaper transitions
0–1 and 1–2. To prevent this, we replace the yi j with direction dependent variables xi j : We
rearrange (8) and combine the first summand with the regulariser from (7) to arrive at the
following equations (for now ignoring |s|):

∑ ∑h∇xi , λsi i − ∑hλsi , ∑ yisj [i < j] − ysji [i > j] i,
s

i

i

j6=i

with [·] denoting Iverson brackets. Let xi j := [yi j ]+ and x ji := [−y ji ]+ , where [·]+ := max(0, ·).
Expanding the gradient (5) we get, per simplex s,

∑ ∑ λsi xvi ([Jv ]+ − [−Jv ]+ ) − hλsi , ∑ (xi j − x ji )i,
i v∈s

(9)

j:i6= j

which we analyse further to achieve non-metric costs. It was observed in [53] that the xi j ∈
Rd can be interpreted as encoding the “label mass” that transitions from label i to label j
in a specific direction. Positivity constraints (by definition) on the xi j avoid the transport of
negative label mass. To anchor transport of mass on the actual mass of label i present at a
vertex, we introduce the variables xii for the mass that remains at label i, and split the above
constraints into two separate sets with the help of additional dual variables θ :
λsi ( ∑ xvi [Jv ]+ − ∑ xsi j ) + θsi ( ∑ xvi [−Jv ]+ − ∑ xsji ) + ∑ δ≥0 (xsi j ).
v∈s

v∈s

j

j

(10)

i, j

Note that this construction is only possible because our elements (simplices) are of strictly
positive volume, in contrast to zero sets in Ω w.r.t. the Lebesgue measure. Finally, we can
write down our discrete energy in the Lagrange basis defined on the simplex mesh M:
min max ∑ ∑ ρvi xvi + δ∆ (xvi ) + ∑ ∑ |s|||xsi j − xsji ||W i j +

xi ,xi j λ i ,θ i v∈V i

∑∑ ∑
s

i< j s∈S

∑ xsji ) + ∑ ∑ λsi (∑ xvi [Jv ]+ − ∑ xsi j ) + ∑ δ≥0 (xsi j ),

θsi ( xvi [−Jv ]+ −
v∈s
j
i

s

i

v∈s

j

(11)

i, j

where we have moved the weighting with |s| from the constraint set to the regulariser, and
denote by δ∆ (·) the indicator function of the unit simplex.

4

Semantic Reconstruction Model

A prime application scenario for our FEM multi-label energy model (11) is 3D semantic
reconstruction. In particular, we focus on an urban scenario and let our labeling functions

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

7

Figure 4: (a): Wulff-shape (red) with isolines. (b): Minkowski sum of two Wulff-shapes.
(c): Simplices are split after inserting new vertices (blue) close to the surface (green). Right:
Initialisation of vertices after refinement. (d): Finite differences on a regular grid ([18]) only
cover constraints in green areas, the P1 basis covers all of the domain Ω.
encode freespace (i = 1), building wall, roof, vegetation or ground. Objects that are not
explicitly modeled are collected in an extra clutter class. We define the data cost ρ at a 3Dpoint x ∈ Ω as in [18]: project x into the camera views c ∈ C, and retrieve the corresponding
depth dˆc (x) and class likelihoods σci (x) from the image space. The σ i are obtained from a
MultiBoost classifier. For the depth we look at the difference between the actual distance
dc (x) to the camera and the observed depth: d(x, c) := dc (x) − dˆc (x). For the freespace label
we always set the cost to 0, for i 6= 1 we define:
ρ i (x) := ∑ σci (x)[(k− 1)ε ≤ d(x, c) ≤ kε] + β [|d(x, c)| ≤ kε] sign(d(x, c)).

(12)

c∈C

This model assumes independence of the per-pixel observations, and exponentially distributed inlier noise in the depthmaps, bounded by a parameter kε (k=3 in practice). It is
essentially a continuous version of [18], see that paper for details. The parameter ε sets a
lower bound for the minimal height of the simplices in the tesselation, and thus defines the
target resolution. The discretisation of the data cost involves a convolution with the respective basis functions, which can be approximated via sampling. Please refer to the supplementary material for details. The Wulff-shapes W i j in (11) are given as the Minkowski sum of the
L2 -Ball, B2κ i j := {x ∈ Rd |kxk2 ≤ κ i j } and an anisotropic shape Ψi j : W i j := Ψi j ⊕ B2κ i j . In the
isotropic part, κ i j contains the neighbourhood statistics of the classes. The anisotropic part
Ψi j models the likelihood of a transition between classes i and j in a certain direction. Fig. 4
(a,b) shows an example. For our case we prefer flat, horizontal surfaces at the following label
transitions: ground-freespace, ground-building, building-roof, ground-vegetation and rooffreespace. A second prior prefers vertical boundaries for the transitions building-freespace
and building-vegetation. More details on the exact form can be found in [18].
The energy (11) is already in primal-dual form, such that we can apply the minimisation
scheme of [9], with pre-conditioning [35]. That numerical scheme requires us to project
onto shapes that are Minkowski sums of convex sets. In our case, the sets are simple and
the projection onto each shape can be performed in closed form. We employ a Dykstralike projection scheme [5], which avoids storing additional variables and proves remarkably
efficient, see supplementary material. We also project the labeling functions xi directly onto
the unit simplex [48]. In order to extract the transition surface, we employ a variant of
marching tetrahedra (triangles) [42], using the isolevel at 0.5 for each label.
We conclude with two interesting remarks. First, note that a tesselation with a regular
grid [18] can be seen as a simplified version of our discretisation in the P1 Lagrange basis. In
Fig. 4d we consider the 2D case of the regular grid used in [18]. Here, variables are defined
at voxel level. In its dual graph, the vertex set consists of the corners of the primal grid cubes,
leading to shifted indicator variables. Per vertex the data term is mainly influenced from the
cost in its Voronoi area. Similarly, [18] evaluates the data cost at grid centers, approximately

8

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES
overall
acc. [%]
Scene 1
Scene 2

Tetra

Octree

MB

84.0
92.5

83.9
92.8

82.5
89

Figure 5: Left: Synthetic 2D scene, colors indicate Table 1: Quantitative comparison
ground (gray), building (red) and roof (yellow). Mid- with octree model [4] and Multidle: Control mesh. Right: Example reconstruction. Boost input data.
corresponding to integration within the respective Voronoi-area of grid cell. Furthermore,
taking finite differences in this regular case corresponds to verifying constraints for only one
of the two triangles (Fig. 4d). The supplementary material includes a more formal analysis.
Second, our formulation is adaptive, in the sense of [4]: Hierarchical refinement of the
tesselation can only decrease our energy. Hence, our scheme is applicable when refining the
model on-the-fly. We again must defer a formal proof to the supplementary material and give
an intuitive, visual explanation in Fig. 4c . Assume that (x∗ , λ ∗ , θ ∗ ) is an optimal solution
for a certain triplet M = {F,V, S}. Then a refined tesselation M̂ = {F̂, V̂ , Ŝ} can be found by
introducing additional vertices, i.e. V ⊂ V̂ (ideally on the label transition surfaces). To define
a new set of simplices, we demand that no faces are flipped, ∀ŝ ∈ Ŝ, ∃s ∈ S with ŝ ∩ s = ŝ.
Then one can find a new variable set and data cost ρ̂ with the same energy: We initialise
the new variables from the continuous solution at the respective location, and find new ρv̂
by integration. Subsequent minimisation in the refined mesh can only decrease the energy.
The argument works in both ways: Vertices that have the same solution as their adjacent
neighbors can be removed without changing the energy. For now we stick to this simple
scheme, future work might explore more sophisticated ideas, e.g. along the lines of [17].

5

Evaluation

Before we present results on challenging real 3D data we evaluate our method in 2D on a
synthetic dataset. All results are obtained with a multi-core implementation, on a 12-core, 3.5
GHz machine with 64GB RAM. For clarity, we only present the Lagrange discretisation. We
refer to the supplementary material for an evaluation of the Raviart-Thomas discretisation.
Input Data. We create a synthetic 2D scene composed of 4 labels: free space, building, ground and roof, surrounded by 17 virtual cameras. To replace depth maps and classlikelihood images, we extract 2D points on the boundary “surface” and assign ground truth
label costs to each point. For the evaluation in 3D, we use three real-world aerial mapping
data sets. Our method requires two types of input data: depthmaps and pixel-wise class
probabilities (cf . Sec. 4). Moreover, we build a control mesh M around the initially predicted surface, to facilitate our FEM discretisation. Ideally, the control mesh enwraps the
true surface, using a finer meshing close to it. We densely evaluate the data cost at the vertices of a regular data cost grid and let each control vertex accumulate the cost of its nearest
neighbours in that grid, to approximate an integration over its Voronoi cell.
2D Lagrange results. Fig. 5 illustrates the result we obtain in a perfect setting. The original 2D image serves as ground truth for our quantitative evaluation. In this baseline setting,
our method achieves 99.8% of overall accuracy and 99.7% of average accuracy, confirming
the soundness of our Lagrange discretisation. In order to evaluate how our model behaves
in a more realistic setting, we conduct a series of experiments where we incrementally add
different types of perturbations. Our algorithm is tested against: (i) noise in the initial 2D
point cloud, respectively depth maps, (ii) wrong class probabilities and (iii) ambiguous class

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

9

Figure 6: (a) Illustration of the control mesh foundation. Dots represent values of the datacost grid and crosses the control vertices. Voronoi cells of the control vertices are depicted
with dashed grey lines and the control mesh with a solid black line. Colors indicate ground
(purple), building (red), roof (yellow), free space (cyan) and no datacost (black). (b) Quantitative evaluation of Lagrange FEM method w.r.t. different degradations of the input data.
probabilities of random subsets, (iv) missing data, e.g. deleting part of a facade to simulate
unobserved areas, (v) sparsity of the initial point cloud. Fig. 6b illustrates the influence of
defective inputs. Under reasonable assumptions on the magnitude of the investigated perturbations, we do not observe a significant loss in accuracy. The reconstruction quality starts to
decrease if more than half of the input data is misclassified or if the input point cloud is excessively sparse, meaning that >50% of the input is wrong or nonexistent. Average accuracy
is naturally more sensitive, due to the larger relative error in small classes.
Influence of the control mesh. Recall from (12) that the data cost of a control vertex
v ∈V is approximately equal to an integral of ρ over its respective Voronoi area (cf . Fig. 6a,
left). Therefore, but also because of the sign change in (12), vertices close to the surface
receive small cost values and are mainly steered by the regulariser, i.e. these vertices realise
a smooth surface. On the other hand, vertices that integrate only over areas with positive or
negative sign determine the inside/outside decision, but are more or less agnostic about the
exact location of the surface. We conclude that a sufficient amount of control vertices should
lie within the band [dˆ − 3ε; dˆ + 3ε] defined by the truncation of the cost function around
the observed depth dˆ (cf . (12) and [18]). Ideally control vertices are equally distributed
along each line-of-sight in front and behind the putative match (cf . Fig. 6a, middle column).
Undersampling within the near-band can lead to smooth, but inaccurate results (cf . Fig. 6a,
top right). Unobserved transitions, e.g. building-ground or roof-building, can also lead to
problems if the affected simplices are too large. To mitigate the effect, we add a few vertices
(e.g., a sparse regular grid) on top of the control mesh (cf . Fig. 6a, bottom row). Finally,
oversampling each line-of-sight in order to increase the resolution of the control mesh is not
recommended, the right spacing is determined by the noise level and ε and k, chosen in (12).
To conclude, it is an important advantage of the FEM framework that additional vertices
can be inserted as required, without changing the energy. In future work we will use this
flexibility to develop smarter control meshes, possibly as a function of the local noise level.
3D Lagrange results. To test our algorithm on real world data, we focus on a dataset
from the city of Enschede. Complementary results for other datasets are shown in the supplementary material. As baseline we use [4], the current state-of-the-art in large-scale semantic

10

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

3D reconstruction. Due to the lack of 3D ground truth, we follow their evaluation protocol
and back-project subsets of the 3D reconstruction to image space, where it is compared to a
manual segmentation. As can be seen in Fig. 7 and Tab. 1, the two results are similar in terms
of quantitative correctness. We note that measuring labeling accuracy in the 2D projection
does not consider the geometric quality of the reconstruction within regions of a single label.
Figs. 1 and 7 show city-modelling results obtained from (nadir and oblique) aerial images. Visually, our models are crisper and less “clay-like”. Compared to axis-aligned discretisation schemes, e.g. [4, 18], our method appears to better represent surfaces not aligned
with the coordinate axis, and exhibit reduced grid aliasing. Both effects are consistent with
the main strength of the FEM framework, to adapt the size and the orientation of the volume
elements to the data. Small tetrahedra, and vertices that coincide with accurate 3D points on
surface discontinuities, favour sharp surface details and crease edges (e.g., substructures on
roofs). Faces that follow the data points rather than arbitrary grid directions mitigate aliasing
on surfaces not aligned with the coordinate axes (e.g., building walls). The freedom of a
local control mesh unleashes the power of the regulariser in regions where the evidence is
weak or ambiguous (e.g., roads, weakly textured building parts).
As already mentioned, our FEM framework can be readily combined with on-the-fly
adaptive computation, as used in the baseline [4]. Compared to their voxel/octree model,
adaptive refinement is straight-forward, due to the flexibility of the FEM framework, which
allows for the introduction of arbitrary new vertices. As a preliminary proof-of-concept, we
have tested the naive refinement scheme described in Sec. 4. We execute three refinement
steps, where we repeatedly reconstruct the scene and subsequently refine simplices that contain surface transitions, while lowering ε by half. Compared to computing everything at the
final resolution, this already yields substantial savings of 89–97% in memory and 82–93% in
computation time, without any loss in accuracy. Targetting ε ≥ √13 (measured w.r.t. a bounding box of 256 units), the runtimes for the tested scenes are 1h04m–1h47m and memory
consumption is 573–764 MB, on a single machine.

Figure 7: Quantitative evaluation of Scene 1 from Enschede. Left: One of the input images.
Middle left: Semantic 3D model. Middle right: Back-projected labels overlayed on the
image. Right: Error map, misclassified pixels are marked in red.

6

Conclusion

We have proposed a novel framework for the discretisation of multi-label problems, and
have shown that, in the context of semantic 3D reconstruction, the increased flexibility of
our scheme allows one to better adapt the discretisation to the data at hand. Our basic idea
is generic and not limited to semantic 3D reconstruction or the specific class of regularisers. We would like to explore other applications where it may be useful to abandon grid
discretisations and move to a decomposition into simplices.
Acknowledgements: Audrey Richard was supported by SNF grant 200021_157101. Christoph Vogel
and Thomas Pock acknowledge support from the ERC starting grant 640156, ’HOMOVIS’.

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

11

Supplementary Material
This document provides supplementary information to support the main paper. It is structured as follows: Sec. A gives more information about the data and pre-processing used in
our experiments, not mentioned in the paper due to lack of space. We hope that the added
details will help readers to better appreciate the experimental results. In Sec. B we show
complementary results obtained with the proposed Lagrange FEM method on other datasets,
as well as the full large-scale reconstruction of the city of Enschede. Sec. C contains technical details and formal proofs that had to be omitted in the paper. Finally, Sec. D discusses
our formalism for the case of the Raviart-Thomas basis (instead of Lagrange P1), leading
to piecewise constant labels. We also show results in 2D and 3D and a comparison to those
obtained with the Lagrange basis.

A

Input Data

For our real-world experiments, we start from aerial images, cf . Fig. 8. To mitigate foreshortening and occlusion, images are acquired in a Maltese cross configuration, with four oblique
views in addition to the classical nadir view. We orient the images with VisualSFM [50],
create depth maps from neighbouring views with Semi-global Matching [6, 19], and predict
pixel-wise class-conditional probabilities with a MultiBoost classifier [3]. The classifier is
trained on a few hand-labeled images, using the same features as [4]: raw RGB-intensities in
a 5 × 5 window, and 19 geometry features (height, normal direction, anisotropy of structure
tensor, etc.) derived from the depth map.

Figure 8: Input data. Left: aerial input images for one position (four oblique views to the
north, east, south, west, and a nadir image). Middle: oriented image block. Right: depth
map and class probabilities (visualised by maximum-likelihood labels).

B

Additional Visualizations

We have tested our semantic reconstruction method on several (synthetic) 2D and (real)
3D datasets. Here we provide additional examples to give the reader an impression of the
variety of cases tested in our evaluation. We apply the same prior models as for our 3D
reconstructions. We prefer flat, horizontal structures in the model for the following label

12

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

Figure 9: Example scenes of our 2D data set and results obtained with our Lagrange FEM
method. Left: Data term at vertices of the control mesh. Colors for the data cost indicate:
free space/empty space (cyan), building (red), ground (pink), roof (yellow), occupied space
(green) and no data cost (black). Middle: Semantic 2D model. Right: Classification result,
misclassified pixels are depicted in red.

transitions: ground-freespace, ground-building, building-roof and roof-freespace. The second prior applies to the transition building-freespace and prefers vertical boundaries. Fig. 9
shows examples for different degradations of the synthetic input (many more cases were
tested). In the top row, we simulate imperfect classifier input by adding noise to the semantic class likelihoods. In our experience, the method is still able to reconstruct the geometry
quite well, but sometimes assigns the wrong label. A closer inspection reveals that, locally,
the roof and building classes are confused in locations where the class likelihoods are significantly wrong. The global geometry and labels in other regions remain unaffected. The
second row gives an example of missing input data, a frequent situation in the real world,
due to occlusions and constraints on camera placement. Fortunately, missing data does not
seem to greatly challenge our method. In fact, our method is specifically designed to work
well for these cases and complete the outline, relying on the prior assumptions about pairwise class transitions and class-specific local shape. In the last row we utilise only a sparse
control mesh, even near the surface. The method can still recover the geometry, but struggles
to determine the correct semantic labeling near the (unobserved) roof-to-building transition.
The adaptive version of our method is designed to avoid exactly that case. It refines the
control mesh near the predicted transitions, effectively increasing the resolution at the most

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

13

Figure 10: Additional datasets. First row: Original aerial images. Second row: Semantic 3D
models obtained with the Lagrange FEM method for Enschede (left), Zürich (middle) and
Dortmund (right). For the latter, light green denotes an additional class grass and agricultural fields.
promising locations.
Fig. 10 shows city models obtained from two additional aerial datasets (Zürich, Switzerland and Dortmund, Germany), and a further patch from Enschede. These results qualitatively illustrate that our method works for different image sets and architectural layouts.
Finally, we show the complete semantic 3D reconstruction of Enschede. Fig. 11 shows
the model rendered in an oblique view, together with the corresponding viewpoint in Google
Earth, to illustrate its accuracy and high level-of-detail.

C

Proofs

In this section we give the technical proofs promised in the main paper, as well as further
details about the optimisation. We start with a discussion of the extension to non-metric
energies, and its consequences on the equivalence of continuous and discrete models.

C.1

Non-metric priors: continuous vs. discrete

The main message of this section is that a non-metric model does not exist in the continuous
view, unless one imposes additional constraints on the function spaces. We briefly explain
why: Let’s look at the transition boundary between two labels i and k. Without additional
constraints, one can always introduce a zero set with label j between the two, i.e., a set
with Lebesgue measure 0 in the domain space. If the transition costs are not metric, then
the cost for the label pair {i, k} is potentially higher than the sum of the costs for {i, j}
and { j, k}. Inserting the zero set will avoid that extra cost and the energy will be underestimated. In other words, let S be a segmentation of Ω into regions Si and Sk , labeled with
i and k respectively. Assume further that their costs do not fulfill the triangle inequality
w.r.t. another label j. Then one can find a sequence of segmentations Sn := {Sin , Snj , Skn } with

14

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

Figure 11: Large-scale semantic 3D reconstruction of Enschede (Netherlands), computed
from aerial images with our Lagrange FEM method. Top: View from Google Earth (not
used during reconstruction). Bottom: Our model from matching viewpoint.

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

15

S = lim (Sn ): the label j disappears in the limit, such that lim inf E(Sn ) < E(S). Hence,
n→∞
n→∞
metric transition costs are a necessary condition for the lower semi-continuity of the energy
functional. Methods that try to resolve the issue with additional constraints on the function
space, for instance by demanding Lipshitz continuity of the labeling functions, are an active
research area, e.g. [7], but are beyond the scope of this work.
The above conceptual problem does have consequences for a practical implementation:
Any discretisation of the domain will ultimately consist only of a finite number of elements
of measurable (> 0) volume. Thus, the label j in the example will not disappear completely
from the solution, and the computed energy matches the solution. In practice, one can simply
prescribe a minimum edge length in the tesselation, since one cannot refine infinitely. Note
that this also constrains the Lipshitz constant of the labeling functions; they are restricted to
values between 0 and 1, such that the Lipshitz constant of functions f ∈ P1 (M) defined on the
mesh M = {V, F, S} is bounded by minv∈s,s∈S ||Jv ||, cf . (13). Because we utilise a Delaunay
triangulation/tetrahedralisation of the domain and also limit the minimal dihedral angle, a
further constraint on the edge length implies a bound on the Lipshitz constant. Note also, our
analysis implies that a discrete solution in the non-metric setting does not have a continuous
counterpart, and consequently investigations of the limiting case, i.e., convergence analysis
after infinite refinement of the tesselation, are futile.

C.2

Gradient in the Lagrange basis

We show that gradients of functions in the P1 (Lagrange) basis are constant per simplex s
and given by:
∇φs = ∑ φv Jv
(13)
v∈s

| fv |
Here, the coefficients φv ∈ R and Jv ∈ Rd denote a vector of length |s|d
, normal to the face
fv opposite to vertex v, and pointing inwards towards the center of the simplex. Recall that
| fv | is the area of face fv and |s| is the volume of simplex s.
The gradient can be obtained with basic algebra. First, notice that the gradient of φs in
(13) has to fulfill hvl − vk , ∇φs i = φvl − φvk , meaning that integration along the edge leads to
the respective change in φs . After collecting a sufficient number of linear equations of this
form, one can directly solve the resulting linear system. Since Jv is, by definition, orthogonal
to all edges that do not involve vertex v, we arrive at (13).
Formally, we pick one vertex v of simplex s and compile for l = 1 . . . d (vl 6= v) equations
of the form hvl − v, ∇φs i = φvl − φv . By construction, hJvl , vk − vi = δk=l . The vector Jvl is
normal to face fvl . The scalar product of the edge (vl − v) and the normal is the ”height”
within the simplex, so with the chosen scaling of Jvl we have hJvl , vl − vk i = 1 for any k 6= l.
Thus multiplying each side of our equation system by a matrix with the vectors Jvl , l =
1 . . . d as columns leads to: ∇φs = ∑l Jvl (φvl − φv ). If we can show that ∑l Jvl = −Jv , then we
arrive at the desired expression (13). For vk , v j 6= v, h∑l Jvl , v j − vk i = hJvl , v j i − hJvk , vk i = 0
and h∑l Jvl , vk − vi = hJvk , vk − vi = 1. All equations are also fulfilled by Jv in place of ∑l Jvl ,
which concludes the proof.

C.3

Data Term for Lagrange basis

We again start from the ideas in the main paper. We have to convert the continuous data
costs ρ i into discrete form (in a practical implementation, “continuous” means that the cost

16

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

can be evaluated at any z ∈ Ω). In our basis representation, we can get discrete cost values
for the basis elements by convolving the continuous cost with the respective basis function.
For simplicity, we consider the P1 basis function here. Thus, we seek a cost per vertex ρvi .
In detail we obtain:
Z
Ω

xi (z)ρ i (z)dz = ∑

Z

s

s

xsi (z)ρ i (z)dz = ∑

Z

Z

∑ φv

∑

s∈N(v) s

v∈V

|

∑ φv p1s,v (z)ρ i (z)dz =

s v∈s

s

!
p1s,v (z)ρ i (z)dz
{z

= hρvi , xvi i.

(14)

}

:=ρvi

To numerically compute ρvi , we sample ρ i at a finite number of locations z ∈ Ω. For
each z we determine into which simplex s it falls, and accumulate the contributions of ρ i (z)
over all i = 1 . . . m, weighted by their barycentric coordinates. The final step is to scale ρvi
by ∑s∈N(v) |s|
d and divide by the sum of weights assigned to vertex v. In other words, we
compute the sample mean and scale it by the area covered by the vertex. In our current
implementation ρ is sampled at regular grid points, without importance sampling. This
simple strategy is indeed very similar to the method employed in [4, 18]. There, the data
cost is evaluated on a regular grid, by reprojecting grid vertices into each image, computing
the data term, and adding its respective contribution to the grid location. Such a “per-voxel
accumulation” is equivalent to integrating the data cost within the respective Voronoi-area
of a vertex in the dual grid: the latter is proportional to the number of regular samples that
fall into a Voronoi-cell and therefore have the respective vertex as nearest neighbour. Hence,
summing the individual contributions directly corresponds to integrating the data term within
the Voronoi region.

C.4

Grid vs. P1

Here, we detail why the grid-based version with finite differences (corresponding to [18])
can be seen as an approximation of our proposed FEM discretisation with P1 basis elements,
if the vertices (cells) are aligned in a regular grid. Without loss of generality we consider a
grid of edge length 1, and note that in this case the gradient for a function f : Ω ⊂ Rd → R
at a grid point x, evaluated with forward differences becomes:
∇ f = ( fx+e1 − fx , . . . , fx+ed − fx )T = ∑ ei fx+ei − ∑ ei fx = ∑ fx+ei Jx+ei + fx Jx ,
i

(15)

i

with ei the unit vector in direction i. We have used the identity ei = Jx+ei , according to the
definition in Sec. C.2, and obtain the last equality from ∑l Jvl = −Jv , cf . Sec. C.2. This is exactly the formula for the gradient of the corresponding P1 function in the simplex defined by
the vertices {x, x + ei }di=1 . Accordingly, if implemented as finite differences, the constraints
on the dual vector field λ , see Eq. (2) from the paper, are only checked within the respective
simplex, but not in the whole domain (1/2 of the domain in 2D; 1/6 in 3D). Note also that,
with grid-aligned vertices, the simplex in question cannot be part of a partition of Ω ⊂ Rd ,
unless d ≤ 2: edges of adjacent faces would intersect.
Fig. 12 illustrates the specific case with d = 2. On the left, the regular grid (green) and
the triangle (simplex, red). Grid centers correspond to vertices in the (triangle-)mesh. The
grid corresponds to the discretisation used in [18], whereas the simplex mesh is used in this

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

17

Figure 12: Left: Grid (green) and simplex mesh (red) cover the same domain, but are offset
against each other. Grid centers correspond to vertex positions. The gradient in a triangle
(middle) corresponds to the gradient computed with forward differences (right).
paper. The gradient of the lower left triangle for the simplex mesh corresponds exactly to
the one computed via forward differences as shown in (15). Consequently, discretisation
via finite differences is a special case of our method, where the elements are layed out on a
regular grid, and the constraints are tested only in the upper right triangle (in the 2D case).

C.5

Adaptiveness

We have stated in the paper that our formulation is adaptive, in the sense that a hierarchical
refinement of the tesselation can only decrease the energy. We have also explained a way to
find the refined tesselation of Ω, by introducing additional vertices and splitting simplices s,
such that no faces f ∈ F are flipped; and we have put forward a procedure to initialise the
new variables. Here, we formally prove that the described scheme is sound.
Let (x∗ , λ ∗ , θ ∗ ) be the solution for a triplet M = {F,V, S}. And let M̂ = {F̂, V̂ , Ŝ} be the
refined mesh with V ⊂ V̂ and ∀ŝ ∈ Ŝ, ∃s ∈ S with ŝ ∩ s = ŝ. Furthermore, we define the sets
V̄ = V̂ \V and S̄ = Ŝ \ S to denote newly introduced vertices and simplices. Our construction
works by induction, i.e. we introduce one vertex v at a time. The vertex is assumed to lie in
d+1
simplex s ∈ S, s = {vk }d+1
/ ∀k 6= l. By
k=1 , which is split into simplices {s̄k }k=1 with s̄k ∩ s̄l = 0,
|s̄k |
d+1 |s̄k |
definition, vertex v̄ has the barycentric coordinates |s| , i.e. v̄ = ∑k=1 |s| vk .
We initialize the labeling variables at new vertices v̄ ∈ V̄ via barycentric interpolation:
s̄k | i
xvk . Dual variables of the new simplices s̄k , and also the transition variables
xv̄i = ∑vk ∈s ||s|
xi j , are simply copied from their enclosing simplex s: λs̄ik := λsi , θs̄ik := θsi and xs̄i kj := xsi j .
The data terms at the new vertex v̄, as well as at the other vertices vk of simplex s, are
(re)computed following (14), see also Fig. 13. We call the new variables x̄, λ̄ , θ̄ and claim
that EM (x∗ , λ ∗ , θ ∗ ) = EM̂ (x̄, λ̄ , θ̄ ) and that x̄, λ̄ , θ̄ is feasible. The latter is trivially the case
(cf . Eq. (11) from the paper); transition variables remain positive by construction, the newly
introduced labeling variables still fulfill the simplex constraints, and they induce the same
gradients in the new simplices as in the parent simplex. The same applies for the cost of the
regulariser in the new simplices: ∑i< j |s|||xsi j − xsji ||W i j = ∑i< j ∑s̄k ∈s |s̄k |||xs̄i kj − xs̄jik ||W i j . More
R
interesting is the data cost. We introduce the notation ρv,s := s p1s,v (z)ρ i (z)dz for the data
cost at vertex v, originating from the integral over simplex s. By induction, we need to only

18

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

Figure 13: Updated data term after adding a new vertex.
verify the following equality for simplex s, which is split into {s̄k }d+1
k=1 :

∑ xvi k ρvk ,s = ∑ xvi k ∑ ρvk ,s̄ j + ∑ xv̄i ρv̄,s̄ j = ∑ xvi k ∑ ρvk ,s̄ j +
k

k

j

j6=k

k

j6=k

|s̄k |
xvi k ρv,s̄ j
|s| ∑
j

(16)

Recall we have ”old” vertices vk , s = {vk }d+1
k=1 and a new vertex v̄. According to (16), we
must verify:
|s̄k |
|s̄k |
∑ ρv,s̄ j ⇔ s ρ(z)p1s,vk (z)dz = s ρ(z) ∑ p1s¯j ,vk (z)+ |s| ∑ p1s̄ j ,v (z)dz.
|s|
j
j
j6=k
j6=k
(17)
It is sufficient to show
Z

Z

ρvk ,s = ∑ ρvk ,s̄ j +

p1s,vk (z) =

∑ p1s¯j ,vk (z) +

j6=k

|s̄k |
p1s̄ j ,v (z) ∀z ∈ s.
|s| ∑
j

(18)

The right hand side represents a linear function for each {s̄k }d+1
k=1 . We can check if both sides
agree on d + 1 points in each simplex, which is easy to verify. The locations we check –
substitute z on both sides of Eq. (18) – are {vk }d+1
k=1 and v. These are the defining vertices of
the d + 1 simplices {s̄k }d+1
.
Left
and
right
hand
side vanish, except for v and vk . Finally, we
k=1
|s̄k |
get |s| for v and 1 for vk on both sides.
In our adaptive version, we directly follow the proof and split simplices with the introduction of a single new vertex. We emphasise again that this splitting schedule is merely a
proof of concept. The FEM discretisation allows for more sophisticated refinement schemes,
e.g., along the lines of [17], or flipping edges according to the energy functional, etc.

C.6

Optimisation

The energy (11) from the paper is given in primal-dual form, optimisation with existing tools
is straight-forward. We apply the minimisation scheme of [9], with pre-conditioning [35].
Internally, that algorithm however requires the projection onto the Wulff-shapes W i j , which
is slightly more involved.
C.6.1

Proxmap for the Minkowski sum of convex sets

Recall that, per label pair {i, j}, our Wulff-shapes are of the form W i j := Ψi j ⊕ B2κ i j . They
are the Minkowski sum of two simple convex sets. Recall that the Ψi j encode the direction
dependent likelihood of a certain label transition. In our case, all Wulff-shapes permit a

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

19

closed form projection scheme, such that we solve the following sub-problem as proximal
step, independently per simplex s:
1
1
arg min ||xi j − xi j ||2 + ||x ji − x ji ||2 + sup wT (xi j − x ji ) + ι≥0 (xi j ) + ι≥0 (x ji ).
2
xi j ,x ji 2
w∈Ψi j ⊕B2

(19)

κi j

For the following derivation we rename the two sets W1 := Ψi j and W2 := B2κ i j . In order to
decouple the argument within the regulariser, we introduce auxiliary variables {yk , zk }2k=0
and additional Lagrange multipliers {µk , λk }2k=0 , and replace xi j and x ji respectively:
min

xi j ,x ji ,yk ,zk

1
1
max ||xi j − xi j ||2 + ||x ji − x ji ||2 +
2
µk ,λk 2
(20)

2
T

ij

T

T

ji

sup w (yk − zk ) + ι≥0 (y0 ) + ι≥0 (z0 ) − ∑ λk (x − yk ) − µk (x − zk ).

∑

k∈{1,2}w∈Wk

k=0

Optimality w.r.t. xi j , x ji implies:
2

2

xi j = xi j + ∑ λk and x ji = x ji + ∑ µk ,
k=0

(21)

k=0

which, after reinserting into (20), leads to:
−1 2
−1 2
|| ∑ λk − xi j ||2 +
|| ∑ µk − x ji ||2 +
yk ,zk µk ,λk 2
2 k=0
k=0

min max

(22)

2
T

sup
w1 ∈W1 ,w2 ∈W2

T

T

T

w1 (y1 − z1 ) + w2 (y2 − z2 ) + ι≥0 (y0 ) + ι≥0 (z0 ) + ∑ λk yk + µk zk .
k=0

Applying Fenchel-duality yields:
max min
µk ,λk zk

−1 2
−1 2
|| ∑ λk − xi j ||2 +
|| ∑ µk − x ji ||2
2 k=0
2 k=0

(23)

2
T

− ιW1 (−λ1 ) − ιW2 (−λ2 ) − ι≤0 (−λ0 ) − ι≤0 (−µ0 ) + ∑ (λk + µk ) zk .
k=1

The latter summand requires λ1 = −µ1 and λ2 = −µ2 :
1 2
1 2
min || ∑ λk −xi j ||2+ || ∑ λk +x ji−µ0 ||2+ιW1 (−λ1 )+ιW2 (−λ2 )+ι≥0 (λ0 )+ι≥0 (µ0 ). (24)
2 k=1
µ0 ,λk 2 k=0
In this last form, we can apply a few iterations of block coordinate descent on the dual
variables and recover the update for xi j , x ji from (21).

D
D.1

Raviart-Thomas basis
Methodology

In this section, we show how to discretise the convex relaxation, Eq. (2) from the paper, for
the case of the Raviart-Thomas basis. For convenience, we restate the energy:
Z

m

i Ω

i=1

min max ∑ ρ i(z)xi(z)+hxi(z), ∇·λ i(z)idz, s.t. λ i(z)−λ j(z) ∈ W i j, ∑ xi(z)=1, xi(z)≥0. (25)
xi

λi

20

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

The Raviart-Thomas basis is chosen as a strong contrast to the (preferred) Lagrange
basis. With Raviart-Thomas functions, we model the dual functions λ in (25), within our trial
space. The Raviart-Thomas RTk (M) basis functions describe a div-conforming polynomial
basis of order k+1, i.e. the divergence of the modeled vector field is continuous across
simplices. We again discretise on a simplex mesh M = {F,V, S} with vertices v ∈ V, v ∈ Rd ;
faces f ∈ F defined by d vertices; and simplices s ∈ S defined by d + 1 vertices, which
partition Ω: ∪k sk = Ω, sl ∩ sk = fk,l ∈ F.
RT 0 (M) := {p : Ω → Rd |φ (x) := ∑ φs (x) with φs (x) := cs x + ds , cs ∈ R, ds ∈ Rd ,
s∈S

if x ∈ s and 0 else, and φs (x) is continuous for x ∈ fv (s) in direction ν sfv }.

(26)

Here, we have used ν sfv to denote the (outward-pointing) normal of face fv of simplex s. By
convention the face fv is located opposite the vertex v. We construct our linear basis with
functions that are defined for each face fv in a simplex s, and can be described in a local
form as:
| fv |
0
φs,v
(x) := (x − v)
if x ∈ s and 0 else,
|s|d
where we again let | fv | denote the area of the face and |s| the volume of the simplex. Let ν sfv
be the normal of face fv in simplex s, then the basis functions fulfill:
0
hφs,u
(x), ν sfv i := [u = v] ∀x ∈ fv ,

(27)

with [·] denoting the Iverson bracket.
These basis functions make up the global function space by enforcing a consistent orientation. For each face f we can distinguish its two adjacent simplices s+ and s− , by analysing
the scalar product of the vector 1 and the normal ν s±
f of the shared face f (by convention
±

±

again pointing outwards of the respective simplex). W.l.o.g., we define sisf := signhν sf ; 1i,
+

−

i.e. s+f ν sf = s−f ν sf . The global basis functions per face fv are then given by:
0
φs,v

:=

(
| fv |
sisfv (x − v) |s|d
0

if x ∈ s
else.

(28)

In each simplex, our vector-field φs (x) ∈ Rd can then be defined in the following manner,
with coefficients φ fv ∈ R:
0
φs (x) := ∑ φ fv φs,v
(x).
v∈s

By construction, cf . (27),(28), the vector-field is continuous along a face x ∈ f in direction
of the face normal ν f (of arbitrary, but fixed orientation), i.e. for neighbouring faces s+ and
s− we have:
±
|f|
hφs± (x), ν f i = sisf hx − v± , ν f i ± φ fv± = φ fv± .
(29)
|s |d
Here, v+ is the vertex in simplex s+ opposite to the shared face, and v− is the vertex in
simplex s− . Thus, our function in φ is a RT function iff for all neighbouring faces s± we
have φ fv+ = φ fv− . In other words, basis coefficients only exist for faces of the simplices.

21

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

The variables we are interested in are the labeling functions xi , which are members of
our test function space, composed of piecewise constant functions per simplex:
U 0 (M) := {u : Ω → R|u(x) := ∑ us (x), with us (x) = us if x ∈ s and 0 else}.

(30)

s∈S

Before we can utilise our new basis to discretise (25) we need a way to enforce the
constraints on our dual variables λ i (z)−λ j (z) ∈ W i j for all z ∈ Ω. It is sufficient to enforce
the constraints on the dual functions λ in (25) only at the face midpoints z fv := 1/d ∑w∈ fv w
of faces fv ∈ s. This ensures the constraints are also valid for any point in the simplex
s. Because the Wulff shapes are convex, it is sufficient to prove that a vector field φ (x) ∈
RT 0 (M), φ (x) := ∑s∈S φs (x) at any point x ∈ s can be written as a convex combination of the
values at the face midpoints:
φs (x) =

∑ αz fv φ (z fv ), with ∑ αz fv = 1.

fv ∈s

fv ∈s

After some elementary algebra it turns out that, if x = αi vi , then αz fv := (1 − d · αi ) encode
this convex combination. Furthermore, the value of φs at a location x ∈ s can be found by
linear combination of basis coefficients at the vertices of s:
φs (x) = ∑ sisfv (x − v)
v∈s

D.2

| fv |
φv .
|s|d

(31)

Discretisation

With these relations, we can discretise the energy (25) for labeling functions xi ∈ U 0 (M) and
dual vector-field λ i ∈ RT 0 (M). First, we convert the continuous data costs ρ i into a cost
per simplex ρsi , which can again be achieved
by convolving
the cost with the respective (per
R
R
simplex constant) basis function: ρsi := s us (z)ρ i (z)dz = s ρ i (z)dz. In practice, the integral
is computed via sampling. Next, we discretise the second part of our energy with the help of
the divergence theorem and (29):
Z

xi (z)∇·λ i (z)dx = ∑

Z

s∈S S

Ω

xsi ∇·λ (z)dx= ∑

Z

xsi hλ (z), ν(z)idz=

s∈S ∂ S

∑ xsi λ fiv | fv |sisfv

(32)

v∈s, s∈S

As shown, we need to verify the constraints only at face midpoints z fv . The vectors λs (z fv )
are linear in the basis coefficients for any z ∈ Ω, and the discretised version of (25) becomes
min max ∑ ρsi xsi +∑ xsi λ fiv | fv |sisfv , s.t. λsi (z fv )−λsj (z fv ) ∈ W i j , xsi ∈ ∆ ∀i < j, v ∈ s, s ∈ S.
xi

λ i s∈S

v∈s,s∈S

(33)
Here, we let ∆ encode the unit simplex. Finally, for every simplex s, we replace the constraint
set ∑i< j λsi (z fv )−λsj (z fv ) ∈ W i j in the same manner as for the Lagrangian basis. We introduce
auxiliary variables and Lagrange multipliers yis,j fv , ∀i < j, and exploit Fenchel-Duality to
obtain
max min ∑ ∑ ||yis,j fv ||W i j − ∑ ∑hλsi (z fv ),
λsi

ij

ys v∈s i< j

max min ∑ ∑
λsi

ij

ys v∈s i< j

v∈s i

||yis,j fv ||W i j −

∑∑

v∈s i

λ fiv

∑ yis,j fv − ∑ ys,ji fv i =

j:i< j

| fv |sisfv
|s|d

j: j<i

"

(34)

!#!
T

∑ (z f¯ −v)
f¯∈s

∑

j:i< j

yis,j f¯ −

∑

j: j<i

ys,ji f¯

.

22

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

Figure 14: Left: Synthetic 2D scene. Colors indicate ground (gray), building (red) and roof
(yellow). Middle: Zoom of the control mesh. Right: Reconstructed semantic 2D model.
Furthermore, recall that we use Neumann conditions at the boundary of Ω, which translates
into coeffients λ fi = 0, ∀ f ∈ ∂ Ω. Combining (33) and (34), we get the (metric) energy for
the Raviart-Thomas discretisation:
min max ∑ ∑ ρsi xsi + ||yis,j fv ||W i j + ι∆ (xsi )

xi ,yi j λ i s∈S i

+∑ ∑
v∈s i

λ fiv | fv |sisfv

xsi −

1
|s|d

"

!#!
T

∑ (z f¯ −v)
f¯∈s

∑

j:i< j

yis,j f¯ −

∑

j: j<i

ys,ji f¯

(35)

To extend it to non-metric pairwise costs, as in the Lagrangian case, we need to impose
additional assumptions. One possibility is to utilise basis functions for the dual variables,
which are continuous in all directions at the faces. In that case, it is only necessary to check
the constraints at the faces and not for each face in each simplex, i.e. the variables for yis+j , f
and yis−j , f merge into one set. Another possibility is to only force the normal component
along the faces of λ to be contained in the Wulff-shapes. In this direction, RT is already
continuous and the Lagrange multipliers y can be merged. This line of attack leads to a
scheme that is remarkably similar to belief propagation on a Markov random field, in the
sense that the discretisation lacks a continuous counterpart to begin with, and may lead to
stronger grid artifacts. We stop at this point and leave an investigation of such models to
future work.

D.3

2D results

Fig. 14 illustrates the result we obtain with the Raviart-Thomas FEM method (RT). We use
the same (perfect) baseline setting as for the Lagrange FEM method (P1) in the main paper.
In that setting, the RT method achieves 97.5% of overall accuracy and 92.8% of average
accuracy. While these results confirm that also the RT method is sound, they also show its
limitations compared to the Lagrange basis. Simplices not aligned with object boundaries,
straddling multiple labels, will necessarily introduce errors in the reconstruction. Note that
we do not used edge information to guide the meshing; especially since such information
is not available for our target application, semantic 3D reconstruction. We refer to the 3D
qualitative comparison (cf . Sec. D.4) for a more detailed analysis of the differences between
the two methods.

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

23

Figure 15: Quantitative evaluation of Raviart-Thomas FEM method w.r.t. different degradations of the input data.

Figure 16: Quantitative evaluation of Scene 1 from Enschede. Left: One of the input images.
Middle left: Semantic 3D model obtained with our Raviart-Thomas FEM method. Middle
Right: Back-projected labels overlayed on the image. Right: Error map, misclassified pixels
are marked in red.
NB: Further to this manuscript, the supplementary material contains a short video, which
shows the diffusion of the indicator function over 1000 iterations for both proposed methods.
We perform also the same series of experiments where we incrementally add different
types of perturbations, cf . Sec. 5 of the main paper. Fig. 15 shows the corresponding behaviour of our RT method. Generally speaking, both models shows a similar sensitivity to
defective inputs, but with a small edge for the Lagrange method, which consistently reaches
higher overall accuracy.

D.4

3D results

Fig. 16 shows a quantitative evaluation of the Raviart-Thomas basis, equal to the one of the
Lagrange basis presented in the main paper. As before, the colors encode building (red),
ground (gray), vegetation (green), roof (yellow) and clutter (blue). We summarise the outcomes in Tab. 2. The differences between the Lagrange basis and octree are vanishingly
small, on the other hand we notice a bigger gap between the Raviart-Thomas basis and octree. We also present a qualitative comparison of the two bases in Fig. 17. The differences
are immediately apparent, which confirms the numbers given in Tab. 2. Raviart-Thomas
labels entire simplices, so the reconstruction consists of piecewise constant elements. On
the contrary, the Lagrangian basis has the advantage that the labeling functions are linear
and can be interpreted as (signed) distance functions, such that a smooth iso-surface can be
extracted, here done with marching tetrahedra. Despite the piecewise constant reconstruction, the RT basis measures metric quantities – in contrast to, for instance, Markov random
fields, where pairwise distances between the simplices would have to be designed explicitly
to achieve similar effects.

24

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

Data set
Scene 1

Error measure
Overall acc. [%]
Average acc. [%]

Tetra P1
84.0
81.1

Tetra RT
81.9
79.1

Octree
83.9
80.6

MB
82.5
81.4

Table 2: Quantitative comparison of our two proposed FEM methods with octree model [4]
and MultiBoost input data [3].

Figure 17: Reconstruction with Raviart-Thomas (left) and with the Lagrange basis (right).
We deliberately select a low resolution, choose flat shading and plot mesh edges, to accentuate the differences. Please refer to the text for details.

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

25

References
[1] Yingze Bao, Manmohan Chandraker, Yuanqing Lin, and Silvio Savarese. Dense object
reconstruction using semantic priors. In CVPR, 2013.
[2] Sören Bartels. Total variation minimization with finite elements: Convergence and
iterative solution. SIAM, 50(3), 2012.
[3] D. Benbouzid, R. Busa-Fekete, N. Casagrande, F-D. Collin, and B. Kégl. MULTIBOOST: a multi-purpose boosting package. JMLR, 2012.
[4] Maros Blaha, Christoph Vogel, Audrey Richard, Jan Dirk Wegner, Thomas Pock,
and Konrad Schindler. Large-scale semantic 3d reconstruction: An adaptive multiresolution model for multi-class volumetric labeling. In CVPR, 2016.
[5] J. P. Boyle and R. L. Dykstra. A method for finding projections onto the intersection of
convex sets in Hilbert spaces. Lecture Notes in Statistics, 1986.
[6] G. Bradski. The OpenCV Library. Dr. Dobb’s Journal of Software Tools, 2000.
[7] Elie Bretin and Simon Masnou. A new phase field model for inhomogeneous minimal
partitions, and applications to droplets dynamics, 2017.
[8] Franco Brezzi and Michel Fortin. Mixed and Hybrid Finite Element Methods. SpringerVerlag New York, 1991.
[9] Antonin Chambolle and Thomas Pock. A first-order primal-dual algorithm for convex
problems with applications to imaging. JMIV, 40(1), 2011.
[10] Antonin Chambolle, Daniel Cremers, and Thomas Pock. A convex approach to minimal partitions. SIAM, 5(4), 2012.
[11] R. Courant. Variational methods for the solution of problems of equilibrium and vibrations. Bull. Amer. Math. Soc., 49(1), 01 1943.
[12] D. Cremers and K. Kolev. Multiview stereo and silhouette consistency via convex
functionals over convex domains. PAMI, 33(6), 2011.
[13] D. Cremers, T. Pock, K. Kolev, and A. Chambolle. Convex Relaxation Techniques for
Segmentation, Stereo and Multiview Reconstruction. In Advances in Markov Random
Fields for Vision and Image Processing. MIT Press, 2011.
[14] Brian Curless and Marc Levoy. A volumetric method for building complex models
from range images. SIGGRAPH, 1996.
[15] Ricardo G. Durán. Mixed Finite Element Methods. Springer Berlin Heidelberg, 2008.
[16] Yasutaka Furukawa and Jean Ponce. Accurate, dense, and robust multi-view stereopsis.
PAMI, 2010.
[17] Eitan Grinspun, Petr Krysl, and Peter Schröder. CHARMS: A Simple Framework for
Adaptive Simulation. SIGGRAPH, 2002.

26

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

[18] Christian Häne, Christopher Zach, Andrea Cohen, Roland Angst, and Marc Pollefeys.
Joint 3d scene reconstruction and class segmentation. In CVPR, 2013.
[19] Heiko Hirschmüller. Stereo processing by semiglobal matching and mutual information. PAMI, 2008.
[20] M. Jancosek and T. Pajdla. Multi-view reconstruction preserving weakly-supported
surfaces. In CVPR, 2011.
[21] Michael Kazhdan, Matthew Bolitho, and Hugues Hoppe. Poisson surface reconstruction. In EUROGRAPHICS, 2006.
[22] K. Kolev, T. Brox, and D. Cremers. Fast joint estimation of silhouettes and dense 3D
geometry from multiple images. PAMI, 2012.
[23] Ilya Kostrikov, Esther Horbert, and Bastian Leibe. Probabilistic labeling cost for highaccuracy multi-view reconstruction. In CVPR, 2014.
[24] Abhijit Kundu, Yin Li, Frank Dellaert, Fuxin Li, and James M. Rehg. Joint semantic
segmentation and 3d reconstruction from monocular video. In ECCV, 2014.
[25] Patrick Labatut, Jean-Philippe Pons, and Renaud Keriven. Efficient Multi-View Reconstruction of Large-Scale Scenes using Interest Points, Delaunay Triangulation and
Graph Cuts. In ICCV, 2007.
[26] L’ubor Ladický, Paul Sturgess, Christopher Russell, Sunando Sengupta, Yalin Bastanlar, William Clocksin, and Philip Torr. Joint optimisation for object class segmentation
and dense stereo reconstruction. In BMVC, 2010.
[27] Mats G. Larson and Fredrik Bengzon. The Finite Element Method: Theory, Implementation, and Applications. Springer Publishing Company, Incorporated, 2013.
[28] Jan Lellmann and Christoph Schnörr. Continuous multiclass labeling approaches and
algorithms. SIIMS, 4(4), 2011.
[29] Shubao Liu and David B. Cooper. Ray Markov random fields for image-based 3d
modeling: Model and efficient inference. In CVPR, 2010.
[30] William E. Lorensen and Harvey E. Cline. Marching cubes: A high resolution 3d
surface construction algorithm. In SIGGRAPH, 1987.
[31] F. Maggi. Sets of Finite Perimeter and Geometric Variational Problems: An Introduction to Geometric Measure Theory. Cambridge Studies in Advanced Mathematics.
Cambridge University Press, 2012.
[32] Richard A. Newcombe, Shahram Izadi, Otmar Hilliges, David Molyneaux, David
Kim, Andrew J. Davison, Pushmeet Kohli, Jamie Shotton, Steve Hodges, and Andrew
Fitzgibbon. Kinectfusion: Real-time dense surface mapping and tracking. In ISMAR,
Washington, DC, USA, 2011.
[33] Matthias Nießner, Michael Zollhöfer, Shahram Izadi, and Marc Stamminger. Real-time
3d reconstruction at scale using voxel hashing. ACM Trans. Graph., 32(6), November
2013.

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

27

[34] Claudia Nieuwenhuis, Eno Töppe, and Daniel Cremers. A survey and comparison of
discrete and continuous multi-label optimization approaches for the potts model. IJCV,
104(3), 2013.
[35] Thomas Pock and Antonin Chambolle. Diagonal preconditioning for first order primaldual algorithms in convex optimization. In ICCV, 2011.
[36] Thomas Pock, Daniel Cremers, Horst Bischof, and Antonin Chambolle. Global solutions of variational models with convex regularization. SIIMS, 3(4), 2010.
[37] J. Reddy. An Introduction to the Finite Element Method. McGraw-Hill Education,
2005.
[38] Walter Ritz. Über eine neue methode zur lösung gewisser variationsprobleme der mathematischen physik. Journal für die reine und angewandte Mathematik, 1909.
[39] Nikolay Savinov, L’ubor Ladický, Christian Häne, and Marc Pollefeys. Discrete optimization of ray potentials for semantic 3d reconstruction. In CVPR, 2015.
[40] F. Steinbruecker, J. Sturm, and D. Cremers. Volumetric 3d mapping in real-time on a
cpu. In ICRA, 2014.
[41] E. Strekalovskiy and D. Cremers. Generalized ordering constraints for multilabel optimization. In ICCV, 2011.
[42] G. Treece. Regularised marching tetrahedra: improved iso-surface extraction. Computers & Graphics, 23(4), August 1999. ISSN 00978493.
[43] Ali Osman Ulusoy, Michael J. Black, and Andreas Geiger. Patches, planes and probabilities: A non-local prior for volumetric 3D reconstruction. In CVPR, 2016.
[44] Ali Osman Ulusoy, Michael J. Black, and Andreas Geiger. Semantic multi-view stereo:
Jointly estimating objects and voxels. In CVPR, 2017.
[45] Vibhav Vineet, Ondrej Miksik, Morten Lidegaard, Matthias Nießner, Stuart Golodetz,
Victor A. Prisacariu, Olaf Kähler, David W. Murray, Shahram Izadi, Patrick Perez, and
Philip H. S. Torr. Incremental dense semantic stereo fusion for large-scale semantic
scene reconstruction. In ICRA, 2015.
[46] George Vogiatzis, Carlos Hernández Esteban, Philip H. S. Torr, and Roberto Cipolla.
Multiview stereo via volumetric graph-cuts and occlusion robust photo-consistency.
PAMI, 29(12), 2007.
[47] H. H. Vu, P. Labatut, J. P. Pons, and R. Keriven. High accuracy and visibility-consistent
dense multiview stereo. PAMI, 34(5), May 2012.
[48] Weiran Wang and Miguel Á. Carreira-Perpiñán. Projection onto the probability
simplex: An efficient algorithm with a simple proof, and an application. CoRR,
abs/1309.1541, 2013.
[49] Martin Weber, Andrew Blake, and Roberto Cipolla. Sparse finite elements for geodesic
contours with level-sets. In Tomás Pajdla and Jiří Matas, editors, ECCV, 2004.

28

RICHARD, VOGEL, ET AL.: SEMANTIC 3D RECONSTRUCTION WITH FEM BASES

[50] Changchang Wu. VisualSFM: A visual structure from motion system, 2011.
[51] C. Zach, T. Pock, and H. Bischof. A globally optimal algorithm for robust TV-L1 range
image integration. In ICCV, 2007.
[52] Christopher Zach. Fast and high quality fusion of depth maps. 3DV, 2008.
[53] Christopher Zach, Christian Häne, and Marc Pollefeys. What is optimized in convex
relaxations for multilabel problems: Connecting discrete and continuously inspired
MAP inference. PAMI, 2014.

