# FaSS-MVS - Fast Multi-View Stereo with Surface-Aware Semi-Global Matching
from UAV-borne Monocular Imagery

Boitumelo Ruf<EMAIL_ADDRESS><EMAIL_ADDRESS>Martin
Weinmann<EMAIL_ADDRESS>Stefan Hinz<EMAIL_ADDRESS>
###### Abstract

With FaSS-MVS, we present an approach for fast multi-view stereo with surface-
aware Semi-Global Matching that allows for rapid depth and normal map
estimation from monocular aerial video data captured by UAVs. The data
estimated by FaSS-MVS, in turn, facilitates online $3$D mapping, meaning that
a $3$D map of the scene is immediately and incrementally generated while the
image data is acquired or being received. FaSS-MVS is comprised of a
hierarchical processing scheme in which depth and normal data, as well as
corresponding confidence scores, are estimated in a coarse-to-fine manner,
allowing to efficiently process large scene depths which are inherent to
oblique imagery captured by low-flying UAVs. The actual depth estimation
employs a plane-sweep algorithm for dense multi-image matching to produce
depth hypotheses from which the actual depth map is extracted by means of a
surface-aware semi-global optimization, reducing the fronto-parallel bias of
SGM. Given the estimated depth map, the pixel-wise surface normal information
is then computed by reprojecting the depth map into a point cloud and
calculating the normal vectors within a confined local neighborhood. In a
thorough quantitative and ablative study we show that the accuracies of the
$3$D information calculated by FaSS-MVS is close to that of state-of-the-art
approaches for offline multi-view stereo, with the error not even being one
magnitude higher than that of COLMAP. At the same time, however, the average
run-time of FaSS-MVS to estimate a single depth and normal map is less than
$14$ % of that of COLMAP, allowing to perform an online and incremental
processing of Full-HD imagery at $1$-$2$ Hz.

###### keywords:

Multi-View Stereo , Plane-Sweep Multi-Image Matching , Semi-Global
Optimization , Surface-Awareness , Online Processing , Oblique Aerial Imagery

[1]organization=Fraunhofer Institute of Optronics, System Technologies and
Image Exploitation (IOSB), Fraunhofer Center for Machine Learning,
city=Karlsruhe, postcode=76133, country=Germany [2]organization=Institute of
Photogrammetry and Remote Sensing, Karlsruhe Institute of Technology (KIT),
city=Karlsruhe, postcode=76133, country=Germany

##  1 Introduction

The image-based estimation of depth maps and geometry by dense image matching
(DIM) and multi-view stereo (MVS) is one of the fundamental tasks in
photogrammetry, remote sensing and computer vision, facilitating a wide range
of high-level applications, such as autonomous navigation, urban planning and
monitoring, simulation and 3D modeling, as well as virtual, mixed and
augmented reality. For all these applications, usually no time constraints
exist, meaning that it does not matter whether the processing time is in the
order of a few minutes, a few hours or sometimes even a couple of days, as
long as the final result is as accurate and complete as possible. Especially
for the task of urban monitoring and planning, a large-scale $3$D
reconstruction and modeling is desired, often facilitated by additional data,
e.g. data captured by means of airborne laser scanning (ALS) or terrestrial
laser scanning (TLS), or manual post-processing.

However, the on-going development and increasing availability of commercial
off-the-shelf (COTS) unmanned aerial vehicles (UAVs) open up new possibilities
and applications for image-based $3$D mapping, in both offline and online
processing. In recent years, for example, the use of COTS UAVs by emergency
forces, such as the fire brigade and medical rescue services, heave steadily
increased, facilitating disaster relief or search-and-rescue missions by
allowing a quick and large-scale assessment of the situation or enabling the
monitoring of areas which are inaccessible for ground forces (Restas, 2015;
Furutani and Minami, 2021). In this, image-based techniques and photogrammetry
based on aerial reconnaissance are a key element in supporting the rescue
workers, provided that the environmental conditions, e.g. whether and daytime,
allow for a visual inspection (Furutani and Minami, 2021). There exists a
large collection of software toolboxes, such as COLMAP (Schönberger and Frahm,
2016; Schönberger et al., 2016), for performing offline photogrammetric $3$D
reconstruction allowing to accurately reconstruct the disaster site from
aerial imagery. Their focus, however, is primarily on offline and accurate
processing. This hinders the use for a rapid $3$D mapping during the image
acquisition, due to the required run-time for high-accurate $3$D modeling and
the consideration of all input images for processing. However, in order to
efficiently aid first responders in their tasks, the run-time of the
algorithms matters which, in turn, raises the need for efficient, fast and
incremental $3$D mapping in order to support a rapid assessment. Here, the
availability of $3$D data, for example, allows to reason on damage caused by
an incidence, or structural integrity of a partly collapsed building, as well
as route planning through areas which are difficult to access or to account
for $3$D geometry when creating an orthographic map. In order to accommodate
these applications, we propose a novel approach for fast multi-view stereo
with surface-aware Semi-Global Matching (SGM), denoted as FaSS-MVS. The
approach,

  * 1.

uses plane-sweep sampling to perform hierarchical dense multi-image matching,

  * 2.

utilizes and extends the widely used Semi-Global Matching algorithm
(Hirschmüller, 2005, 2008) to favor not only fronto-parallel surfaces in the
computation of dense depth maps, by incorporating a surface-aware
regularization based on local surface normals,

  * 3.

efficiently computes dense depth, normal and confidence maps from image
sequences, allowing to facilitate the task of incremental UAV-borne $3$D
mapping,

  * 4.

is quantitatively evaluated on two public datasets for dense MVS with accurate
ground truth and that is demonstrated on a use-case specific dataset.

FaSS-MVS combines and extends our previous work presented in Ruf et al. (2017,
2019), by

  * 1.

a more detailed description of the employed algorithms,

  * 2.

extending the plane-sweep multi-image matching for the use of non-fronto-
parallel plane orientations,

  * 3.

improving the surface-aware regularization of the SGM algorithm,

  * 4.

using a different confidence measure for estimation of the confidence map,

  * 5.

a thorough evaluation and ablation study with respect to different aspects and
configurations of the approach,

  * 6.

providing a detailed discussion with respect to support of rescue workers by
aerial reconnaissance.

Even though this approach is proposed with the above-mentioned use-case in
mind, it is not restricted to airborne data and can also be used to perform an
incremental and online $3$D mapping of an environment captured by a ground-
based robot or sensor system.

In the scope of this work, different terminologies with respect to the run-
time capabilities and processing rate of the considered algorithms are used.
Here, _real-time processing_ is used to denote a fast and low latency
processing. This means that the processing rate of the algorithm is high
enough, in order to allow an immediate reaction based on the calculated
results. For example, in the case of reactive collision avoidance based on
depth data from a stereo camera, the calculation of the scene depth needs to
be fast enough to still be able to initiate an appropriate maneuver to avoid
imminent collision. Real-time processing is deliberately not defined by a
minimum processing rate, since the available reaction time depends on various
factors, such as flight speed. In contrast, _online processing_ is used to
denote a fast processing of the algorithm without setting hard time
constraints. Ideally, an algorithm for online processing should be able to
keep up with the frame rate of the input data. At the same time, however, it
is not critical if the algorithm has a high latency and if the results are
only available a couple of frames after the input of corresponding reference
frame. This, for example, is the case when performing depth estimation from
two or more images that are captured by a single moving camera, i.e. monocular
multi-view stereo. Here, enough input images with appropriate baseline have
first to be collected before the actual processing can be started. Both real-
time and online processing refer to a computation during the acquisition or
receiving of the input data. Thus, corresponding algorithms only have a
confined set of input data available during execution. _Offline_ processing,
on the other hand, is done fully disconnected from the actual acquisition of
the input data and, thus, it is assumed that algorithms, which are executed
offline, have access to all available input data.

###  1.1 Paper Outline

This paper is structured as follows: In Section 1.2, the related work on
incremental image-based $3$D mapping for online processing as well as modern
approaches for learning-based DIM and MVS are briefly summarized. In this, it
is also delineated, how the presented approach differs from those presented in
the related work. In Section 2, the overall processing pipeline of the
presented approach is illustrated and outlined with a quick overview. This is
followed by a detailed description on the implementation and methodology of
the individual steps of the processing pipeline. The approach is
quantitatively and qualitatively evaluated on two public and two private
datasets. The datasets, the error metrics as well as the results of the
conducted experiments are presented in Section 3. Subsequently, the findings
are discussed and put into context of the considered use-case in Section 4,
before providing a summary and concluding remarks as well as a short outlook
on future work in Section 5.

###  1.2 Related Work

Due to the ever-increasing demand for detailed $3$D models, the research in
the fields of photogrammetry, remote sensing and computer vision has brought
up a number of software suites and applications, that focus on the estimation
of accurate and dense depth and geometry information from a large set of input
images, by means of DIM and MVS. Prominent and widely used representatives of
such applications are MVE (Goesele et al., 2007), PMVS (Furukawa and Ponce,
2010), SURE (Rothermel et al., 2012; Wenzel et al., 2013b), COLMAP
(Schönberger et al., 2016) and OpenMVS111http://cdcseacave.github.io/openMVS,
to name a few. These approaches, however, are designed for offline processing,
aiming at the accuracy and completeness of the resulting $3$D model, while
assuming that all input data is available at the time of reconstruction and
that no critical constraints on the computation time or hardware resources are
set. In contrast, the aim of FaSS-MVS is to extract dense depth and geometry
information from image sequences, while they are acquired. Or at least while
the image data stream is received if a direct processing is not possible due
to the acquisition by a small UAV and its limited hardware resources for
example. Thus, the focus lies in the incremental and online processing of the
input data by DIM and MVS. In the following, we give a brief overview of the
related work on incremental camera-based mapping for online processing in
Section 1.2.1. Due to the advancements of deep-learning-based approaches in
nearly all fields of computer vision, we also provide a short overview of the
related work on learning-based approaches for DIM and MVS in Section 1.2.2.

####  1.2.1 Incremental Camera-Based Mapping for Online Processing

Early work on incremental and online camera-based mapping of the local
environment was mainly driven by robotic and augmented reality (AR)
applications (Klein and Murray, 2007; Davison et al., 2007; Eade and Drummond,
2006). Here, the main goal was to robustly localize the camera pose, and in
turn the sensor carrier, with respect to its surrounding, in order to navigate
through the environment or enhance the camera images with additional
information. Since the focus of these so-called simultaneous localization and
mapping (SLAM) algorithms is the estimation of the camera pose and trajectory,
the detailed and dense mapping of the environment was rather of secondary
interest. Thus, these approaches mainly relied on point features for the
tracking and mapping rather than direct pixel matching. However, in order to
provide a convincing AR experience, a dense and detailed model of the
environment is essential. Subsequent works (Newcombe and Davison, 2010;
Newcombe et al., 2011) have proposed a dense mapping simultaneous to the
acquisition of the image data and the localization of the camera, resulting in
a detailed reconstruction of a small AR workspace. Since these approaches,
however, aim to reconstruct rather small-scale environments, they make use of
short baseline video clips for the image matching, which in turn allows to
rely on dense optical flow methods to find dense pixel correspondences
(Newcombe and Davison, 2010). In contrast, as input to the approach presented
in this work, it is assumed to have image data captured by a UAV, which is
typically flying several tens of meters away from the object of interest.
Thus, the presented approach is rather aimed to densely map a large-scale
environment, which in turn hinders to track pixel-wise correspondences between
consecutive frames, but requires a wide-baseline image matching instead.
However, the presented approach is not solely restricted to large-scale
environments and a wide baseline image matching, as experiments with respect
to the employed multi-image matching, done in previous work (Ruf et al.,
2017), show.

Early work on camera-based mapping and reconstruction of urban surroundings
was done by Gallup et al. (2007) and Pollefeys et al. (2008), who employed the
plane-sweep algorithm for true multi-image matching, first proposed by Collins
(1996), to map and reconstruct building facades from images captured by a
vehicle-mounted camera in real-time. In this, they rely on vanishing points,
which are detected in the input images, and on data from an additional
inertial measurement unit (IMU) to recover the orientations of the building
facades and the ground plane relative to the camera. To find the optimal plane
configuration for each pixel and, in turn, extract a depth map from the
results of the DIM, Pollefeys et al. (2008) employ a Bayesian formulation with
a subsequent selection of the winner-takes-it-all (WTA) solution, while Gallup
et al. (2007) minimize a formulated energy functional. In (Pollefeys et al.,
2008), the estimation of the camera poses is done by using a Kanade-Lucas-
Tomasi (KLT) feature tracker. A big advantage in the mapping and
reconstruction of urban areas is that most objects in such scenery can be
approximated well by planar structures, which is why plane-sweep DIM is well-
suited for this task. Other approaches for urban reconstruction from ground-
based imagery, like those from Furukawa et al. (2009), Sinha et al. (2009) and
Gallup et al. (2010), perform a piece-wise planar reconstruction by fitting
multiple, differently oriented planes into the scene and optimizing
photometric consistency. In this, they minimize an energy functional by a
graph-cut algorithm, which takes a couple of minutes on a commodity CPU.
Algorithms having a couple of minutes run-time to estimate a single depth map
do not seem to be suitable for fast and online processing at first sight.
However, depending on their ability to be parallelized and optimized for the
execution on a GPU, they might be useful after all.

Around the same time as the previously mentioned work was released,
Hirschmüller (2005, 2008) proposed the so-called SGM algorithm, which evolved
into one of the most widely used approaches for both online and offline DIM,
due its efficiency and convincing results. It has been deployed on both
desktop (Spangenberg et al., 2014; Banz et al., 2011) and embedded (Hernandez-
Juarez et al., 2016; Zhao et al., 2020; Ruf et al., 2021) hardware and is used
in a wide range of applications, such as advanced driver assistance systems
(ADAS) (Spangenberg et al., 2014), real-time obstacle detection and collision
avoidance on-board UAVs (Barry et al., 2015) and urban mapping and
reconstruction from aerial imagery (Rothermel et al., 2012; Wenzel et al.,
2013b; Haala et al., 2015). In their work, Sinha et al. (2014) combine the
plane-sweep multi-image matching with the SGM algorithm to estimate dense and
highly accurate disparity maps. In contrast to the presented approach, Sinha
et al. (2014) use local slanted planes, which are extracted from feature
correspondences, to create disparity hypotheses and employ the SGM algorithm
to recover a disparity map. They, evaluate their approach on a high-resolution
stereo benchmark and achieve significant improvement over the standard SGM
algorithm in both run-time and accuracy. The improvement in terms of run-time
is attributed to the fact that the local plane-sweep allows to test a locally
confined part of the complete disparity range for each pixel, thus reducing
the computational complexity of the optimization within the SGM algorithm.
Similar improvements to overcome the problem of high computational complexity
due to a large disparity range, which is inherent to oblique aerial imagery,
were done by Haala et al. (2015), by embedding the SGM into a hierarchical
coarse-to-fine processing.

Even though a large number of urban environments can be well abstracted by
piecewise planar reconstructions, not all structures are fronto-parallel,
meaning that their surface orientations are not parallel to the image plane.
In order to account for slanted surfaces, Kuschk and Cremers (2013), for
example, have incorporated a second-order smoothness assumption into their
energy function. The initial formulation of the SGM algorithm, however, only
models a first-order smoothness term and thus favors fronto-parallel surfaces,
leading to stair-casing artifacts when reconstructing slanted surfaces.
Especially when aiming for a visually appealing reconstruction of the
environment, this is to be avoided. While Hermann et al. (2009) and Ni et al.
(2018) propose to incorporate a second-order smoothness assumption into the
formulation of the SGM energy function, Scharstein et al. (2017) propose a
more simplistic and yet effective improvement to address this issue. More
specifically, plane priors are used, which, for example, can be recovered from
normal maps or point correspondences, to adjust the zero-cost transition
within the path aggregation of the SGM, thus penalizing deviations from the
surface orientation represented by the prior. The major advantage over the
other approaches is that the pixel-wise offset for the zero-cost transition
can be calculated in advance and is in its magnitude the same for opposite
aggregation paths, making its use very efficient.

####  1.2.2 Learning of Dense Image Matching and Multi-View Stereo
Reconstruction

With the advancements and success of deep-learning-based methods in other
topics of computer vision and photogrammetry, such as object detection,
classification or image segmentation, it was just a matter of time when the
first learning-based approaches for the task of DIM and MVS, that would
outperform state-of-the-art model-based approaches, would be presented. Early
works (Han et al., 2015; Zbontar and LeCun, 2016; Hartmann et al., 2017) use
deep convolutional neural networks (CNNs) to learn similarity measures between
image patches and, in turn, build up a $3$D cost volume from which a disparity
or depth map is extracted by conventional methods, e.g. SGM (Hirschmüller,
2005, 2008).

Early approaches to perform actual MVS with deep learning are the so-called
MVSNet (Yao et al., 2018) and DeepMVS (Huang et al., 2018). Both use the
plane-sweep algorithm to match the pixels of multiple input images based on
learned features and similarity measures and to build up a cost volume, just
like the conventional approaches. To regularize the computed cost volume and
to extract the depth map, both use a $3$D U-Net (Ronneberger et al., 2015).
$3$D CNNs, such as the $3$D U-Net, use a great amount of memory and are
computationally not very efficient, which is why other approaches, such as the
ones presented in (Yao et al., 2019; Yan et al., 2020), exchange the $3$D
U-Net by a cascade of $2$D CNNs. Further approaches (Cheng et al., 2020; Gu et
al., 2020) remedy the high memory consumption by the use of hierarchical
coarse-to-fine processing, as also done in our work. In the construction of
the cost volume there also exist other strategies, such using as gated
convolution (Yi et al., 2020) or reprojecting the image data into a $3$D voxel
grid (Ji et al., 2017).

What all these approaches have in common, however, is that they are trained in
a supervised manner, requiring datasets with appropriate ground truth. Most of
them are using, among others, the DTU MVS benchmark (Jensen et al., 2014),
which also serves as evaluation dataset in the scope of our work. The
availability and versatility of appropriate datasets, however, is not very
high, especially with respect to real-world scenarios, which still greatly
hinders the practical use of deep-learning-based MVS approaches. To overcome
this problem, recent approaches, such as the ones presented in (Khot et al.,
2019; Huang et al., 2021) try to train models in an unsupervised, or sometimes
also denoted as self-supervised, manner. But again, their practical use and
ability for generalization still needs more studies (Khot et al., 2019). These
limitations are the reasons why learning-based approaches for the task of MVS
are not yet practical for the considered use-case, namely to reliably assist
emergency forces in the incremental and online mapping of the operational
area.

In summary, the presented approach adopts a plane-sweep algorithm similar to
the one presented by Pollefeys et al. (2008) to perform efficient dense multi-
image matching and employs an improved implementation of the SGM algorithm to
extract the depth map from the results of the DIM. The use of a plane-sweep
algorithm for the task of DIM is mainly motivated by its ability to create
depth hypotheses by matching an arbitrary number of input images as well as
the fact that it can efficiently be optimized for the massively parallel
execution on GPUs, making particularly suitable for online processing. In the
improved implementation of the SGM algorithm, we, among others, adopt the
approach presented by Scharstein et al. (2017) to account for non-fronto-
parallel surfaces by adjusting the zero-cost transition based on surface
information stored inside a normal map. Very similar to the approach presented
in this work seems to be the approach from Roth and Mayer (2019). They also
rely on the improvements proposed by Scharstein et al. (2017) and combine the
SGM with a plane-sweep DIM. However, their work focuses on the estimation of
disparity images from ground-based stereo image pairs and was only evaluated
on synthetic scenes so far. Moreover, we also propose to reduce the fronto-
parallel bias of the SGM algorithm, by adjusting the zero-cost transition in
the path aggregation based on the gradient of the minimum cost path. The
complete depth estimation pipeline is embedded in a hierarchical processing
scheme, just as proposed by Haala et al. (2015), in order to reduce the
computational complexity induced by a large scene depth inherent to oblique
aerial imagery.

##  2 Materials and Methods

The following sections first give an overview on the full processing pipeline
for an efficient MVS with plane-sweep multi image matching and surface-aware
SGM (Section 2.1). This is followed by a detailed description of the dense
multi-image matching with plane-sweep sampling (Section 2.2) and of the
proposed extensions of the SGM algorithm for surface-aware depth map
extraction (Section 2.3). In addition to the depth map, the presented approach
also computes a normal map (Section 2.4) and a confidence map (Section 2.5).
In a final post-processing step (Section 2.6), further outliers are removed by
masking out untextured regions and checking for geometric consistency.

###  2.1 Processing Pipeline for Fast Multi-View Stereo using Plane-Sweep
Multi-Image Matching and Surface-Aware Semi-Global Optimization

Figure 1: Overview of the proposed approach for incremental MVS with plane-
sweep multi-image matching and surface-aware SGM optimization. Given a bundle
of images and corresponding camera poses
$\left(\mathcal{I},\mathbf{P}\right)_{k}$ of an input sequence, a hierarchical
MVS estimation is performed to recover a depth, normal and confidence map
$\left(\mathcal{D},\mathcal{N},\mathcal{C}\right)$.

An overview of the processing pipeline of our approach for incremental MVS,
based on plane-sweep sampling and a surface-aware SGM optimization (FaSS-MVS)
is illustrated by Figure 1. Given an input bundle
$\left(\mathcal{I},\mathbf{P}\right)_{k}\in\Omega$, consisting of $k$ input
images $\mathcal{I}$, extracted from an image sequence in sequential order, as
well as corresponding camera poses $\mathbf{P}$, our approach computes depth,
normal and confidence maps $\left(\mathcal{D},\mathcal{N},\mathcal{C}\right)$
for a defined reference image $\mathcal{I}_{\mathrm{ref}}$, which is typically
the middle one of the input bundle $\Omega$. In this, we assume that the input
has been calibrated, i.e. that the images are free of lens distortion and that
the full projection matrix
$\mathbf{P}_{k}=\mathbf{K}\left[\mathbf{R}_{k}^{\intercal}\ \
-\mathbf{R}_{k}^{\intercal}\mathrm{C}_{k}\right]$ of each image is known.
Here, $\mathrm{C}_{k}\in\mathbb{R}^{3}$ denotes the locations of the camera
centers with respect to a reference coordinate system
$\mathrm{O}_{\mathrm{ref}}$, while the column vectors of
$\mathbf{R}^{\intercal}_{k}\in SO(3)$ hold the normalized coordinate axes of
the camera coordinate system $\mathrm{O}_{\mathrm{cam}}$, as seen from
$\mathrm{O}_{\mathrm{ref}}$. The intrinsic calibration matrix $\mathbf{K}$ is
equal for all cameras, since it is assumed that the images are captured from a
single camera. Typical for an MVS approach, it is assumed that the input
images depict the scene which is to be reconstructed from slightly different
viewpoints.

Before any processing, a Gaussian image pyramid with $n$ pyramid levels is
computed for each image of the input bundle, allowing a subsequent
hierarchical processing. The lowest pyramid levels ($l=0$) hold the input
images with their original image size. While moving up the pyramid, the images
are first blurred with a Gaussian filter of size of $3\times 3$ pixels and
with $\sigma=1$, before being scaled down by a factor of $0.5$ in both image
directions. Furthermore, the intrinsic calibration matrices are also adjusted
by halving the focal length and the coordinates of the principal point, in
order to account for the reduced image size. This results in an augmentation
of the input bundle $\Omega$ by $n-1$ additional sets. In the following, a
superscript is used to mark the results and processes at a specific pyramid
level. The pipeline is initialized at the coarsest pyramid level $l=(n-1)$
with the smallest image size, executing three subsequent computational parts
at each pyramid level, thus resulting in a coarse-to-fine processing.

The first part of the actual processing, namely the depth estimation, computes
a depth map $\mathcal{D}^{l}$ and is, in turn, subdivided into a plane-sweep
multi-image matching creating depth hypotheses and the SGMx optimization
extracting the optimal depth from the set of hypotheses. The latter one adopts
the SGM approach, which was first presented by Hirschmüller (2005, 2008), to
the plane-sweep matching, and extends it to account for non-fronto-parallel
surface structures. Section 2.2 gives a detailed description on the employed
plane-sweep algorithm, while the extension of the SGM algorithm is discussed
in Section 2.3. A concluding depth refinement and median filter with a kernel
size of $5\times 5$ pixels is used to filter small outliers in the resulting
depth map.

The second part of the hierarchical processing estimates a normal map
$\mathcal{N}^{l}$ from the previously computed depth map $\mathcal{D}^{l}$.
Apart from being an additional output of FaSS-MVS, the normal map is also used
as part of the SGMx optimization to account for the surface orientation in the
next hierarchical iteration. The normal map is regularized by an appearance-
based weighted Gaussian smoothing in order to smooth out small variations
while preserving discontinuities. Section 2.4 describes the details on how the
normal map is extracted from a single depth map.

In the third part, a confidence map $\mathcal{C}^{l}$ is computed, holding
pixel-wise confidence scores in the interval of $\left[0,1\right]$ with
respect to the depth estimates in $\mathcal{D}^{l}$. The final confidence
scores are computed based on the surface orientation at the considered pixel.
Details on the computation of $\mathcal{C}^{l}$ are given in Section 2.5.

Inherent to a hierarchical coarse-to-fine processing, as long as the lowest
levels of the image pyramids are not yet reached, i.e. while $l>0$, the depth
map $\mathcal{D}^{l}$ and normal map $\mathcal{N}^{l}$ computed at level $l$
are used to initialize the depth map estimation at the next pyramid level
$l-1$. In this, $\mathcal{D}^{l}$ and $\mathcal{N}^{l}$ are upscaled with
nearest neighbor interpolation to the image size of the next pyramid level,
yielding $\bar{\mathcal{D}}^{l}$ and $\bar{\mathcal{N}}^{l}$. Then,
$\bar{\mathcal{D}}^{l}$ is first used to compute the pixel-wise sampling range
$\Gamma^{\,l-1}_{\mathrm{p}}$ of the multi-image plane-sweep algorithm at the
next pyramid level. Here, the $\Gamma^{\,l-1}_{\mathrm{p}}$ is computed for
each pixel $\mathrm{p}$ separately, based on the previous depth estimate
$\bar{d}^{\,l}_{\,\mathrm{p}}=\bar{D}^{\,l}(\mathrm{p})$ and a predefined
window with a radius of $\Delta d$ around $\bar{d}^{\,l}_{\,\mathrm{p}}$:

$\displaystyle\Gamma^{\,l-1}_{\mathrm{p}}$
$\displaystyle=\left[d^{\,l-1}_{\,\mathrm{p,min}},d^{\,l-1}_{\,\mathrm{p,max}}\right],\text{with}$
(1) $\displaystyle d^{\,l-1}_{\,\mathrm{p,min}}$
$\displaystyle=\bar{d}^{\,l}_{\,\mathrm{p}}-\Delta d,$ $\displaystyle
d^{\,l-1}_{\,\mathrm{p,max}}$
$\displaystyle=\bar{d}^{\,l}_{\,\mathrm{p}}+\Delta d.$

In the first iteration, the sampling range is set equally for all pixels and
is parameterized by minimum and maximum scene depth:
$\Gamma=\left[d_{\,\mathrm{min}},d_{\,\mathrm{max}}\right]$. The upscaled
normal map $\bar{\mathcal{N}}^{l}$ is used by one of the proposed SGM
extensions to account for the surface orientation within the scene. The final
depth, normal and confidence maps are the outcome of the processing at the
lowest pyramid level. They are denoted as $\mathcal{D}$, $\mathcal{N}$ and
$\mathcal{C}$ respectively, and have the same image size as the input images.

In a final post-processing step, we use a Difference-of-Gaussian (DoG) filter
(Wenzel, 2016) to unmask image regions, which do not have distinctive texture
information that allows to perform a reliable matching, as well as a geometric
consistency check. Details on the final post-processing are given in Section
2.6.

###  2.2 Real-Time Dense Multi-Image Matching with Plane-Sweep Sampling

Figure 2: Overview of the algorithm for plane-sweep multi-image matching. A
scene is sampled by a plane $\Pi=(\mathrm{n},\delta)$, with $\mathrm{n}$ being
the normal vector of the plane and $\delta$ being the orthogonal distance of
the plane from $\mathrm{C}_{\mathrm{ref}}$, that is swept along its normal
vector between two bounding planes $\Pi_{\mathrm{max}}$ and
$\Pi_{\mathrm{min}}$ through space. For each distance $\delta$ of $\Pi$ the
reference pixel $\mathrm{p}^{\mathrm{ref}}$ is projected by the plane induced
homography $\mathbf{H}_{\mathrm{ref}\rightarrow k}$ into arbitrary number of
viewpoints, where it is matched against the corresponding pixel in
$\mathcal{I}_{k}$.

Dense image matching refers to the process of finding a dense correspondence
field between pixels of two or more images from which depth hypotheses are
extracted. In contrast to a sparse correspondence field, which is typically
only established for discriminative feature points, the aim of dense image
matching is to find correspondences for a majority of pixels in the reference
image, usually relying on a photometric similarity measure. If the relative
camera poses of input images are known, these correspondences can further be
transformed into a depth map, which encodes the pixel-wise scene depth from
the vantage point of the reference camera for which the correspondence field
and depth map are computed. Furthermore, if intrinsic information of the
reference camera as well as an absolute position with respect to a reference
coordinate system are known, the depth map can be transformed into a point
cloud in three-dimensional space, revealing the structure of the scene
(Remondino et al., 2013).

When considering the case of a calibrated stereo setup with two views, the
fundamental matrix of the epipolar geometry allows for each pixel in the
reference image to reduce the search space for the pixel correspondence in the
matching image to the corresponding one-dimensional epipolar line, greatly
reducing the complexity. The epipolar constraint is a key characteristic of
the two-view geometry. A similar characteristic exists also for the three- and
four-view geometry. Analogously to the fundamental matrix, the trifocal and
quadrifocal tensors allow to establish point-line correspondences between
three and four views, respectively. Apart from the increasing complexity in
the computation of the tensors with the number of views due to more degrees of
freedom, their use, however, does not go beyond four views (Hartley and
Zisserman, 2004).

Besides the fundamental matrix, there exists a second projective relationship
between two views, which can be extended equally to an arbitrary number of
views. Again, we consider the case of a two-camera setup with known intrinsic
and extrinsic parameters. This time, however, an additional scene plane
$\Pi=\left(\mathrm{n},\delta\right)$ is positioned the field of view of both
cameras. Here, the scene plane $\Pi$ parameterized by its normal vector
$\mathrm{n}$ and its distance $\delta$ from the first camera, i.e. the
reference camera. Given this setup, the image point
$\mathrm{p}^{\mathrm{ref}}$ in reference camera image can directly be mapped
onto the image point $\mathrm{p}^{k}$ in the other camera image via the
homography $\mathbf{H}$ induced by the plane $\Pi$ according to

$\displaystyle\mathrm{p}^{k}$
$\displaystyle=\mathbf{H}\left(\Pi,\mathbf{P}_{\mathrm{ref}},\mathbf{P}_{k}\right)\cdot\mathrm{p}^{\mathrm{ref}}\
,\ \ \text{with}$ (2)
$\displaystyle\mathbf{H}\left(\Pi,\mathbf{P}_{\mathrm{ref}},\mathbf{P}_{k}\right)$
$\displaystyle=\mathbf{K}_{k}\cdot\frac{\mathbf{R}-\mathrm{t}\mathrm{n}^{\intercal}}{\delta}\cdot\mathrm{K}_{\mathrm{ref}}^{\text{-}1}.$

Here, $\mathbf{K}_{\mathrm{ref}}$ and $\mathbf{K}_{k}$ denote the intrinsic
matrices of both cameras, and $\left[\mathbf{R}\ \mathrm{t}\right]$ denotes
the relative transformation matrix of the neighbor pose $\mathbf{P}_{k}$ with
respect to reference pose $\mathbf{P}_{\mathrm{ref}}$. Equation 2 is
geometrically interpreted by casting a viewing ray through the image point
$\mathrm{p}^{\mathrm{ref}}$ and intersecting it with the scene plane $\Pi$,
yielding a scene point $\mathrm{P}_{\Pi}$, which is then projected into the
second camera, resulting in the image point $\mathrm{p}^{k}$ (Hartley and
Zisserman, 2004).

####  2.2.1 The Hierarchical Plane-Sweep Algorithm for Real-Time Multi-Image
Matching

Data: a calibrated image bundle $\Omega^{l}$ at the pyramid level $l$, and a
set of planes $\Pi$ with a normal vector $\mathrm{n}$ and varying distances
$\delta$ as well as a local depth sampling range
$\Gamma^{\,l}_{\mathrm{p}}=[d^{\,l}_{\,\mathrm{p,min}},d^{\,l}_{\,\mathrm{p,max}}]$.

Result: three-dimensional cost volume $\mathcal{S}$, holding the pixel-wise
matching score for each pixel
$\mathrm{p}^{\mathrm{ref}}\in\mathcal{I}^{l}_{\mathrm{ref}}$ and plane $\Pi$.

1

2determine bounding planes $\Pi_{\mathrm{min}}$ and $\Pi_{\mathrm{max}}$
located at $\delta_{\mathrm{min}}$ and $\delta_{\mathrm{max}}$, so that the
local depth range $\Gamma^{\,l}_{\mathrm{p}}$ is completely sampled (cf.
Section 2.2.2).

3foreach  _pixel $\mathrm{p}^{\mathrm{ref}}\in\mathcal{I}^{l}_{\mathrm{ref}}$
and distance
$\delta\in\left[\delta_{\mathrm{min}},\delta_{\mathrm{max}}\right]$_ do

4

5 configure scene plane $\Pi=(\mathrm{n},\delta)$.

6 determine pixels $\mathrm{p}^{\,k}$ in all matching images
$\mathcal{I}^{l}_{k}\in\Omega^{l}\setminus\mathcal{I}^{l}_{\mathrm{ref}}$ :
$\mathrm{p}^{\,k}=\mathbf{H}\left(\Pi,\mathbf{P}^{l}_{\mathrm{ref}},\mathbf{P}^{l}_{k}\right)\cdot\mathrm{p}^{\mathrm{ref}}.$

7 warp local image patches $\mathcal{P}^{l}_{k}\in\mathcal{I}^{l}_{k}$ around
$\mathrm{p}^{\,k}$, with the same size as the support region of the matching
cost function $C(\cdot)$, into $\mathcal{I}^{l}_{\mathrm{ref}}$ :
$\tilde{\mathcal{P}}^{l}_{k}=\mathbf{H}\left(\Pi,\mathbf{P}^{l}_{\mathrm{ref}},\mathbf{P}^{l}_{k}\right)^{\text{-}1}\cdot\mathcal{P}^{l}_{k}.$

8 compute the matching cost $s\left(\mathrm{p},\Pi\right)$ between reference
patch $\mathcal{P}^{l}_{\mathrm{ref}}\in\mathcal{I}^{l}_{\mathrm{ref}}$ and
$\tilde{\mathcal{P}}^{l}_{k}$ for left and right subset of cameras separately
:
$s^{\mathrm{L}}\left(\mathrm{p},\Pi\right)=\sum_{k<\mathrm{ref}}C\left(\mathcal{P}^{l}_{\mathrm{ref}},\tilde{\mathcal{P}}^{l}_{k}\right),$
$s^{\mathrm{R}}\left(\mathrm{p},\Pi\right)=\sum_{k>\mathrm{ref}}C\left(\mathcal{P}^{l}_{\mathrm{ref}},\tilde{\mathcal{P}}^{l}_{k}\right).$

9 store the minimum of left and right matching cost (accounting for occlusions
as described by Kang et al. (2001)) into three-dimensional cost volume
$\mathcal{S}$ :
$\mathcal{S}^{l}\left(\mathrm{p},\Pi\right)=\min\\{s^{\mathrm{L}}\left(\mathrm{p},\Pi\right),s^{\mathrm{R}}\left(\mathrm{p},\Pi\right)\\}.$

10 end foreach

Algorithm 1 Plane-sweep multi-image matching executed at a specific pyramid
level $l$ of the proposed hierarchical processing scheme.

Based on the relationship between two cameras and a scene plane, Collins
(1996) proposed an algorithm for true multi-image matching. This algorithm
samples the scene-space between two bounding planes $\Pi_{\mathrm{min}}$ and
$\Pi_{\mathrm{max}}$, located at $\delta_{\mathrm{min}}$ and
$\delta_{\mathrm{max}}$, by sweeping a plane along its normal vector
$\mathrm{n}$ through space and matching the input images according to Equation
2 for each distance $\delta\leavevmode\nobreak\ \in\leavevmode\nobreak\
\left[\delta_{\mathrm{min}},\leavevmode\nobreak\ \delta_{\mathrm{max}}\right]$
of the plane relative to the reference camera. For each position of the plane,
an arbitrary number of matching images are warped by the plane-induced
homography $\mathbf{H}^{\text{-}1}_{\mathrm{ref}\rightarrow k}$ into the view
of the reference camera, where they are matched against the reference image.
If the scene plane is close to a three-dimensional structure, then the
corresponding image regions of the warped matching images overlap with the
projection in the reference image, allowing to deduce the scene depth of the
corresponding object from the parameterization of the corresponding plane (cf.
Figure 2). First denoted as space-sweep algorithm, it was adopted by numerous
studies on multi-image matching and MVS (Gallup et al., 2007; Pollefeys et
al., 2008; Sinha et al., 2014), eventually denoting it as plane-sweep
algorithm. The presented approach for hierarchical multi-image matching is
based on the plane-sweep algorithm presented by Pollefeys et al. (2008) and it
is described in Algorithm 1.

As part of the actual image matching, the Hamming distance of the census
transform (CT) (Zabih and Woodfill, 1994) as well as a negated, truncated and
scaled form of the normalized cross-correlation (NCC) (Scharstein et al.,
2017; Sinha et al., 2014) are employed and evaluated as cost function
$C(\cdot)$. And since the approach considers a bundle of input images with an
equal number of matching images to either sides of the reference image, the
approach presented by Kang et al. (2001) is adopted, using the minimum
aggregated matching cost of the left and right subset of the matching images
in order to account for occlusions. The resulting three-dimensional cost
volume $\mathcal{S}^{\,l}$ is of size $w^{\,l}\times
h^{\,l}\times|\delta^{\,l}|$, with $w^{\,l}$ and $h^{\,l}$ being the width and
height of the reference image and $|\delta^{\,l}|$ being the number plane
positions at which the matching is performed, all with respect to the current
pyramid level $l$. In this, the cost volume $\mathcal{S}^{\,l}$ is implemented
as a dynamic cost volume (Haala et al., 2015) for all but the topmost pyramid
level, since the sampling range $\Gamma^{\,l}_{\mathrm{p}}$ is determined
independently for each pixel $\mathrm{p}$. Nonetheless, the complete set of
plane distances $\delta\in[\delta_{\mathrm{min}},\delta_{\mathrm{max}}]$,
deduced from $\Gamma$, are precomputed for each pyramid level $l$ and the same
for all pixels. This, in turn, allows to precompute the homographic mappings
for all planes $\Pi$. In the following, we describe how we find the bounding
planes $\Pi_{\mathrm{min}}$ and $\Pi_{\mathrm{max}}$ and the corresponding
distances $\delta_{\mathrm{min}}$ and $\delta_{\mathrm{max}}$ (Section 2.2.2),
as well as how we use the cross-ratio to find appropriate sampling steps
within $\Gamma$ (Section 2.2.3).

####  2.2.2 Determining the Bounding Planes Corresponding to the Given Depth
Range

As previously described, it is assumed that two bounding planes, namely
$\Pi_{\mathrm{min}}$ and $\Pi_{\mathrm{max}}$ with corresponding distances
$\delta_{\mathrm{min}}$ and $\delta_{\mathrm{max}}$, in between which the
scene is to be sampled, are known. In case of a fronto-parallel sampling
strategy, i.e. $\mathrm{n}=(0\ 0\ \text{-}1)^{\intercal}$ with respect to the
local camera coordinate system, the distances $\delta_{\mathrm{min}}$ and
$\delta_{\mathrm{max}}$ are equal to the minimum and maximum depth, namely
$d_{\mathrm{min}}$ and $d_{\mathrm{max}}$. This, however, does not hold for
non-fronto-parallel plane orientations. To find the bounding planes for
slanted plane orientations, a view-frustum, corresponding to the reference
camera, is first constructed for which the depth is to be estimated. This
view-frustum is represented by a pyramid that resembles the field-of-view of
the camera and that is truncated by two fronto-parallel near and far planes
that are located at $d_{\mathrm{min}}$ and $d_{\mathrm{max}}$. Given the four
corner points of the view-frustum on the near plane
$\mathrm{X}^{\mathrm{near}}_{i}$ and the four on the far plane
$\mathrm{X}^{\mathrm{far}}_{i}$, the minimum and maximum distance
$\delta_{\mathrm{min}}$ and $\delta_{\mathrm{max}}$ are found according to:

$\displaystyle\delta_{\mathrm{min}}$
$\displaystyle=\min\limits_{i}(|\mathrm{n}^{\intercal}\cdot\mathrm{X}^{\mathrm{near}}_{i}|)\
,\ \ \text{and}$ (3) $\displaystyle\delta_{\mathrm{max}}$
$\displaystyle=\min\limits_{i}(|\mathrm{n}^{\intercal}\cdot\mathrm{X}^{\mathrm{far}}_{i}|).$

In order to avoid an orientation flip of the image data, all camera centers
$\mathrm{C}_{i}$ need to lie in front of $\Pi_{\mathrm{min}}$ with respect to
the sweeping direction, thus for all camera centers
$\mathrm{n}^{\intercal}\cdot\mathrm{C}_{i}+\delta_{\mathrm{min}}>0$ must hold.

####  2.2.3 Finding the Sampling Steps by Utilizing the Cross-Ratio

As stated by Equation 2, the sampling planes $\Pi$ of the plane-sweep
algorithm are parameterized by two parameters, namely the normal vector
$\mathrm{n}$, denoting the orientation and sweeping direction of the plane,
and the orthogonal distance $\delta$ from the optical center of the reference
camera $\mathrm{C}_{\mathrm{ref}}$. While $\mathrm{n}$ allows to adjust the
warping of the images and, thus, the image matching to the surface orientation
within the scene, the second parameter $\delta$ determines the step-size with
which the scene is sampled. In this, a straight-forward approach would be to
select the step-size in such a way that the scene is sampled with a desired
resolution, i.e. sweeping the planes at equidistant unit intervals through the
scene space. However, it is not guaranteed that a thorough sampling of the
scene with a small step-size results in a higher accuracy. If the step-size is
not chosen in accordance with the camera positions of the input images and the
baseline between the cameras, the matching results of two or more consecutive
plane positions might not reveal enough difference and, thus, introduce
ambiguities between multiple plane hypotheses. Furthermore, in terms of
efficiency, it is important to vary the sampling rate in scene space with
respect to the distance of the plane relative to the reference camera, since
the perspective projection requires an increasingly smaller step-size as the
plane moves closer to the camera.

Data: two cameras with full projection matrices $\mathbf{P}_{\mathrm{ref}}$
and $\mathbf{P}_{k}$, an image point $\mathrm{p}^{\mathrm{ref}}$ inducing
largest disparity when warped from $\mathcal{I}_{\mathrm{ref}}$ to
$\mathcal{I}_{k}$, as well as two bounding planes $\Pi_{\mathrm{min}}$ and
$\Pi_{\mathrm{max}}$.

Result: list of orthogonal plane distances $\delta$ relative to
$\mathrm{C}_{\mathrm{ref}}$, such that the maximum pixel-displacement between
the warped images of two consecutive planes is less than or equal to $1$.

1 calculate the viewing ray $\mathrm{V}^{\mathrm{ref}}_{\mathrm{p}}$, going
through $\mathrm{C}_{\mathrm{ref}}$ and $\mathrm{p}^{\mathrm{ref}}$, and
intersect it with $\Pi_{\mathrm{min}}$ and $\Pi_{\mathrm{max}}$, yielding the
scene points $\mathrm{P}_{\mathrm{min}}$ and $\mathrm{P}_{\mathrm{max}}$.

2project the optical center $\mathrm{C}_{\mathrm{ref}}$, as well as
$\mathrm{P}_{\mathrm{min}}$ and $\mathrm{P}_{\mathrm{max}}$ onto the image
plane of the second camera, yielding the epipole
$\mathrm{e}^{k}_{\mathrm{ref}}$ and the two image points
$\mathrm{p}^{k}_{\mathrm{min}}$ and $\mathrm{p}^{k}_{\mathrm{max}}$, all lying
on the epipolar line $\mathrm{l}^{k}_{\mathrm{p}}$.

3determine the unit vector
$\mathrm{k}=\frac{\mathrm{p}^{k}_{\mathrm{min}}-\mathrm{p}^{k}_{\mathrm{max}}}{||\mathrm{p}^{k}_{\mathrm{min}}-\mathrm{p}^{k}_{\mathrm{max}}||}$,
being the normalized direction of $\mathrm{l}^{k}_{\mathrm{p}}$ and pointing
from $\mathrm{m}^{k}_{\mathrm{max}}$ to $\mathrm{m}^{k}_{\mathrm{min}}$.

4for  _$\mathrm{p}^{k}_{i}\leftarrow\mathrm{p}^{k}_{\mathrm{max}}$ to
$\mathrm{p}^{k}_{\mathrm{min}}$ by
$\mathrm{p}^{k}_{i+1}=\mathrm{p}^{k}_{i}+\mathrm{k}$_ do

5 given the viewing rays $\mathrm{V}^{k}_{\mathrm{e}_{\mathrm{ref}}}$,
$\mathrm{V}^{k}_{\mathrm{p}_{\mathrm{min}}}$,
$\mathrm{V}^{k}_{\mathrm{p}_{i}}$ and
$\mathrm{V}^{k}_{\mathrm{p}_{\mathrm{max}}}$ going through the optical center
of $\mathrm{C}_{k}$ and $\mathrm{e}^{k}_{\mathrm{ref}}$,
$\mathrm{p}^{k}_{\mathrm{min}}$, $\mathrm{p}^{k}_{i}$ and
$\mathrm{p}^{k}_{\mathrm{max}}$ respectively, apply Equation 23 and Equation
24 to compute $\mathrm{P}_{i}\in\mathrm{V}^{\mathrm{ref}}_{\mathrm{p}}$
according to:  $\displaystyle
Q(\mathrm{V}^{k}_{\mathrm{e}_{\mathrm{ref}}},\mathrm{V}^{k}_{\mathrm{p}_{\mathrm{min}}},\mathrm{V}^{k}_{\mathrm{p}_{i}},\mathrm{V}^{k}_{\mathrm{p}_{\mathrm{max}}})$
$\displaystyle=\frac{\sin(\alpha(\mathrm{V}^{k}_{\mathrm{e}_{\mathrm{ref}}},\mathrm{V}^{k}_{\mathrm{p}_{i}}))\cdot\sin(\alpha(\mathrm{V}^{k}_{\mathrm{p}_{\mathrm{min}}},\mathrm{V}^{k}_{\mathrm{p}_{\mathrm{max}}}))}{\sin(\alpha(\mathrm{V}^{k}_{\mathrm{e}_{\mathrm{ref}}},\mathrm{V}^{k}_{\mathrm{p}_{\mathrm{max}}}))\cdot\sin(\alpha(\mathrm{V}^{k}_{\mathrm{p}_{\mathrm{min}}},\mathrm{V}^{k}_{\mathrm{p}_{i}}))}$
$\displaystyle=\frac{\Delta(\mathrm{C}_{\mathrm{ref}},\mathrm{P}_{i})\cdot\Delta(\mathrm{P}_{\mathrm{min}},\mathrm{P}_{\mathrm{max}})}{\Delta(\mathrm{C}_{\mathrm{ref}},\mathrm{P}_{\mathrm{{max}}})\cdot\Delta(\mathrm{P}_{\mathrm{min}},\mathrm{P}_{i})}.$

6 since
$Q(\mathrm{C}_{\mathrm{ref}},\mathrm{P}_{\mathrm{min}},\mathrm{P}_{i},\mathrm{P}_{\mathrm{max}})=Q(\mathrm{C}_{\mathrm{ref}},\delta_{\mathrm{min}},\delta,\delta_{\mathrm{max}})$,
derive $\delta$ relative to $\mathrm{C}_{\mathrm{ref}}$ according to:
$\displaystyle\frac{\delta\cdot\left(\delta_{\mathrm{max}}-\delta_{\mathrm{min}}\right)}{\delta_{\mathrm{max}}\cdot\left(\delta-\delta_{\mathrm{min}}\right)}$
$\displaystyle=\frac{\sin(\alpha(\mathrm{V}^{k}_{\mathrm{e}_{\mathrm{ref}}},\mathrm{V}^{k}_{\mathrm{p}_{i}}))\cdot\sin(\alpha(\mathrm{V}^{k}_{\mathrm{p}_{\mathrm{min}}},\mathrm{V}^{k}_{\mathrm{p}_{\mathrm{max}}}))}{\sin(\alpha(\mathrm{V}^{k}_{\mathrm{e}_{\mathrm{ref}}},\mathrm{V}^{k}_{\mathrm{p}_{\mathrm{max}}}))\cdot\sin(\alpha(\mathrm{V}^{k}_{\mathrm{p}_{\mathrm{min}}},\mathrm{V}^{k}_{\mathrm{p}_{i}}))}.$

7 end for

Algorithm 2 Finding plane distances $\delta$ by utilizing the cross-ratio.

Thus, a common approach is to select the sampling positions of the planes
according to the disparity change induced by two consecutive planes. In this,
the pixel-wise motion between the warped images of two consecutive planes
should not exceed an absolute value of $1$ (Szeliski and Scharstein, 2004;
Pollefeys et al., 2008). To ensure that the maximum disparity change between
the warped images of two consecutive planes is less or equal to $1$, Pollefeys
et al. (2008) evaluate the displacements occurring at the boundaries of the
most distant images for a set of planes with predefined parameters and only
select those that fulfill the stated criterion for the actual image matching.
Thus, for each plane within the predefined set, an additional test is
performed, involving the warping of image points at the boundaries, in order
to determine whether the plane is suitable or not.

In contrast, with this approach, we aim to directly derive the distances of
the sampling planes from correspondences in image space, implicitly leading to
an inverse depth sampling in scene space. Given an image point
$\mathrm{p}^{\mathrm{ref}}\in\mathcal{I}_{\mathrm{ref}}$ in the reference
image, an intuitive approach would be to select multiple sampling points
$\mathrm{p}^{k}\in\mathcal{I}_{k}$ on the corresponding epipolar line
$\mathrm{l}^{k}_{\mathrm{p}}$ in one of the other cameras, and find the
corresponding plane distances by triangulation between
$\mathrm{p}^{\mathrm{ref}}$ and $\mathrm{p}^{k}$. The effort of triangulation,
however, can be avoided by relying on the cross-ratio which is invariant under
perspective projection (cf. A).

Given the projection matrices $\mathbf{P}_{\mathrm{ref}}$ and $\mathbf{P}_{k}$
of two cameras, as illustrated in Figure 3, we center the coordinate system in
the optical center of the reference camera $\mathrm{C}_{\mathrm{ref}}$, so
that $\mathbf{P}_{\mathrm{ref}}=\mathbf{K}\left[\mathbf{I}\ 0\right].$ Thus,
we are aiming to find the plane distances relative to
$\mathrm{C}_{\mathrm{ref}}$. In case of MVS with multiple cameras, for
$\mathrm{C}_{k}$, we select the camera which will induce the largest image
offset, and thus giving an upper bound on the disparity range. As already
noted by Pollefeys et al. (2008), this is typically the camera which is most
distant from the reference camera.

Figure 3: Determination of the orthogonal distance parameter of the sampling
planes of the plane-sweep multi-image matching by using the cross-ratio and
epipolar geometry. Here, $\mathrm{C}_{\mathrm{ref}}$ and $\mathrm{C}_{k}$
represent the positions of the optical centers of the two cameras.

Furthermore, we again assume to have two bounding planes $\Pi_{\mathrm{max}}$
and $\Pi_{\mathrm{min}}$, which limit the sweep space of the planes. Just as
in the work of Pollefeys et al. (2008), it is ensured that the sweep space
does not intersect the convex hull of the cameras, in order to avoid
inversions in the process of image matching (cf. Section 2.2.2). We choose
$\mathrm{p}^{\mathrm{ref}}$ as the pixel which induces the largest disparity
when warped from $\mathcal{I}_{\mathrm{ref}}$ to $\mathcal{I}_{k}$ via
$\mathbf{H}\left(\Pi_{\mathrm{min}},\mathbf{P}_{\mathrm{ref}},\mathbf{P}_{k}\right)$,
typically one of the four corners, in order to guarantee a maximum disparity
change between successive planes and find a set of sampling planes, with
distances $\delta$ relative to $\mathrm{C}_{\mathrm{ref}}$ according to
Algorithm 2 and illustrated by Figure 3.

This approach is computationally efficient and not restricted to a fronto-
parallel orientation of the sampling planes as long as the optical axis of the
reference camera intersects with the planes and the sweeping vector has a
component that is parallel to the optical axis. Furthermore, in order to
accommodate for all possible setups of $\mathrm{C}_{\mathrm{ref}}$ and
$\mathrm{C}_{k}$, it is important to use
$Q(\mathrm{V}^{k}_{\mathrm{e}_{\mathrm{ref}}},\mathrm{V}^{k}_{\mathrm{p}_{\mathrm{min}}},\mathrm{V}^{k}_{\mathrm{p}_{i}},\mathrm{V}^{k}_{\mathrm{p}_{\mathrm{max}}})$
in Algorithm 2, since $\mathrm{e}^{k}_{\mathrm{ref}}$ would flip to the side
of $\mathrm{p}^{k}_{\mathrm{max}}$ if the focal plane of reference camera is
behind $\mathrm{C}_{k}$.

###  2.3 Depth Map Computation with Surface-Aware Semi-Global Matching

The hierarchical plane-sweep algorithm for multi-image matching, as described
in Section 2.2, produces at each pyramid level a three-dimensional cost volume
$\mathcal{S}^{l}(\mathrm{p},\Pi)$, which holds matching costs for each pixel
$\mathrm{p}\in\mathcal{I}_{\mathrm{ref}}$ induced by a given plane $\Pi$
located at distance $\delta$ orthogonal to the optical center
$\mathrm{C}_{\mathrm{ref}}$ of the reference camera. In the second stage of
the depth estimation within FaSS-MVS, the cost volume is regularized by a
semi-global optimization scheme, yielding a dense depth map $\mathcal{D}^{l}$.
Building upon the original SGM approach (Hirschmüller, 2005, 2008), we propose
three different optimization schemes (SGMx). Apart from a straight-forward
adaptation of the SGM approach to the plane-sweep sampling, we also adopt the
approach of Scharstein et al. (2017) to also favor slanted surfaces by
considering surface information available in the form of surface normals.
Furthermore, we investigate a third extension, which penalizes deviations from
the gradient of the minimum cost path within the SGM optimization scheme. In
the following, we first recap the SGM algorithm in Section 2.3.1, before
describing the proposed extension in detail in Section 2.3.2.

####  2.3.1 Semi-Global Matching

In early studies on dense image matching (DIM), Scharstein and Szeliski (2002)
have grouped the stereo algorithms into three groups based on their
optimization strategy, i.e. the approach on how to extract the final disparity
or depth map from the depth hypotheses produced by the image matching. The
first group contains all algorithms, which employ a local optimization scheme.
In this, a locally confined support region is used to find the pixel-wise
optimal depth, i.e. the WTA solution, from the cost volume. While local
algorithms are computationally very efficient, they only model an implicit
smoothness assumption within the local window. Global methods on the other
hand, which make up the second category of the taxonomy presented by
Scharstein and Szeliski (2002), formulate explicit smoothness assumptions as
part of an energy function, which is to be minimized in a global optimization
scheme. This, in turn, however, makes such algorithms computationally more
expensive, while at the same time producing more accurate results than local
methods.

The third category, to which the Semi-Global Matching (SGM) algorithm of
Hirschmüller (2005, 2008) belongs to, is a subgroup to the global methods. It
contains algorithms, which explicitly formulate the stereo problem in a global
energy function, yet employ dynamic programming to approximate the computation
of the optimal solution. In dynamic programming, the original optimization
problem is broken down into a number of smaller, less complex problems, which
can be processed and solved independently of each other. The solution of the
initial problem is then the sum of the solutions of the sub-problems. This
makes such algorithms computationally very efficient, while at the same time
resulting in a globally consistent and accurate solution.

In his work, Hirschmüller (2005, 2008) proposed the SGM algorithm for the task
of disparity estimation as part of the stereo normal case. In this, the stereo
problem is formulated as a two-dimensional Markov Random Field (MRF) with the
energy function $E(\mathcal{U})$, which is to be minimized in order to find
the disparity map $\mathcal{U}$:

$\begin{split}E(\mathcal{U})=\sum_{\mathrm{p}}\Big{(}\mathcal{S}(\mathrm{p},u_{\mathrm{p}})&+\sum_{\mathrm{q}\in\mathcal{W}_{\mathrm{p}}}\varphi_{1}\cdot[|u_{\mathrm{p}}-u_{\mathrm{q}}|=1]\\\
&+\sum_{\mathrm{q}\in\mathcal{W}_{\mathrm{p}}}\varphi_{2}\cdot[|u_{\mathrm{p}}-u_{\mathrm{q}}|>1]\Big{)}.\end{split}$
(4)

here, $[\cdot]$ denotes the Iverson bracket, penalizing an absolute deviation
of $1$ between disparity $u$ of pixel $\mathrm{p}$ and a neighboring pixel
$\mathrm{q}\in\mathcal{W}_{\mathrm{p}}$ by the penalty $\varphi_{1}$, and a
deviation greater than $1$ by $\varphi_{2}$, with $\varphi_{1}<\varphi_{2}$.

The minimization of Equation 4, however, is efficiently approximated by the
utilization of dynamic programming. For each pixel, the matching costs inside
the cost volume $\mathcal{S}$ are aggregated along numerous concentric one-
dimensional paths with direction $\mathrm{r}$, which can be processed fully
independently of each other. In this, the matching costs for a pixel
$\mathrm{p}$ and a disparity $u\in\Gamma=[u_{\mathrm{min}},u_{\mathrm{max}}]$
are recursively aggregated along the path according to:

$\displaystyle
L_{\mathrm{r}}(\mathrm{p},u)=\mathcal{S}(\mathrm{p},u)+\min\limits_{u^{\prime}}\Big{(}L_{\mathrm{r}}(\mathrm{p-r},u^{\prime})+V(u,u^{\prime})\Big{)}.$
(5)

The unary data term of the energy functional, which holds the corresponding
matching costs for pixel $\mathrm{p}$ and disparity $u$ inside the cost volume
$\mathcal{S}$ is denoted as $\mathcal{S}(\mathrm{p},u)$. The smoothness term,
which penalizes deviations in the disparity $u$ of the currently considered
pixel $\mathrm{p}$ and the disparity $u^{\prime}$ of the previous pixel along
the path is formulated by:

$\displaystyle V(u,u^{\prime})$ $\displaystyle=\begin{cases}0&,\ \text{if}\
u=u^{\prime}\\\ \varphi_{1}&,\ \text{if}\ \left|u-u^{\prime}\right|=1\\\
\varphi_{2}&,\ \text{if}\ \left|u-u^{\prime}\right|>1.\end{cases}$ (6)

After all paths for pixel $\mathrm{p}$ have been processed, the aggregated
costs of the individual paths are summed up and stored inside the aggregated
cost volume $\bar{\mathcal{S}}$:

$\displaystyle\bar{\mathcal{S}}(\mathrm{p},u)$
$\displaystyle=\sum\limits_{r}{L_{\mathrm{r}}(\mathrm{p},u)}.$ (7)

From this, the final disparity map is extracted by computing the pixel-wise
WTA solution:

$\displaystyle\mathcal{U}(\mathrm{p})=\operatorname*{arg\,min}\limits_{u}{\bar{\mathcal{S}}(\mathrm{p},u)}.$
(8)

####  2.3.2 Extension of Semi-Global Matching to Plane-Sweep Matching and
Different Surface Orientations

In the following, we describe in detail, how we have adopted the proposed SGM
algorithm to the creation of depth hypotheses with the help of plane-sweep
multi-image matching, as well as the extensions we have employed to account
for non-fronto-parallel surfaces. The proposed extensions SGMx only affect the
aggregation of the matching costs along the concentric paths. Thus, the
extraction of the depth map $\mathcal{D}$ is done analogously to Equation 7
and Equation 8 of the SGM algorithm, with the disparity being substituted by
depth. If a fronto-parallel plane orientation is considered during the plane-
sweep, the depth can be directly extracted from the plane parameterization.
However, for non-fronto-parallel orientations, the depth map $\mathcal{D}$ is
computed by a pixel-wise intersection of the viewing rays with the
corresponding WTA solutions.

Figure 4:  Illustration of the three presented SGMx path aggregations along
one path direction $\mathrm{r}$. Column 1: Reference image and normal map of a
building. Illustrated area marked with yellow line. Column 2: SGMΠ path
aggregation. The blue and pink lines represent the blue and pink surface
orientations on the building facade. Aggregating the path costs for pixel
$\mathrm{p}$ at plane $\Pi$, SGMΠ will incorporate the previous costs at the
same plane position (green) without additional penalty. The previous path
costs at $\Pi\,\pm 1$ (yellow) will be penalized with $\varphi_{1}$. The
previous path costs located at $\Pi\text{+}2$ (red), which is actually located
on the corresponding surface, will be penalized with the highest penalty
$\varphi_{2}$. Column 3: SGM${}^{\Pi\text{-}\mathrm{sn}}$ uses the normal
vector $\mathrm{n_{p}}$, encoding the surface orientation at pixel
$\mathrm{p}$, and computes a discrete index jump $\Delta i_{\mathrm{sn}}$,
which ideally adjusts the zero-cost transition, causing the previous path
costs at $\Pi\text{+}2$ to not be penalized. Column 4: Similar to
SGM${}^{\Pi\text{-}\mathrm{sn}}$, SGM${}^{\Pi\text{-}\mathrm{pg}}$ adjusts the
zero-cost transition. However, the discrete index jump $\Delta
i_{\mathrm{pg}}$ is derived from the running gradient $\nabla\mathrm{r}$ of
the minimum cost path.

##### Resolving Plane Hypotheses with Semi-Global Matching

Since the plane-sweep algorithm does not compute hypotheses on disparities,
but rather pixel-wise plane distances relative to the reference camera and
thus depth, the first SGM extension we propose is a straight-forward adaption
of the standard SGM algorithm to a multi-view plane-sweep sampling. In this,
the formulation of the SGM path aggregation is modified to

$\displaystyle
L_{\mathrm{r}}(\mathrm{p},\Pi)=\mathcal{S}(\mathrm{p},\Pi)+\min\limits_{\delta^{\prime}}\Big{(}L_{\mathrm{r}}(\mathrm{p-r},\Pi^{\prime})+V_{\Pi}(\Pi,\Pi^{\prime})\Big{)},$
(9)

where $\Pi$ denotes the sampling plane at distance $\delta$. The smoothness
term $V_{\Pi}$ now penalizes the selection of different planes between
adjacent pixels along the path $L_{\mathrm{r}}$, instead of disparities. It is
formulated as:

$\displaystyle V_{\Pi}(\Pi,\Pi^{\prime})$ $\displaystyle=\begin{cases}0&,\
\text{if}\ I(\Pi)=I(\Pi^{\prime})\\\ \varphi_{1}&,\ \text{if}\
\left|I(\Pi)-I(\Pi^{\prime})\right|=1\\\ \varphi_{2}&,\ \text{if}\
\left|I(\Pi)-I(\Pi^{\prime})\right|>1,\end{cases}$ (10)

with $I(\cdot)$ being a function that returns the index of $\Pi$ within the
set of sampling planes (cf. Figure 4, Column 2). We denote this extension as
_plane-wise_ SGM (SGMΠ). In our previous publication, (Ruf et al., 2019) we
have referred to this extension as _fronto-parallel_ SGM (SGMfp), since we
have only considered a fronto-parallel sweeping direction so far. However, the
extension is not restricted to a fronto-parallel plane orientation in the
plane-sweep sampling and is also evaluated with slanted planes in the scope of
this work. Given a pixel-wise WTA plane parameterization, the corresponding
depth is extracted by intersecting the viewing ray through pixel
$\mathrm{p}=(p_{x}\ p_{y})^{\intercal}$ with the corresponding plane:

$\displaystyle
d_{\mathrm{p}}=\frac{-\delta}{\mathrm{n}^{\intercal}\cdot\mathbf{K}^{\text{-}1}\cdot(p_{x}\
p_{y}\ 1)^{\intercal}}.$ (11)

##### Incorporating Surface Normals to Adjust the Zero-Cost Transition

The smoothness term of the initial SGM algorithm is formulated with discrete
disparity differences (cf. Equation 6), penalizing discrete disparity jumps
between neighboring pixels. In its optimization scheme, it does not consider
any subpixel disparity and thus favors fronto-parallel surface structures,
leading to stair-casing artifacts if no post-processing is employed
(Scharstein et al., 2017). The same holds for our first extension, SGMΠ. Even
though the plane-sweep sampling also supports non-fronto-parallel plane
orientations, the smoothness term of SGMΠ (cf. Equation 10) does not, and
strongly penalizes index jumps in the sampling planes of more than 1. While
this is desired if the plane-orientation coincides with the surface
orientation, it will still lead to stair-casing artifacts if the surface and
plane orientations do not align. In order to overcome the favoring of fronto-
parallel structures and adjust the smoothness term of SGM to surfaces that are
slanted with respect to the sampling direction, Scharstein et al. (2017)
suggest to add an offset to the smoothness term. This offset can be extracted
from additional information on the surface orientation, e.g. surface normals,
that will make the zero-cost transition coincide with the surface orientation.
We have adopted this approach as part of our second extension and thus call it
_surface normal_ SGM (SGM${}^{\Pi\text{-}\mathrm{sn}}$).

In our hierarchical approach, we extract the normal vectors from the normal
map $\mathcal{N}^{l+1}$, which was estimated in the previous level of the
pyramid (cf. Figure 1). The pixel-wise normal vectors
$\mathrm{n}_{\mathrm{p}}=\mathcal{N}^{l+1}(\mathrm{p})$ indicate the surface
orientation at the scene point $\mathrm{P}$, which is computed by intersecting
the viewing ray through $\mathrm{p}$ with the plane $\Pi$. From this, the
discrete index jump ${\Delta i}_{\mathrm{sn}}$ through the set of sampling
planes can be calculated which is caused by the tangent plane to
$\mathrm{n}_{\mathrm{p}}$. Since the plane-sweep sampling is not restricted to
fronto-parallel plane orientations, the index jump ${\Delta i}_{\mathrm{sn}}$
needs to be computed based on the difference between the tangent plane at
$\mathrm{P}_{\Pi}$ and the orientation of the sampling planes in the direction
$\mathrm{r}$ of the currently considered aggregation path. With ${\Delta
i}_{\mathrm{sn}}$, the smoothness term used by our extension
SGM${}^{\Pi\text{-}\mathrm{sn}}$ is adjusted according to

$\displaystyle
V_{\Pi\text{-}\mathrm{sn}}(\Pi,\Pi^{\prime})=V_{\Pi}(\Pi+{\Delta
i}_{\mathrm{sn}},\Pi^{\prime}).$ (12)

This allows to align the zero-cost transition of the SGM path aggregation to
the surface orientation of the scene (cf. Figure 4, Column 3). The pixel-wise
discrete index jumps can be computed once for each pixel $\mathrm{p}$ and each
path direction $\mathrm{r}$, as also noted by Scharstein et al. (2017),
providing little computational overhead.

##### Penalizing Deviations from the Gradient of the Minimum Cost Path

Instead of relying on additional information, e.g. normal vectors, the third
of our proposed extension computes the running gradient $\nabla\mathrm{r}$ in
scene space, which corresponds to the minimal path costs, in order to adjust
the zero-cost transition in the aggregation of the path costs for non-fronto-
parallel surface orientations. Hence, it is denoted as _path gradient_ SGM
(SGM${}^{\Pi\text{-}\mathrm{pg}}$).

The gradient vector $\nabla\mathrm{r}=\mathrm{P}-\mathrm{P^{\prime}}$ in scene
space is dynamically computed while traversing along the path $\mathrm{r}$. In
this, $\mathrm{P}$ again denotes the scene point that is found by intersecting
the viewing ray through $\mathrm{p}$ with $\Pi$, while $\mathrm{P^{\prime}}$
denotes the scene point, which is parameterized by $\mathrm{p^{\prime}}$ and
the plane $\hat{\Pi}^{\prime}$. Here, $\mathrm{p^{\prime}}=\mathrm{p+r}$
represents the predecessor of $\mathrm{p}$ along the path $\mathrm{r}$ and
$\hat{\Pi}^{\prime}$ denotes the plane at distance
$\hat{\delta}=\operatorname*{arg\,min}_{\delta}{L_{\mathrm{r}}(\mathrm{p^{\prime}},\Pi)}$
associated with the previous minimal path costs.

From this, a discrete index jump ${\Delta i}_{\mathrm{pg}}$ is computed, which
is again used to account for possibly slanted surfaces in scene space by
adjusting the zero-cost transition of the smoothness term according to

$\displaystyle
V_{\Pi\text{-}\mathrm{pg}}(\Pi,\Pi^{\prime})=V_{\Pi}(\Pi+{\Delta
i}_{\mathrm{pg}},\Pi^{\prime}).$ (13)

This implicitly penalizes deviations from the running gradient between two
scene points corresponding to two consecutive pixels on the aggregation path
$\mathrm{r}$ (cf. Figure 4, Column 4).

####  2.3.3 Adaptive Smoothness Penalties Within Semi-Global Matching

In the original publication of the SGM algorithm, Hirschmüller (2005, 2008)
suggests to use an adaptive adjustment of the second penalty $\varphi_{2}$
according to the image gradient along path $\mathrm{r}$. This should enforce a
preservation of depth discontinuities at object boundaries. In this work, the
adaptive adjustment of $\varphi_{2}$ is based on the absolute intensity
difference ($\Delta\mathcal{I}_{\mathrm{pq}}$) between two neighboring pixels
$\mathrm{p}$ and $\mathrm{q}$. It is formulated by

$\displaystyle\varphi_{2}=\varphi_{1}\cdot\left(1+\alpha\cdot\exp\left(-\frac{\Delta\mathcal{I}_{\mathrm{pq}}}{\beta}\right)\right)$
(14)

with
$\Delta\mathcal{I}_{\mathrm{pq}}=|\mathcal{I}(\mathrm{p})-\mathcal{I}(\mathrm{q})|$.
Moreover, $\alpha=8$ and $\beta=10$ are set according to Scharstein et al.
(2017). Since the presented approach uses multiple matching images, the SGM
penalties are multiplied with the number of input images inside the left and
right subsets with respect to $\mathcal{I}_{\mathrm{ref}}$, since the matching
costs are summed up within these image sets.

####  2.3.4 Depth Refinement

In a final optimization step in the process of depth map estimation, the
pixel-wise estimates are refined to lie in between the actual sampling steps
of the plane-sweep sampling. Here, the simple and yet effective approach to
fit a parabola through the WTA solution and its two neighbors is employed.
This is done analogously to implementing a disparity refinement based on curve
fitting for the stereo normal case (Scharstein and Szeliski, 2002). But since
the sampling points are not equidistantly spaced between each other, the depth
differences between the WTA solution and its two neighbors needs to be
considered in the optimization and the finding of the curves’ minimum.

###  2.4 Extraction of Surface Normals from Depth Maps

From the estimated depth map $\mathcal{D}$ our approach computes a normal map
$\mathcal{N}$, holding the local surface orientations in the form of three-
dimensional normal vectors. The surface normal vectors $\mathrm{n_{p}}$ are
calculated for each pixel $\mathrm{p}$ by reprojecting $\mathcal{D}$ into a
three-dimensional point cloud, based on the intrinsic camera parameters and
the corresponding depth estimate. Then, the cross-product
$\mathrm{n_{p}}=\mathrm{h_{p}}\times\mathrm{v_{p}}$ is computed, with
$\mathrm{h_{p}}$ being the difference vector between the scene points of two
neighboring pixels to $\mathrm{p}$ in horizontal direction, and
$\mathrm{v_{p}}$ being the difference vector in vertical direction.

Solely using the cross-product to compute the surface orientation does not
incorporate any local smoothness assumption, which is why we apply an
a-posteriori smoothing to the normal map. In this, an appearance-based
weighted Gaussian smoothing in a local two-dimensional window
$\mathcal{W}_{\mathrm{p}}$ around $\mathrm{p}$ is employed, which adjusts the
smoothing strength depending on the intensity difference between
$\mathrm{q}\in\mathcal{W}_{\mathrm{p}}$ and $\mathrm{p}$:

$\displaystyle\mathcal{N}(\mathrm{p})=\frac{\bar{\mathrm{n}}_{\mathrm{p}}}{\left|\bar{\mathrm{n}}_{\mathrm{p}}\right|},$
(15)

with

$\bar{\mathrm{n}}_{\mathrm{p}}=\mathrm{n_{p}}+\sum\limits_{\mathrm{q}\in\mathcal{W}_{\mathrm{p}}}\mathrm{n_{q}}\cdot\frac{1}{\sqrt{2\pi\sigma^{2}}}\cdot\exp\left(-\frac{\left(\mathrm{q}-\mathrm{p}\right)^{2}}{2\sigma^{2}}-\frac{\Delta\mathcal{I}_{\mathrm{pq}}}{\beta}\right).$
(16)

Similar as in Equation 14, $\beta$ is set to $10$, while $\sigma$ is fixed to
the radius of $\mathcal{W}_{\mathrm{p}}$.

###  2.5 Estimation of Confidence Measures based on Surface Orientation

Apart from the depth map $\mathcal{D}$ and the normal map $\mathcal{N}$, the
presented approach also computes confidence measures with respect to the depth
estimates in the range of $\left[0,\ 1\right]$ and stores them inside a
confidence map $\mathcal{C}$. Such confidence measures allow a subsequent
reasoning on the certainty of the corresponding estimates and, in turn,
improve further processing. Thus, confidence maps are helpful byproducts for
subsequent steps, such as depth map fusion or scene interpretation.
Furthermore, they allow to gain more insight on the effects of different
configurations of the presented approach, based on receiver operating
characteristic (ROC) curve analysis (cf. Figure 8). In our previous work (Ruf
et al., 2019), the computation of the pixel-wise confidence values was based
on the data stored inside the aggregated cost volume $\bar{\mathcal{S}}$
computed by the SGM optimization. In particular, a confidence value based on
the results of the path aggregation, adopted from the observation of Drory et
al. (2014), as well as based on the uniqueness of the WTA solution was used.
Further, experiments and evaluations relying on the analysis of a ROC curve
have shown that these confidence measures do not provide suitable prediction
on the certainty of the depth estimates. Thus, in this paper, the computation
of the pixel-wise confidence measures relies on the geometric characteristic
of the estimated depth map and is deduced from the normal vectors stored
inside the normal map $\mathcal{N}$ and the plane orientations of the plane-
sweep sampling.

In particular, the geometric confidence measure is based on the enclosed
angles between the local surface orientation stored inside the normal map
$\mathrm{n_{p}}=\mathcal{N}(\mathrm{p})$, the orientation of the sampling
plane $\mathrm{n}_{\Pi}$ and the reverted viewing direction $\mathrm{v}$. This
is adopted from the geometric weighting factor proposed by Kolev et al.
(2014). They argue that a depth estimate is more accurate if the surface
orientation of the observed geometry is fronto-parallel to the image plane of
the camera, and less accurate if the camera is observing slanted surfaces.
This correlation is modeled by the scalar product between the surface
orientation and the reverted viewing direction. Furthermore, since the image
warping, as part of the image matching, can be aligned to the surface
orientation by adjusting the normal vector of the plane-sweep algorithm, the
plane orientation is also considered. Thus, the geometry-based weighting
factor is computed according to:

$\displaystyle\mathcal{C}(\mathrm{p})$
$\displaystyle=\begin{cases}\frac{\langle\mathrm{n}_{\mathrm{p}},\mathrm{n}_{\Pi}\rangle\langle\mathrm{n}_{\Pi},\mathrm{v}\rangle\
-\ \cos\rho}{1\ -\ \cos\rho}&,\ \text{if}\
\\{\sphericalangle(\mathrm{n}_{\mathrm{p}},\mathrm{n}_{\Pi})\land\sphericalangle(\mathrm{n}_{\Pi},\mathrm{v})\\}\leq\rho\\\
0&,\ \text{otherwise}\ .\end{cases}$ (17)

All of the above vectors are assumed to be normalized and given with respect
to the local coordinate system of the camera, thus $\mathrm{v}=(0\ 0\
\text{-}1)^{\intercal}$. Just as in the work of Kolev et al. (2014), a
critical angle $\rho=60^{\circ}$ is used to mark the measurements, for which
the enclosed angles exceed this threshold, as unreliable. The additional
consideration of $\mathrm{n}_{\Pi}$ in Equation 17 implicitly models the
indirect matching of the input images via the plane-induced homography.

###  2.6 Post-Processing and Depth Map Filtering

In a final post-processing step, remaining outliers are removed from the
depth, normal and confidence maps by applying a Difference-of-Gaussian (DoG)
filtering (Section 2.6.1) and a geometric consistency check (Section 2.6.2).

####  2.6.1 Difference-of-Gaussian Filtering

Data: unfiltered depth, normal and confidence maps ($\mathcal{D}$,
$\mathcal{N}$ and $\mathcal{C}$) as well as corresponding reference image
$\mathcal{I}_{\mathrm{ref}}$.

Result: filtered $\mathcal{D}$, $\mathcal{N}$ and $\mathcal{C}$, in which all
estimates corresponding to weakly textured areas in
$\mathcal{I}_{\mathrm{ref}}$ are removed.

1 use a Gaussian filter with a kernel of $7\times 7$ pixels to smooth the
reference frame $\mathcal{I}{{}_{\mathrm{ref}}}$, yielding
$\mathcal{I}^{\mathrm{smooth}}_{\mathrm{ref}}$.

2compute the DoG image depicting local image gradients, according to:
$\mathcal{I}^{\mathrm{DoG}}_{\mathrm{ref}}=\mathcal{I}_{\mathrm{ref}}-\mathcal{I}^{\mathrm{smooth}}_{\mathrm{ref}}$.

3apply a binary threshold to compute the DoG mask
$\mathcal{M}^{\mathrm{DoG}}$, marking all image areas in which the intensity
change is greater than $0.5$.

4remove activation areas smaller than $7$ pixels in
$\mathcal{M}^{\mathrm{DoG}}$ by applying a speckle filter.

5dilate $\mathcal{M}^{\mathrm{DoG}}$ with a kernel size of $3\times 3$ pixels
to fill small holes in activation areas.

6remove deactivation areas smaller than $21$ pixels by applying a speckle
filter to the inverted DoG mask
$\mathcal{M}^{\mathrm{DoG}\text{-}\mathrm{inv}}=1-\mathcal{M}^{\mathrm{DoG}}$.

7invalidate pixels in $\mathcal{D}$, $\mathcal{N}$ and $\mathcal{C}$ for which
$\mathcal{M}^{\mathrm{DoG}}=1$.

Algorithm 3 Overview on the Difference-of-Gaussian filter to invalidate all
image areas belonging to weakly textured areas.

As proposed by Wenzel (2016), the DoG filter allows to remove estimates from
$\mathcal{D}$, $\mathcal{N}$ and $\mathcal{C}$ by masking out pixels in image
regions that only provide little textural information (e.g. unsharp or
overexposed areas). Here, it is assumed that in such regions the image
matching is ambiguous and that it leads to less accurate results. The DoG
filter is used to detect weakly textured areas inside the reference image
$\mathcal{I}_{\mathrm{ref}}$ and build up a binary image mask, which, in turn,
is used to remove the estimates from the corresponding maps. Algorithm 3
provides an overview on the implementation of the employed DoG filter, which
is similar to the one proposed in (Wenzel, 2016).

####  2.6.2 Geometric Consistency based on Mutual Reprojection Error

In the stereo normal case, in which the depth is constructed from only two
views, a simple and yet very effective approach to enforce geometric
consistency and perform occlusion detection is the so-called left-right
consistency check. In this, the pixel-wise disparities stored inside the
disparity image of the reference image are compared to the approximated
disparities inside the disparity map of the matching image. If they differ by
a certain threshold, the corresponding estimates inside the disparity map of
the reference image are invalidated. This can also be formulated as the mutual
reprojection error, in which each pixel $\mathrm{p}^{\mathrm{ref}}$ of
$\mathcal{D}_{\mathrm{ref}}$, having a depth estimate
$d^{\,\mathrm{ref}}_{\mathrm{p}}$, is projected into the view of a second
depth map $\mathcal{D}^{k}$, according to $d^{\,\mathrm{ref}}_{\mathrm{p}}$
and the corresponding projection matrices $\mathbf{P}_{\mathrm{ref}}$ and
$\mathbf{P}_{k}$, yielding the image point $\mathrm{p}^{k}$. Given
$\mathrm{p}^{k}$ and the corresponding depth $d^{\,k}_{\mathrm{p}}$ from
$\mathcal{D}_{k}$, the image point $\mathrm{p}^{k}$ is projected back into the
view of $\mathcal{D}_{\mathrm{ref}}$, resulting in
$\tilde{\mathrm{p}}^{\mathrm{ref}}$. Finally, if the Euclidean distance
between $\mathrm{p}^{\mathrm{ref}}$ and $\tilde{\mathrm{p}}^{\mathrm{ref}}$
exceeds a given threshold $\eta_{\mathrm{r}}$ the estimate at
$\mathrm{p}^{\mathrm{ref}}$ is invalidated. Schönberger et al. (2016)
formulate this reprojection error for pixel $\mathrm{p}$ in a reference view
and a neighboring map $k$ as
$\epsilon^{k}_{\mathrm{r}}(\mathrm{p})=\|\mathrm{p}-\mathbf{H}^{k}_{\mathrm{p}}\cdot\mathbf{H}_{\mathrm{p}}\cdot\mathrm{p}\|$,
with $\mathbf{H}_{\mathrm{p}}$ being the forward projection into the view $k$
according to $d^{\,\mathrm{ref}}_{\mathrm{p}}$, and
$\mathbf{H}^{k}_{\mathrm{p}}$ being the corresponding backward projection
according to $d^{\,k}_{\mathrm{p}}$.

We have adopted this approach to perform a final geometry-based filtering
between a number of depth maps within a sliding window. This is not part of
the actual hierarchical processing pipeline, but rather a separate post-
processing step, since it requires the results of other image bundles of the
input sequence. If possible, the center-most depth map of the sliding window
$\Psi$ is chosen as a reference view for which the filtering is performed. At
the beginning or the end of the sequence, where the sliding window would
exceed the boundaries, the window is shifted to either side of the reference
view so that it is always inside the boundaries of the sequence and that no
depth map is filtered multiple times. Apart from thresholding the reprojection
error $\epsilon^{k}_{\mathrm{r}}$, another criterion is introduced to evaluate
the geometric consistency, namely the number of neighboring views for which
the reprojection error is within the threshold, i.e. the number of hits:
$\epsilon_{\textrm{h}}(\mathrm{p})=\sum_{k}[\epsilon^{k}_{\mathrm{r}}(\mathrm{p})<\eta_{\mathrm{r}}]$,
with $[\cdot]$ being the Iverson bracket. Algorithm 4 gives an overview on
geometric consistency and the according filtering of the depth, normal and
confidence maps. In the scope of this work, the size of the sliding window is
set to $|\Psi|=5$, the threshold for the reprojection error to
$\eta_{\mathrm{r}}=10$ and the consistency threshold to $\eta_{\mathrm{h}}=3$.

Data: depth, normal and confidence maps ($\mathcal{D}_{k}$, $\mathcal{N}_{k}$
and $\mathcal{C}_{k}$) within a sliding window $\Psi$ of the input sequence as
well as corresponding projection matrices $\mathbf{P}_{k}$.

Result: filtered $\mathcal{D}_{\mathrm{ref}}$, $\mathcal{N}_{\mathrm{ref}}$
and $\mathcal{C}_{\mathrm{ref}}$ of reference view, in which all estimates
that are not geometrically consistent are removed.

1select $\mathcal{D}_{\mathrm{ref}}$, $\mathcal{N}_{\mathrm{ref}}$ and
$\mathcal{C}_{\mathrm{ref}}$ corresponding to the center-most view within the
sliding window $\Psi$.

2foreach  _pixel $\mathrm{p}^{\mathrm{ref}}\in\mathcal{D}_{\mathrm{ref}}$ and
neighboring view $k\in\Psi$_ do

3

4 calculate number of hits for which reprojection error is below threshold:
$\epsilon_{\textrm{h}}(\mathrm{p})=\sum\limits_{k}[\epsilon^{k}_{\mathrm{r}}(\mathrm{p})<\eta_{\mathrm{r}}],\
\text{with}\ \
\epsilon^{k}_{\mathrm{r}}(\mathrm{p})=\|\mathrm{p}-\mathbf{H}^{k}_{\mathrm{p}}\cdot\mathbf{H}_{\mathrm{p}}\cdot\mathrm{p}\|.$

5 if $\epsilon_{\textrm{h}}<\eta_{\mathrm{h}}$, invalidate pixel $\mathrm{p}$
in $\mathcal{D}_{\mathrm{ref}}$, $\mathcal{N}_{\mathrm{ref}}$ and
$\mathcal{C}_{\mathrm{ref}}$, by setting it to $0$.

6 end foreach

Algorithm 4 Overview on the geometric consistency filter.

##  3 Results

The following sections present the results of experiments conducted in the
scope of this work. In these, the effects of different configurations of the
presented approach are evaluated and analyzed with respect to individual
aspects, such as accuracy, efficiency and application-specific usability.
First, the datasets and evaluation metrics used to investigate the potential
of the presented approach are introduced in Section 3.1 and Section 3.2,
respectively. Then, the need for a hierarchical processing is evaluated and
the optimal number of input images, i.e. the optimal size of the input bundle,
is empirically found in Section 3.3. This is followed by a short comparison in
Section 3.4, where the focus is set on the two similarity metrics and cost
functions that are implemented in the scope of this work. In Section 3.5, the
ability of the three SGM extensions to reconstruct non-fronto-parallel surface
structures is evaluated and compared to the effects of using a non-fronto-
parallel plane orientation within plane-sweep sampling. An evaluation of the
improvements gained by post-filtering is presented in Section 3.6, before
comparing the results of the best configurations to those achieved by
approaches for offline MVS, e.g. COLMAP (Schönberger et al., 2016), in Section
3.7. Finally, before discussing the presented results, the outcomes of use-
case-specific experiments are shown and qualitatively illustrated in Section
3.8.

The complete processing pipeline of the presented approach, except the
generation of the Gaussian image pyramids and the parameterization of the
plane-sweep algorithm, is implemented in CUDA and is thus optimized for
massively parallel computing by general purpose computation on a GPU (GPGPU),
which, in turn, is embedded in a C++ application. All experiments, and thus
all timing measurements, were conducted using a NVIDIA Titan X GPU and an
Intel XEON CPU E5-2650 running with 2.20 GHz. Even though the CPU is designed
for a server architecture, only a small part of our approach is run on the
CPU, and thus its superiority over commodity desktop hardware is
insignificant.

###  3.1 Evaluation Datasets

The presented approach is quantitatively evaluated on two public datasets,
which also provide an appropriate ground truth, namely the DTU Robot MVS
dataset (Jensen et al., 2014; Aanæs et al., 2016) and the dataset from the
3DOMcity Benchmark (Özdemir et al., 2019). In order to provide appropriate
data for a quantitative evaluation, these two datasets rely on images of
scale-modeled buildings and an urban scenery from which an accurate ground
truth is acquired. For a qualitative evaluation and a discussion regarding the
usability of the presented approach for online dense image matching and 3D
reconstruction, two privately captured datasets of real-world scenes are used,
henceforth referred to as the TMB and the FB dataset. In the following, the
characteristics of these datasets are briefly introduced. In particular, which
portions of the data sets are used and what kind of ground truth is available
for the evaluation.

Figure 5:  Overview of the four datasets used for performance evaluation in
the scope of this paper. Column 1: Two building models from the DTU Robot MVS
dataset. Column 2: Example images in oblique and nadir view from the 3DOMcity
Benchmark dataset. Column 3: Excerpt of the privately acquired TMB dataset.
Column 4: Use-case-specific dataset acquired during an exercise of the local
fire brigade.

##### DTU Robot MVS Dataset

The DTU Robot MVS dataset (Figure 5, Column 1) is comprised of 124 different
table top scenes. For each scene, input images are provided, which were
captured under eight different lighting conditions from 49 locations. These
locations are the same for all scenes, since the camera is mounted on an
industrial robot arm, which could repeatedly be moved to the set camera poses.
Furthermore, for each scene, a ground truth is provided in the form of a
detailed point cloud captured by a structured light scanner.

For the quantitative evaluation of the performance of the presented approach,
21 scans of different building models were selected, since these scenes are
closest to the targeted use-case, namely the online reconstruction of urban
scenes and man-made structures. In this, the already undistorted images with
an image resolution of $1600\times 1200$ pixels, together with the
corresponding intrinsic and extrinsic camera data, were used as input data to
the approach. The benchmark provides an evaluation routine, and a
corresponding score table, for the assessment of a reconstructed point cloud
with respect to the ground truth captured by the structured light scanner.
However, the focus of this work lies in the estimation of depth and normal
maps only, without the subsequent fusion into a three-dimensional point cloud.
Thus, in order to evaluate the accuracy of the estimated depth and normal
maps, corresponding ground truth data is extracted from the structured light
scans by rendering them from the view points of the corresponding input
images, given the intrinsic and extrinsic camera data.

##### 3DOMcity Benchmark Dataset

Within the DTU benchmark, the camera is moved along a circular trajectory,
resembling an orbital flight of a UAV, with the camera focusing on the object
of interest. This kind of camera movement or flight trajectory is typical for
the case when an object is to be fully reconstructed with a maximum precision
(Wenzel et al., 2013a). However, depending on the aircraft and the surrounding
constraints, such a flight is not always feasible or desired. Typical for the
mapping of a larger area, a grid flight is performed, where the aircraft is
flying linearly over the area of interest, with the camera orientation being
fixed with respect to the sensor carrier.

Such a configuration is simulated by the data provided as part of the 3DOMcity
Benchmark (Özdemir et al., 2019) (Figure 5, Column 2). In this, images of a
scale-modeled urban scene, comprised of differently sized and shaped
buildings, as well as streets and vegetation, are captured with a DSLR camera
moved along a rigid bar in parallel lines over the model. The distorted and
undistorted images are provided with a maximum resolution of $6016\times 4016$
pixels from oblique and nadir vantage points, with a forward and sideways
overlap of $85/75\,\%$ in case of the oblique views and $80/65\,\%$ between
the nadir views. While the benchmark internally uses a point cloud of two
exemplary building models, captured by a multi-stripe triangulation-based
laser scanner, to assess the accuracy of submitted DIM point clouds, it
publicly provides a reference point cloud for the task of classification,
which is computed by a semi-global DIM algorithm (Özdemir et al., 2019).

To use the data of the 3DOMcity Benchmark for the quantitative evaluation of
the performance of the presented approach, the already undistorted images are
first down-scaled to a size of $1798\times 1200$ pixels, preserving the
initial aspect ratio, before estimating the intrinsic camera parameters with
the help of COLMAP (Schönberger and Frahm, 2016). The extrinsic camera data is
extracted from the reference that is provided as part of the benchmark. For
the accuracy assessment of the depth maps, the reference data is computed by
rendering the reference DIM point cloud of the whole model, provided for the
task of point cloud classification, from the viewpoints of the input images,
just as in the case of the dataset of the DTU benchmark.

##### Real-World Use-Case-Specific TMB and FB Dataset

The strength and aim of the datasets from the DTU and 3DOMcity benchmarks are
their small scale and the associated ability to record or compute accurate
reference data, which in turn facilitates a quantitative evaluation of the
accuracy of the assessed algorithms. These datasets, however, are recorded
under controlled environments and do not fully accommodate for the use-case
aimed at with the presented approach, namely the online MVS from monocular
video data captured by commodity UAVs, flying at altitudes below $100$ m. In
order to account for the named use-case and perform a qualitative evaluation
on real-world data, appropriate test data was collected into a private
dataset.

This dataset is two-fold. The first part, namely the TMB dataset (Figure 5,
Column 3), consists of four sequences captured by a DJI Phantom $3$
professional, flying around a free-standing house and containers at altitudes
between $8$ m to $15$ m. For the second part, which is denoted as the fire
brigade (FB) dataset (Figure 5, Column 4), images were acquired during a fire
brigade exercise around a big industrial building. The data was captured using
a DJI Matrice 200 with a Zenmuse XT2 sensor flying linearly over the area on
which the exercise was performed. For all sequences, the images were captured
at a frame-rate of about $1$ FPS and down-sampled to an image size of
$1920\times 1080$ pixels. However, due to the varying velocities, the
distances between the frames are not always the same and thus images that are
not appropriate as input to the presented approach, e.g. by providing too
little offset, are discarded. As reference data, detailed point clouds of each
sequence were computed by means of Structure-from-Motion (SfM) and MVS using
COLMAP (Schönberger and Frahm, 2016; Schönberger et al., 2016). In order to
have a metric reference, given the GNSS metadata provided by the UAVs for each
input image, the reference data was transformed into a local east-north-up
(ENU) frame. The undistorted images as well as the intrinsic and extrinsic
camera data produced by COLMAP serve as input to the presented approach for
the evaluation with respect to the TMB and FB dataset.

###  3.2 Error Measures

For a direct quantification of the error between the estimated depth map
$\mathcal{D}_{\mathrm{est}}$ and the corresponding ground truth
$\mathcal{D}_{\mathrm{gt}}$, absolute and relative L1 measures are used:

$\displaystyle\text{L1-abs}(\mathcal{D}_{\mathrm{est}},\mathcal{D}_{\mathrm{gt}})=\frac{1}{|\mathcal{V}|}\sum\limits_{\mathrm{p}\in\mathcal{V}}{|\mathcal{D}_{\mathrm{est}}(\mathrm{p})-\mathcal{D}_{\mathrm{gt}}(\mathrm{p})|}\
,\ \ \text{and}$ (18)
$\displaystyle\text{L1-rel}(\mathcal{D}_{\mathrm{est}},\mathcal{D}_{\mathrm{gt}})=\frac{1}{|\mathcal{V}|}\sum\limits_{\mathrm{p}\in\mathcal{V}}{\frac{|\mathcal{D}_{\mathrm{est}}(\mathrm{p})-\mathcal{D}_{\mathrm{gt}}(\mathrm{p})|}{\mathcal{D}_{\mathrm{gt}}(\mathrm{p})}}\
.$ (19)

In this, $\mathcal{V}$ denotes the set of pixels for which both
$\mathcal{D}_{\mathrm{est}}$ and $\mathcal{D}_{\mathrm{gt}}$ have valid depth
measurements. While L1-abs provides an absolute and, in turn, interpretable
insight on the mean error of the estimated depth map, it is rather unsuitable
for comparing the results across multiple datasets with different depth
ranges. This is because the error of depth measurements typically increases
with increasing depth, which leads to a higher absolute error for datasets
with a larger scene depth. In order to compensate for this effect, the
relative L1-rel measure normalizes the absolute difference by the depth stored
at the corresponding ground truth pixel. This reduces the effect that
erroneous pixels in distant areas of the scene have on the error score, while
at the same time increasing the weight of the pixels that are close to the
camera.

The two error measures introduced above provide a simple strategy to evaluate
the error of the estimates. However, they do not allow to reason about the
completeness and density of the estimated depth map. Thus, since the focus of
this work is on dense MVS, it is also of great interest to know how many
pixels of $\mathcal{D}_{\mathrm{est}}$ are actually filled with correct
estimates. To do so, two tightly coupled error scores, namely the accuracy
($\text{Acc}_{\theta}$) and completeness ($\text{Cpl}_{\theta}$), are also
used. These scores are typically used to evaluate classification tasks, but
have also been used to evaluate range measurements in recent years (Schöps et
al., 2017; Knapitsch et al., 2017). On the one hand, the accuracy
$\text{Acc}_{\theta}$ indicates the amount of pixels within the estimated
depth map $\mathcal{D}_{\mathrm{est}}$, for which the corresponding depth
value is within a given threshold $\theta$ to the ground truth:

$\text{Acc}_{\theta}(\mathcal{D}_{\mathrm{est}},\mathcal{D}_{\mathrm{gt}})=\frac{1}{|\mathcal{E}|}\sum\limits_{\mathrm{p}\in\mathcal{V}}\left[\max\left(\frac{\mathcal{D}_{\mathrm{est}}(\mathrm{p})}{\mathcal{D}_{\mathrm{gt}}(\mathrm{p})},\frac{\mathcal{D}_{\mathrm{gt}}(\mathrm{p})}{\mathcal{D}_{\mathrm{est}}(\mathrm{p})}\right)<\theta\right]\
.$ (20)

The completeness $\text{Cpl}_{\theta}$, on the other hand, indicates the
fraction of the ground truth pixels, for which estimates exist, which are
within the given distance threshold to the reference:

$\text{Cpl}_{\theta}(\mathcal{D}_{\mathrm{est}},\mathcal{D}_{\mathrm{gt}})=\frac{1}{|\mathcal{G}|}\sum\limits_{\mathrm{p}\in\mathcal{V}}\left[\max\left(\frac{\mathcal{D}_{\mathrm{est}}(\mathrm{p})}{\mathcal{D}_{\mathrm{gt}}(\mathrm{p})},\frac{\mathcal{D}_{\mathrm{gt}}(\mathrm{p})}{\mathcal{D}_{\mathrm{est}}(\mathrm{p})}\right)<\theta\right]\
.$ (21)

Again, $\mathcal{V}$ holds the set of pixels for which both
$\mathcal{D}_{\mathrm{est}}$ and $\mathcal{D}_{\mathrm{gt}}$ have valid depth
measurements. Similarly, $\mathcal{E}$ denotes the pixel set with valid
estimates, while $\mathcal{G}$ holds the pixels with valid ground truth
values. In both Equation 20 and Equation 21, the operator $[\cdot]$ refers to
the Iverson bracket. The threshold $\theta$ is given as the percentage of the
corresponding ground truth value. For example, $\text{Acc}_{1.25}$ and
$\text{Cpl}_{1.25}$ hold the fraction of pixels with respect to the
$\mathcal{D}_{\mathrm{est}}$ and $\mathcal{D}_{\mathrm{gt}}$, for which the
difference between the estimate and ground truth is smaller than $25$ % of the
corresponding ground truth depth. These two measures can further be summarized
into a combined score, namely the $\text{F}_{\theta}$-score, which is the
harmonic mean between the $\text{Acc}_{\theta}$ and $\text{Cpl}_{\theta}$:

$\text{F}_{\theta}(\mathcal{D}_{\mathrm{est}},\mathcal{D}_{\mathrm{gt}})=2\cdot\frac{\text{Acc}_{\theta}\cdot\text{Cpl}_{\theta}}{\text{Acc}_{\theta}+\text{Cpl}_{\theta}}\
.$ (22)

Thus, a high $\text{F}_{\theta}$-score indicates a good trade-off between the
achieved accuracy of the depth map and its completeness with respect to the
ground truth.

###  3.3 Need for Hierarchical Processing and Finding the Optimal Number of
Input Images

In the first experiment, the need and importance of the hierarchical
processing scheme within the presented approach as well as the best
configuration on the size of the input bundle are evaluated. In this, a couple
of aspects are considered in order to find the best configuration for the
succeeding experiments. The objective is to find the appropriate number $n$ of
Gaussian pyramid levels and the size $|\Omega|$ of the input bundle, providing
a good trade-off between

  * 1.

the error of the resulting depth maps, measured by L1-abs and L1-rel,

  * 2.

the sampling density of the scene and the entailed resources needed for the
computation,

  * 3.

as well as the resulting processing run-time.

For this experiment, a fronto-parallel plane orientation is used as part of
the plane-sweep image matching and the NCC with a support region of $5\times
5$ pixels is set as similarity measure. The optimization of the cost volume
and the extraction of the optimal depth map is done by employing the SGMΠ
scheme, which is the adoption of the standard SGM optimization to the use of
plane-sweep image matching (cf. Section 2.3.2). The smoothness penalty within
the SGM optimization is set to $\varphi_{1}=100$, while the adaptive
$\varphi_{2}$ penalty as described in Section 2.3.3 is used. This, together
with the $5\times 5$ sized NCC as matching cost, was chosen in accordance with
the work of Scharstein et al. (2017). To find the appropriate height of the
Gaussian pyramids, the size of the input bundle, i.e. the number of input
images, is first set to $|\Omega|=3$.

Table 1: Mean L1 errors achieved on the DTU and 3DOMcity dataset for a different number $n$ of Gaussian pyramid levels as part of the hierarchical processing scheme. The error metrics used are the absolute L1-abs, measured in mm, as well as the relative L1-rel measure. Both are averaged over all evaluated depth maps within each dataset. |  | $\boldsymbol{n=1}$ | $\boldsymbol{n=2}$ | $\boldsymbol{n=3}$ | $\boldsymbol{n=4}$ | $\boldsymbol{n=5}$
---|---|---|---|---|---|---
DTU | L1-abs | 26.394 | 26.221 | 23.473 | 25.045 | 29.676
(in mm) |  $\pm$24.262 |  $\pm$23.835 |  $\pm$19.656 |  $\pm$19.298 |  $\pm$19.436
L1-rel | 0.036 | 0.036 | 0.032 | 0.034 | 0.041
|  $\pm$0.032 |  $\pm$0.032 |  $\pm$0.026 |  $\pm$0.026 |  $\pm$0.026
3DOMcity | L1-abs | 12.789 | 14.936 | 21.801 | 32.458 | 47.422
(in mm) |  $\pm$6.916 |  $\pm$6.754 |  $\pm$8.010 |  $\pm$9.408 |  $\pm$22.292
L1-rel | 0.010 | 0.012 | 0.017 | 0.026 | 0.037
|  $\pm$0.006 |  $\pm$0.006 |  $\pm$0.007 |  $\pm$0.009 |  $\pm$0.014

Table 1 lists the mean errors of the estimated depth maps when evaluated with
different number of pyramid levels on the datasets of both the DTU and
3DOMcity benchmark. In this, the absolute and relative L1 measures are used,
averaged over all depth maps within each dataset. It is to be expected, that
the omission of any hierarchical processing, i.e. the use of only one pyramid
level and thus no coarse-to-fine processing, would lead to the smallest error
between the estimate and ground truth. However, the results reveal that in
case of the DTU dataset the smallest mean error, even if it is only slightly
smaller, is achieved when setting $n=3$, while the best result in case of the
3DOMcity dataset is achieved at $n=1$.

As described in Section 2.2.3, the plane distances within the plane-sweep
sampling, and thus the sampling points, are selected in such a way that two
consecutive planes induce a maximum disparity difference of $1$ pixel.
Depending on the capturing setup, i.e. the relative poses between the images
and their obliqueness and, in turn, the range of the scene depth, this can
lead to a very high number of sampling points and with it to a large memory
consumption, as the dimensions of the three-dimensional cost volume need to be
set accordingly. Thus, in order to not exceed the memory limit, the maximum
number of sampling points for the highest pyramid level is restricted to $256$
in the implementation of the approach. In case of the camera setup of the DTU
dataset and the configuration of this experiment, i.e. having a bundle size of
$|\Omega|=3$, a pyramid height of $3$ is the smallest height at which the
number of sampling points at the highest level does not reach or exceed the
set limit, as Table 2 shows. Comparing Table 1 and Table 2 further reveals
that on both datasets the best results are achieved, when the computation is
initialized at the highest pyramid level with a maximum of $128$ sampling
planes.

Table 2: Processing run-time measured for different configurations of the pyramid height on the DTU and 3DOMcity dataset. In addition, the maximum number of sampling planes with which the scene was sampled at the highest pyramid level is stated. |  | $\boldsymbol{n=1}$ | $\boldsymbol{n=2}$ | $\boldsymbol{n=3}$ | $\boldsymbol{n=4}$ | $\boldsymbol{n=5}$
---|---|---|---|---|---|---
DTU | Run-time | 2365 | 1315 | 386 | 220 | 187
(in ms) |  $\pm$15 |  $\pm$10 |  $\pm$2 |  $\pm$2 |  $\pm$1
max. # planes | 256 | 256 | 128 | 64 | 32
3DOMcity | Run-time | 613 | 431 | 225 | 196 | 192
(in ms) |  $\pm$3 |  $\pm$3 |  $\pm$1 |  $\pm$1 |  $\pm$1
max. # planes | 128 | 64 | 32 | 16 | 8

Another criterion which is used to deduce the best configuration on the height
of the Gaussian pyramid is the run-time needed to estimate a single depth map.
Table 2 additionally lists the corresponding measurements taken, i.e. the
number of milliseconds it takes to estimate a single depth map given a certain
number of pyramid levels, as well as the number of planes used for sampling
the scene space at the highest pyramid level. The measurements again show,
that up to $n=3$ in case of the DTU dataset, the number of sampling planes at
the highest pyramid level is equal to the limit of 256 and that with a smaller
amount of sampling points the run-time decreased drastically. Furthermore, the
significant drop of one second in run-time between using a pyramid height of 2
and 3 suggests that the decreasing use of processing resources on the GPU
increases the processing speed and that going from $n=2$ to $n=3$ makes a
significant improvement in its efficiency. Since the use of a higher number of
pyramid levels does not only reduce the amount of sampling points, but also
the image size at the highest pyramid level and with it the amount of pixels
that need to be matched. Thus, depending on the camera setup, a hierarchical
processing is very important in order to ensure a high sampling density of the
scene space, while at the same time efficiently utilizing the processing
hardware and, in turn, alleviating high processing speeds. In case of the DTU
dataset, this experiment shows that the best number of pyramid levels to be
used is $n=3$, which will thus be set for the successive experiments. In case
of the 3DOMcity dataset, Table 1 suggests that the best configuration is to
use the original image size. A hierarchical processing scheme is needed,
however, in order to use SGM${}^{\Pi\text{-}\mathrm{sn}}$, the extension of
the SGM algorithm to consider local surface orientations in order to account
for slanted surfaces. Thus, in case of the 3DOMcity dataset, the successive
experiments will be executed with $n=2$, which induces only a slightly higher
mean error compared to the best configuration.

Table 3: Mean errors achieved on the DTU and $3$DOMcity dataset for different input bundle sizes $|\Omega|$, i.e. number of images. In addition, the differences in run-time, with respect to the measurements of the first part (i.e. $|\Omega|=3$), are stated. |  | $\boldsymbol{|\Omega|=3}$ | $\boldsymbol{|\Omega|=5}$ | $\boldsymbol{|\Omega|=7}$
---|---|---|---|---
DTU | L1-abs | 23.473 | 19.832 | 21.843
(in mm) |  $\pm$19.656 |  $\pm$16.225 |  $\pm$21.605
L1-rel | 0.032 | 0.027 | 0.031
(in mm) |  $\pm$0.026 |  $\pm$0.021 |  $\pm$0.031
$\Delta$ Run-time |  | +271 | +302
(in ms) |  |  |
3DOMcity | L1-abs | 14.936 | 14.615 | 16.514
(in mm) |  $\pm$6.754 |  $\pm$6.254 |  $\pm$7.569
L1-rel | 0.012 | 0.012 | 0.014
(in mm) |  $\pm$0.006 |  $\pm$0.007 |  $\pm$0.009
$\Delta$ Run-time |  | +360 | +410
(in ms) |  |  |

In the second part of this experiment, the effects of a different number of
input images and, in turn, the optimal size $|\Omega|$ of the input bundle are
evaluated. Here, the settings for the plane-sweep image matching and the
subsequent SGM optimization are kept the same as before. The height of the
Gaussian pyramids is fixed to $n=3$ in case of the DTU dataset and $n=2$ in
case of the data from 3DOMcity dataset. Table 3 lists the mean errors achieved
on both datasets with different number of input images, as well as the
difference in run-time with respect to the best configuration of the first
part of the experiment. The results reveal that the best accuracies are
achieved, when five input images are used for image matching, even though, in
case of the 3DOMcity dataset, it is only a marginal improvement. As expected,
the utilization of more input images in the process of image matching also
leads to an increase in run-time, since more pixels are matched. At the same
time, however, there is more time available to keep up with the image
acquisition as discussed in Section 4.3. In conclusion, in the subsequent
experiments, the size of the input bundle is set to $|\Omega|=5$, while the
height of the Gaussian image pyramids is set to $n=3$ and $n=2$ in case of the
DTU and 3DOMcity dataset, respectively.

###  3.4 Effects of Different Similarity Measures in the Process of Dense
Multi-Image Matching

As part of the plane-sweep multi-image matching, this approach comprises two
different similarity measures and cost functions: The Hamming distance of the
census transform (CT) as well as a truncated and scaled form of the normalized
cross-correlation (NCC). While the CT is computationally less expensive than
the NCC and is thus more suitable for real-time or online processing, it is
less discriminative, which might result in a more ambiguous set of matched
pixel correspondences. When working with a stereo normal case, in which the
input images suffer only from little perspective distortion due to homographic
transformations, the CT outperforms the NCC in both run-time and accuracy (Ruf
et al., 2021). However, as the results in Table 4 show, the perspective
distortion, resulting from the warping of images from converging cameras by
means of the plane-induced homography within the plane-sweep algorithm, leads
to a significant increase in error when using the CT as similarity measure
instead of the NCC.

Table 4: Mean errors achieved on the DTU and 3DOMcity dataset when using different similarity measures and cost functions with different support regions. |  | $\boldsymbol{\text{CT}_{5\times 5}}$ | $\boldsymbol{\text{CT}_{9\times 7}}$ | $\boldsymbol{\text{NCC}_{5\times 5}}$ | $\boldsymbol{\text{NCC}_{9\times 9}}$
---|---|---|---|---|---
DTU | L1-abs | 42.136 | 42.305 | 19.832 | 19.667
(in mm) |  $\pm$37.958 |  $\pm$36.394 |  $\pm$16.225 |  $\pm$16.453
L1-rel | 0.056 | 0.057 | 0.027 | 0.027
(in mm) |  $\pm$0.048 |  $\pm$0.046 |  $\pm$0.021 |  $\pm$0.021
3DOMcity | L1-abs | 22.128 | 26.005 | 14.615 | 13.789
(in mm) |  $\pm$14.218 |  $\pm$14.106 |  $\pm$6.254 |  $\pm$5.962
L1-rel | 0.019 | 0.022 | 0.012 | 0.011
(in mm) |  $\pm$0.014 |  $\pm$0.014 |  $\pm$0.007 |  $\pm$0.006

Apart from the two different similarity measures, the effects of different
support regions are also evaluated in the scope of this experiment. In this,
for each similarity measure, the two most commonly used configurations were
tested, with a support region of a size of $5\times 5$ pixels being a good
trade-off between uniqueness and computational complexity, while, in case of
the CT, a support region of a size of $9\times 7$ pixels is the biggest size
for which the bit-string still fits into a single $64$-bit integer. The
configuration of the plane-sweep algorithm and the SGM optimization is set in
accordance with the values from the first experiment (cf. Section 3.3). In
terms of the SGM penalties, $\varphi_{1}$ is set to $100$ for both
$\text{NCC}_{5\times 5}$ and $\text{NCC}_{9\times 9}$, since the maximum
matching cost of the NCC is normalized to $255$, independent of the support
region. For $\text{CT}_{5\times 5}$ and $\text{CT}_{9\times 7}$, however,
$\varphi_{1}$ is set to $9$ and $24$, respectively, which is equivalent to the
configuration for NCC, when considering the ratio between $\varphi_{1}$ and
the maximum matching cost.

Even though the NCC with a support region of $9\times 9$ pixels achieves the
best results, $\text{NCC}_{5\times 5}$ is selected for further experiments,
since the rise in error is only little but its computational complexity is
less and, in turn, its throughput higher than that of $\text{NCC}_{9\times 9}$
as measured by Ruf et al. (2021).

###  3.5 Ability to Reconstruct Non-Fronto-Parallel Surface Structures

Table 5: Quantitative comparison of the results achieved by different implementations and adaptations of the SGM algorithm in combination with a fronto-parallel sweeping direction in order to also account for non-fronto-parallel surfaces. |  | SGMΠ | SGM${}^{\Pi\text{-}\mathrm{sn}}$ | SGM${}^{\Pi\text{-}\mathrm{pg}}$
---|---|---|---|---
DTU | L1-abs | 19.832 | 19.768 | 19.684
(in mm) |  $\pm$16.225 |  $\pm$16.192 |  $\pm$16.154
L1-rel | 0.027 | 0.027 | 0.027
(in mm) |  $\pm$0.021 |  $\pm$0.021 |  $\pm$0.021
3DOMCity | L1-abs | 14.615 | 14.673 | 15.074
(in mm) |  $\pm$6.254 |  $\pm$6.229 |  $\pm$6.133
L1-rel | 0.012 | 0.012 | 0.012
(in mm) |  $\pm$0.007 |  $\pm$0.007 |  $\pm$0.006

As described in Section 2.3.2, apart from the straight-forward combination of
SGM with the plane-sweep sampling (SGMΠ), this work comprises two further
extensions of the SGM algorithm that allow to account for non-fronto-parallel
surface structures: namely the incorporation of surface normals to adjust the
zero-cost transition in the SGM path aggregation
(SGM${}^{\Pi\text{-}\mathrm{sn}}$) and the penalization of deviations from the
gradient of the minimum cost path (SGM${}^{\Pi\text{-}\mathrm{pg}}$). In
addition, by selecting an appropriate normal vector, and with it a
corresponding sweeping direction, it is also possible to adjust the sampling
plane orientations of the plane-sweep sampling to the scene structures. In the
following, the results achieved by SGM${}^{\Pi\text{-}\mathrm{sn}}$ and
SGM${}^{\Pi\text{-}\mathrm{pg}}$, in combination with a fronto-parallel
sampling plane orientation, are first evaluated and compared to those achieved
by SGMΠ, before the effects of different non-fronto-parallel plane
orientations are examined. The configuration of the other hyper-parameters is
as described and evaluated above: the size of the input bundle is set to
$|\Omega|=5$, the NCC with a support region of $5\times 5$ pixels is used. The
pyramid height is set to $n=3$ and $n=2$ for the DTU and 3DOMcity dataset,
respectively.

Figure 6: Qualitative comparison of the results achieved by the three
different SGM implementations on the DTU dataset. Row 1: Reference data from
the dataset, i.e. the ground truth depth and normal map, as well as the
reference image for which the data us computed. Rows 2 - 4: Data, i.e. depth,
normal and confidence maps, computed by SGMΠ, SGM${}^{\Pi\text{-}\mathrm{sn}}$
and SGM${}^{\Pi\text{-}\mathrm{pg}}$, respectively. As well as difference maps
holding the pixel-wise absolute difference between the estimated depth map and
the ground truth. The color encoding reaches from dark blue (low error) via
green to yellow (high error). The estimated maps are masked according to the
ground truth.  Figure 7: Qualitative comparison of the results achieved by the
three different SGM implementations on the 3DOMcity dataset. Row 1: Reference
data from the dataset, i.e. the ground truth depth and normal map, as well as
the reference image for which the data us computed. Rows 2 - 4: Data, i.e.
depth, normal and confidence maps, computed by SGMΠ,
SGM${}^{\Pi\text{-}\mathrm{sn}}$ and SGM${}^{\Pi\text{-}\mathrm{pg}}$,
respectively. As well as difference maps holding the pixel-wise absolute
difference between the estimated depth map and the ground truth. The estimated
maps are masked according to the ground truth. For visualization in this
figure, the resulting images have been rotated counterclockwise by
$90^{\circ}$. Thus, the color encoding of the normal maps differs from that
used in the other figures. Here, red represents an upwards orientation, while
green represents an orientation to the left.

The quantitative results displayed in Table 5 only reveal minor differences in
the L1 error between the different implementations of the SGM optimization.
While, in case of the DTU dataset, the best results are achieved by the
SGM${}^{\Pi\text{-}\mathrm{pg}}$ implementation, on the 3DOMcity dataset, the
standard adaptation of the SGM optimization to the plane-sweep sampling, i.e.
SGMΠ, reaches the lowest error. The relative L1 error does not reveal any
difference. This is due to the fact that the individual L1-rel scores only
start to differ from the forth decimal place on-wards. Nonetheless, the
ranking with respect to L1-rel score corresponds to that of the L1-abs score.
In a qualitative comparison, Figure 6 reveals that
SGM${}^{\Pi\text{-}\mathrm{sn}}$ leads to a seemingly smoother depth and
normal map (e.g. on the ground plane), while at the same time, however,
loosing small details and amplifying unwanted depth discontinuities in some
areas, such as the building facade. When closely comparing the normal maps
between SGMΠ and SGM${}^{\Pi\text{-}\mathrm{pg}}$, slightly smaller stair-
casing artifacts can be noticed in case of SGM${}^{\Pi\text{-}\mathrm{pg}}$,
which also supports the slightly lower error in Table 5. A qualitative
comparison between the results on the 3DOMcity dataset in Figure 7, however,
does not reveal any noticeable differences between the different
implementations. The reason for the small L1-abs error achieved by SGMΠ on the
3DOMcity dataset is assumed to be due to the fact, that the 3DOMcity dataset
also contains a subset of nadir images, in which only a small number of
slanted surfaces are existing and the fronto-parallel orientation of the
sampling planes coincides with most of the scene structure.

Figure 8:  ROC curves illustrating the error rate achieved by the three
different SGM implementations as a function of increasing density of the
estimated depth map.

To further quantify the strength and weaknesses of the three different SGM
aggregation strategies, three ROC curves, one for each extension, are plotted
for each dataset in Figure 8. In this, the curves illustrate the error rate
achieved by the corresponding SGM extension as a function of increasing
density of the estimated depth map. The density of the depth map is varied by
sampling the number of pixels in steps of $5$ % based on their ordered
confidence stored in $\mathcal{C}$, going from a high to a low confident
estimate. The mean error rate is quantified by $1\text{-}\text{Acc}_{1.05}$
(cf. Equation 20) and states the number of sampled pixels in $\mathcal{D}$
whose absolute difference to the ground truth exceeds $5$ % of the ground
truth value. Thus, at a low density of $\mathcal{D}$, i.e. a high confidence
threshold, the error rate should ideally be at its minimum and then increase
with increasing density, reaching the overall error of $\mathcal{D}$ at a
density of $100$ %. The plots start off at a density of $5$ %, since the error
rate at a density of $0$ % is undefined. The curves in Figure 8 further
support the superiority of the surface-aware SGM extensions over the standard
SGM adaptation to plane-sweep sampling, as for both datasets the curves of
SGM${}^{\Pi\text{-}\mathrm{sn}}$ and SGM${}^{\Pi\text{-}\mathrm{pg}}$ are
below that of SGMΠ, illustrating smaller error rates. The fact that most of
the ROC curves start of on a high error rate at a density of $5$ % and then
drop down before increasing again, suggests that the estimated confidence
values do not represent the certainty of the depth estimates appropriately.
The reasons for this are manifold and this is further discussed in Section
4.4.

Table 6: Run-time and error measurements of the three SGM extensions with $8$ and $4$ aggregation paths each, conducted on the DTU dataset with an image size of $1600\times 1200$ pixels. |  | SGMΠ | SGM${}^{\Pi\text{-}\mathrm{sn}}$ | SGM${}^{\Pi\text{-}\mathrm{pg}}$
---|---|---|---|---
$8$-Path SGM | L1-abs | 19.832 | 19.768 | 19.684
(in mm) |  $\pm$16.225 |  $\pm$16.192 |  $\pm$16.154
L1-rel | 0.027 | 0.027 | 0.027
(in mm) |  $\pm$0.021 |  $\pm$0.021 |  $\pm$0.021
Run-time | 640 | 895 | 2079
(in ms) |  $\pm$3 |  $\pm$2 |  $\pm$3
$4$-Path SGM | L1-abs | 21.072 | 21.091 | 20.908
(in mm) |  $\pm$16.634 |  $\pm$16.632 |  $\pm$16.599
L1-rel | 0.029 | 0.029 | 0.029
(in mm) |  $\pm$0.022 |  $\pm$0.022 |  $\pm$0.022
Run-time | 413 | 546 | 1132
(in ms) |  $\pm$2 |  $\pm$2 |  $\pm$3

Since the presented approach aims for incremental and online processing,
meaning that the computation should ideally keep up with the input stream, the
run-time of each SGM extension is also of interest. In Table 6, the run-time
and error of the complete approach, with above mentioned parameterization and
with respect to each of the three SGM implementations is listed. The
measurements were conducted on the dataset of the DTU benchmark. In addition
to the standard use of $8$ aggregation paths, which achieves the lowest error,
the run-time and error when using only $4$ aggregation paths is also listed.
In the latter case, the diagonal aggregation paths are omitted within the SGM
aggregation. This is motivated by the fact that a number of studies (Banz et
al., 2010; Hernandez-Juarez et al., 2016; Ruf et al., 2021) show that a
reduction in the number of aggregation paths from $8$ to $4$ greatly decreases
the computation time of the SGM aggregation, while only marginally increasing
its error, which is also supported by the numbers in Table 6. Furthermore, the
measurements reveal that especially the SGM${}^{\Pi\text{-}\mathrm{pg}}$
extension introduced a great computational complexity compared to SGMΠ and
SGM${}^{\Pi\text{-}\mathrm{sn}}$. However, the reduction of the aggregation
paths has a great impact on the run-time, reducing it by up to $45$ %, but
only marginally affecting the error. Whether the listed run-time is sufficient
for online processing is further discussed in Section 4.3.

In the second part of this experiment, the use of non-fronto-parallel plane
orientations within the plane-sweep sampling is investigated. In this, an
additional horizontal and vertical orientation, both with respect to the
reference coordinate system of the scene, were selected and compared to the
fronto-parallel sampling direction. For this, the DTU dataset is split into
two subsets, one for the horizontal sampling, in which the camera is looking
in a more downwards direction, as well as one for the vertical sampling, where
the camera pitch is smaller. As reference, for both subsets, results with a
fronto-parallel sampling were computed separately. The quantitative results
reveal a major increase in error when non-fronto-parallel plane orientations
are used for sampling, as also illustrated by the excerpt in Figure 9.
However, the Figure 9 also reveals that in areas where the surface structure
coincides with the sampling direction, e.g. the ground plane, the depth map is
very smooth and consistent.

Figure 9: Qualitative comparison between the use of a fronto-parallel and non-
fronto-parallel sampling direction in combination with SGMΠ. Columns 2 & 4:
Corresponding estimated depth map. Columns 3 & 5: Difference map holding the
pixel-wise absolute difference between the estimated depth map and the ground
truth. The color encoding reaches from dark blue (low error) via green to
yellow (high error). The estimated depth maps and the difference maps are
masked according to the ground truth.

###  3.6 Improvements Gained by Post-Filtering and Geometric Consistency

In a final ablative experiment, the effects of the implemented post-filtering
methods, to remove remaining outliers and supposedly wrong estimates by means
of Difference-of-Gaussian (DoG) filtering (cf. Section 2.6.1) and enforcing
geometric consistency (cf. Section 2.6.2) are studied. While the latter one
relies on the actual estimates, the DoG filter is based on the assumption that
image regions with low texture might lead to ambiguities in the image matching
and, in turn, wrong estimates, masking out corresponding regions. This,
however, might result in the false removal of good or even correct estimates.

Instead of using the absolute and relative L1 metric to quantitatively assess
the results achieved when employing post-filtering, the effects are evaluated
using the accuracy $\text{Acc}_{\theta}$ (cf. Equation 20) and completeness
$\text{Cpl}_{\theta}$ measure (cf. Equation 21). This is because they
indirectly include information on the density of the resulting depth maps,
which should ideally be as high as possible. Since the individual sequences of
the 3DOMcity dataset are made up of too little images in order to perform a
geometric consistency check with the parameterization mentioned in Section
2.6.2, this experiment is only conducted on the dataset of the DTU benchmark.
Figure 10 depicts the results of different post-filtering strategies, i.e. DoG
filtering, geometric consistency filtering, as well as a combination of both,
executed in combination with the three different SGM extensions and a fronto-
parallel sampling. In addition, the accuracy-completeness curve resulting from
the corresponding configurations without post-filtering is also displayed as
reference. In the construction of the curves, the threshold $\theta$ (i.e. the
threshold used to calculate the accuracy and completeness rates) is varied
within the list of $1.25,1.20,1.15,1.10,1.05\ \text{and}\ 1.01$. Note that
with a decreasing threshold, the accuracy and completeness rates drop. Thus,
the highest values are achieved with $\theta=1.25$.

Figure 10:  Accuracy-completeness curves of different post-filtering
strategies, i.e. DoG filtering, geometric consistency filtering, as well as a
combination of both, executed in combination with the three different SGM
extensions and a fronto-parallel sampling. In the construction, the threshold
$\theta$ is varied within the list of $1.25,1.20,1.15,1.10,1.05\ \text{and}\
1.01$. With a decreasing threshold, the accuracy and completeness rates drop.

Most evidently, Figure 10 again reveals that there is not much difference in
the overall error between the three SGM implementations. However, the
accuracy-completeness curves clearly show the differences that the individual
post-filtering strategies make. Unsurprisingly, the reference configuration
where no filtering is employed reaches the highest completeness, as no
estimates are removed from the predicted depth map, which, in turn, also leads
to the lowest accuracy. The use of a DoG filtering clearly improves this, as
it is likely to remove quite a number of wrong estimates originating from low-
textured areas. However, as expected, the DoG filter presumably also removes a
number of correct estimates, as the use of the filtering based on geometric
consistency achieves similar completeness, while at the same time reaching a
higher accuracy. Especially when considering the scores for $\theta=1.01$,
i.e. the lower left end of each curve, the use of a geometric consistency
filter achieves an increase in completeness of approximately $5$ %, while
exceeding the accuracy of the use of a DoG filter by more than $10$ %. A
distinct recommendation on which filter to use, however, cannot be concluded,
since both filtering strategies have their strengths and weaknesses,
especially with respect to online processing as discussed in Section 4.3. A
combination of both filters is not motivated. Even though the accuracy
slightly increases, the completeness drops partly by over $20$ %. Moreover,
this effect can also be achieved by decreasing the threshold of the
reprojection error $\eta_{\mathrm{r}}$ in the geometric consistency check,
which probably will increase the accuracy even more.

Lastly, to directly compare the different SGM extensions in combination with
the geometry-based filtering, which achieves the best results, the
corresponding F-scores for each evaluated $\theta$ (cf. Equation 22) are
listed in Table 7. Just as the results displayed in Table 5, the F-scores
reveal the superiority of SGM${}^{\Pi\text{-}\mathrm{pg}}$ over the other two
implementations, since for all $\theta$, except one,
SGM${}^{\Pi\text{-}\mathrm{pg}}$ reaches the highest F-score.

Table 7: F-scores achieved by the three SGM extensions in combination with the post-filtering based on geometric consistency. | $\boldsymbol{\text{{F}}_{1.25}}$ | $\boldsymbol{\text{{F}}_{1.20}}$ | $\boldsymbol{\text{{F}}_{1.15}}$ | $\boldsymbol{\text{{F}}_{1.10}}$ | $\boldsymbol{\text{{F}}_{1.05}}$ | $\boldsymbol{\text{{F}}_{1.01}}$
---|---|---|---|---|---|---
| (in %) | (in %) | (in %) | (in %) | (in %) | (in %)
SGMΠ | $74.2$ | $74.1$ | $74.0$ | $73.7$ | $71.5$ | $51.9$
SGM${}^{\Pi\text{-}\mathrm{sn}}$ | $74.1$ | $74.1$ | $74.0$ | $73.6$ | $71.5$ | $52.3$
SGM${}^{\Pi\text{-}\mathrm{pg}}$ | $75.6$ | $75.5$ | $75.4$ | $75.1$ | $72.9$ | $51.4$