# Detecting and approximating decision boundaries in low dimensional spaces

M. Grajewski corresponding author. Email<EMAIL_ADDRESS>FH Aachen
University of Applied Sciences, Faculty of Medical Engineering and
Technomathematics, Heinrich-Mußmann-Str. 1, 52428 Jülich, Germany  Institute
for Data-Driven Technologies, FH Aachen University of Applied Sciences,
Heinrich-Mußmann-Str. 1, 52428 Jülich, Germany  A. Kleefeld  Forschungszentrum
Jülich GmbH, Jülich Supercomputing Centre, Wilhelm-Johnen-Str., 52425 Jülich,
Germany  FH Aachen University of Applied Sciences, Faculty of Medical
Engineering and Technomathematics, Heinrich-Mußmann-Str. 1, 52428 Jülich,
Germany

###### Abstract

A method for detecting and approximating fault lines or surfaces,
respectively, or decision curves in two and three dimensions with guaranteed
accuracy is presented. Reformulated as a classification problem, our method
starts from a set of scattered points along with the corresponding
classification algorithm to construct a representation of a decision curve by
points with prescribed maximal distance to the true decision curve. Hereby,
our algorithm ensures that the representing point set covers the decision
curve in its entire extent and features local refinement based on the
geometric properties of the decision curve. We demonstrate applications of our
method to problems related to the detection of faults, to Multi-Criteria
Decision Aid and, in combination with Kirsch’s factorization method, to
solving an inverse acoustic scattering problem. In all applications we
considered in this work, our method requires significantly less pointwise
classifications than previously employed algorithms.

keywords: fault detection, fault approximation, inverse scattering problem,
MCDA

##  1 Introduction

Let us consider a piecewise constant function
$f:\mathbb{R}^{m}\supset\Omega\to\\{1,2,\ldots,n\\}$ with $\Omega$ being
compact, simply connected and equipped with a piecewise smooth boundary,
$m\in\\{2,3\\}$. Such $f$ subdivides $\Omega$ into mutually disjoint subsets
$\Omega_{i}:=f^{-1}(i)$ with $\Omega=\cup_{i=1}^{n}\overline{\Omega_{i}}$ and
$\Omega_{i}\cap\Omega_{j}=\emptyset$, if $i\neq j$. We assume that each
$\Omega_{i}$ features a piecewise smooth boundary. We are interested in
approximating $\Gamma_{i,j}:=\overline{\Omega_{i}}\cap\overline{\Omega_{j}}$
relying on as few evaluations of $f$ as possible and present in this work an
algorithm for this task. We choose this quantity as a measure of efficiency
because evaluating $f$ can be arbitrarily costly in applications and dominates
the runtime in such a case. Our problem can be immediately understood as a
classification problem, such that we identify $\Omega_{i}$ with a class $i$
and interpret the curves or surfaces of discontinuity $\Gamma_{i,j}$ as
decision curves or surfaces.

One field of application is economics and operations research. Multicriteria
Decision Aid (MCDA) methods can help a decision maker to choose the best one
from a finite number of alternatives based on different, even conflicting
criteria. MCDA methods assume that a decision depends on quantifiable
parameters $(x_{1},\ldots,x_{m})$ (“input factors”) and is drawn
deterministically. For an overview over various MCDA approaches, applications
and case studies, we refer to [14, 25] among many others and the references
cited therein. In this context, $\Omega_{i}$ is the set of all input factors
that lead to the decision for the $i$-th alternative in the MCDA method.
Analysing the decision process with respect to the input factors means
consistently describing all $\Omega_{i}$ based upon evaluating the MCDA method
for arbitrary combinations of input factors. This can be achieved by
approximating all $\Gamma_{i,j}$.
Computing a reconstruction of an obstacle in three dimensions is a field of
application in acoustic scattering theory. More precisely, one wants to
determine the support of an inhomogeneous object (its boundary to be exact)
from measured far-field data which typically is a desired task in non-
destructive testing. The far-field data are obtained for different incident
waves and measured points on the unit sphere. Several reconstruction
algorithms to find the boundary of the unknown inhomogeneity are available
such as iterative methods [22], decomposition methods [31] and sampling/probe
methods (see [26] for a detailed overview). The latter ones can be further
categorised into the linear sampling method, the generalized linear sampling
method, the factorization method, the probe method and variants of it (refer
to [12, 6, 20, 18, 26], respectively). However, here we will focus on the
classical factorization method for the acoustic transmission problem, refer
also to [4], with which one can decide if a given point is located inside or
outside the obstacle. Therefore, the factorization method transfers the
reconstruction of an obstacle to a classification problem and thus into a
field of application of our method. Note that the far-field data within [4]
has also been used in [8] and [17].

Finding and approximating the sets $\Gamma_{i,j}$, sometimes called fault
lines, is important in the exploration of natural resources. The presence and
the location of $\Gamma_{i,j}$ can provide useful insights for exploring and
later on exploiting of oil reservoirs [16] and play a significant role in some
geophysical applications [15]. The underlying mathematical problem is closely
related to ours, albeit not the same, as usually, the function $f$ considered
does not provide integer values as in our case. Therefore, an additional
algorithm for detecting fault lines is required then, and the classification
of a single point may be not trivial anymore as it is in our case. Moreover,
algorithms for fault detection may need to deal with noisy data (e.g. [10]),
whereas we consider certain data only.

Many algorithms for detecting and approximating the sets $\Gamma_{i,j}$ have
been proposed, like [5, 16, 15, 1, 10] among many others, which all feature
strengths and weaknesses. However, the vast majority of these approaches
restrict to the 2D case, whereas we present a method for 2D and 3D. The
algorithm proposed in [1] and the work cited therein was the starting point
for our research. Classification is one of the standard problems in Machine
Learning. There are a lot of powerful and versatile algorithms available which
could readily applied to our problem; we refer to [9] for an overview. These
methods are however designed for uncertain data in high dimensional spaces,
whereas we consider secure data in low dimensional spaces, a completely
different use case.

Our method approximately describes the $\Gamma_{i,j}$ by providing a set of
points with a guaranteed maximal normal distance to $\Gamma_{i,j}$. These
points are intended for constructing a polygon (2D) or a surface triangulation
in 3D. While there are more sophisticated and elegant ways of describing these
sets, it allows us to (approximately) replace an actual classification by a
simple and fast point-in-polygon test.

This article is organised as follows: We describe our algorithm for 2D and 3D
in Section 2 and elaborate on the direct and inverse acoustic scattering
problem, one of our applications, in Section 3. We present results and
applications of our algorithm to the detection of faults, to decision modeling
and to inverse scattering in Section 4 and finally conclude in Section 5.

##  2 The algorithm

###  2.1 Detection and Approximation of $\Gamma_{i,j}$ in 2D

In this section, we present our algorithm for approximating
$\Gamma_{i,j}\neq\emptyset$ for $m=2$ by sufficiently many well distributed
and ordered points sufficiently close $\Gamma_{i,j}$. This implicitly
constitutes a polygonal description of $\Gamma_{i,j}$. For convenience, we
provide flow charts (Figs. 1 and 2) and describe selected building blocks in
detail. In the flow charts, actual building blocks are typed in monospace
lettering.

Figure 1: General flow chart.

Figure 2: Flow charts of the Algorithms initialise, fill, expand and adapt.
Note that Building block barymeans inside initialise is employed twice, albeit
to different sets.

######  Algorithm 2.1 (initialise).

This algorithm aims at providing initial approximations to any $\Gamma_{i,j}$;
we refer to Fig. 1 (right) for an overview. From now on, we assume
$\Gamma_{i,j}\neq\emptyset$. In Building block initialset, we sample $f$ on
$\Omega$ rather coarsely and obtain an initial point set $X$ along with the
corresponding classification information. Building block barymeans creates
additional sampling points in the vicinity of any $\Gamma_{i,j}$. Following
Allasia et al. [1], we employ a $k_{\mathrm{near}}$-nearest neighbour
approach: For any $x\in X$, let $N(x)$ be the set of the $k_{\mathrm{near}}$\-
nearest neighbours of $x$ in $X$. If $N(x)\cap\Omega_{i}\neq\emptyset$ and
$N(x)\cap\Omega_{j}\neq\emptyset$, we consider $x$ close to $\Gamma_{i,j}$.
Let $N_{\ell}(x)=N(x)\cap\Omega_{\ell}$ for some $\ell$. Hence,
$N(x)=\bigcup_{i=1}^{r}N_{c_{i}}(x)$ for certain indices $c_{1},\ldots,c_{r}$,
$N_{c_{i}}(x)\neq\emptyset$. We compute the barycentres $b_{c_{i}}$ of all
$N_{c_{i}}(x)$, $1\leq i\leq r$, and then their arithmetic means
$y_{c_{i},c_{j}}=0.5(b_{c_{i}}+b_{c_{j}}),1\leq i<j\leq r$. Let $M(x)$ be the
set of all $y_{c_{i},c_{j}}$ generated from $N(x)$. If
$N(x)\subset\Omega_{\ell}$ for some $\ell$, we set $M(x)=\emptyset$. We end up
with $\mathcal{M}=\bigcup_{x\in X}M(x)$. By definition of $\mathcal{M}$, there
are no duplicate points; however, a practical implementation requires removing
duplicates. After classifying the points in $\mathcal{M}$, we repeat barymeans
on $\mathcal{M}$ obtaining sets $M^{2}(x)$ for any $x\in\mathcal{M}$, then
$\mathcal{M}^{2}=\bigcup_{x\in\mathcal{M}}M^{2}(x)$ and ultimately a further
enriched set of sampling points
$\overline{X}=X\cup\mathcal{M}\cup\mathcal{M}^{2}$.
Building block iniapprox computes initial approximations for all
$\Gamma_{i,j}$. For any $x\in(\mathcal{M}\cup\mathcal{M}^{2})\cap\Omega_{i}$,
we search the nearest point $x^{\prime}\in\overline{X}\cap\Omega_{j},j>i$. We
use $x$ and $x^{\prime}$ as starting points for a bisection algorithm on the
line $\overline{x^{\prime}x}$. If the bisection algorithm is successful, we
end up with a point pair $x^{(i)}\in\Omega_{i}$ and $x^{(j)}\in\Omega_{j}$
with $\|x^{(i)}-x^{(j)}\|\leq 2\varepsilon_{b}$, where $\varepsilon_{b}$ is a
user-prescribed threshold. Then, the distance of
$x^{(i,j)}=0.5(x^{(i)}+x^{(j)})$ to $\Gamma_{i,j}$ is at most
$\varepsilon_{b}$. We subsume the bisection process up to $\varepsilon_{b}$
and computing $x^{(i,j)}$ from its results in Building block bisection. From
now on, we consider point triplets for approximating $\Gamma_{i,j}$ only; for
any such triplet $x$, the superscript $(i)$ denotes the point in $\Omega_{i}$,
the superscript $(j)$ its counterpart in $\Omega_{j}$ and the superscript
$(i,j)$ the arithmetic mean of the two points. We end up with a set of
triplets $\widetilde{S}_{i,j}$. We moreover set
$\widetilde{S}_{i,j}^{(i)}=\\{x^{(i)}\>|\>x\in\widetilde{S}_{i,j}\\}$ and
$\widetilde{S}_{i,j}^{(j)}=\\{x^{(j)}\>|\>x\in\widetilde{S}_{i,j}\\}$ It may
occur that some triplets in $\widetilde{S}_{i,j}$ are tightly clustered. We
thin such clusters as they add to complexity but not to accuracy by removing
appropriate triplets (Fig. 3). After cluster removal, initialise provides sets
of triplets $S_{i,j}$.

######  Remark 2.2.

Building block initialise detects $\Gamma_{i,j}\neq\emptyset$ by
$\widetilde{S}_{i,j}\neq\emptyset$. However, depending on $X$, it may happen
that $\widetilde{S}_{i,j}=\emptyset$ even if $\Gamma_{i,j}\neq\emptyset$.
Reliably detecting all non-empty $\Gamma_{i,j}$ depends on a sufficiently
large $X$ and thus ultimately on the user.

######  Test problem 2.3.

For $\Omega=[0,1]^{2}$, we consider the following partition: Let
$\Omega_{3}:=\\{(x-1)^{6}+(y-0.5)^{6}<0.005\\}\cap\Omega$,
$\Omega_{1}^{\prime}:=\\{y\leq 0.7+0.1\sin(10\pi x^{1.5})\\}\cap\Omega$ and
$\Omega_{2}^{\prime}:=\Omega\setminus\Omega_{1}^{\prime}$. Then, we set
$\Omega_{1}:=\Omega_{1}^{\prime}\setminus\Omega_{3}$ and
$\Omega_{2}:=\Omega_{2}^{\prime}\setminus\Omega_{3}$ and study the partition
$\Omega=\Omega_{1}\cup\Omega_{2}\cup\Omega_{3}$ (Fig. 4, left).

######  Example 2.4.

For Test problem 2.3, we choose an initial sampling set $X$ consisting of 50
Halton-distributed points (Fig. 4, left) and set $k_{\mathrm{near}}=10$ and
$\varepsilon_{b}=0.001$. The set $\mathcal{M}$ consists of 47 additional
points represented as black dots; $\mathcal{M}^{2}$ contains 49 points
displayed in blue. $\widetilde{S}_{1,2}$ contains 43 triplets,
$\widetilde{S}_{1,3}$ contains 26, and $\widetilde{S}_{2,3}$ contains 7 (Fig.
4, right). After cluster removal, we obtain $|S_{1,2}|=23$, $|S_{1,3}|=13$ and
$|S_{2,3}|=4$ (Fig. 3 and Fig. 4, right).

Figure 3: Cluster removal by removeclusters in Example 2.4 for
$\widetilde{S}_{1,2}$ (greyed out); points in $S_{1,2}^{(1)}$ are displayed in
red. This figure is an excerpt of Fig. 4, right.

Figure 4: Left: Partition of $\Omega=[0,1]^{2}$ according to Example 2.4
indicated by grey solid lines with the initial point set $X$ displayed as grey
crosses. We moreover show $\mathcal{M}$ (black points) and $\mathcal{M}^{2}$
(blue points); Right: $S_{1,2}^{(1)}$ (red points), $S_{1,3}^{(1)}$ (blue
points), and $S_{2,3}^{(2)}$ (cyan-coloured points).

The triplets in $S_{i,j}$ usually provide an incomplete approximation of
$\Gamma_{i,j}$ only, feature gaps and lack ordering (e.g. Fig. 4, right).

######  Remark 2.5.

The bisection-based approach in Building block iniapprox implicitly assumes
that $\overline{x^{\prime}x}$ intersects $\Gamma_{i,j}$ and therefore may fail
if this does not hold (Fig. 5). However, such failure modes occur only rarely
in practical computations, and we implemented fallbacks in that case.

Figure 5: Building block bisection can fail: Finding points near
$\Gamma_{1,3}$ fails for the starting points $x\in\Omega_{3}$ and
$x^{\prime}\in\Omega_{1}$ shown, as $\overline{x^{\prime}x}$ does not
intersect $\Gamma_{1,3}$.

Algorithm fill (for an overview, we refer to the flow chart in Fig. 2)
provides triplets with a maximal user-prescribed distance $\varepsilon_{b}$ to
$\Gamma_{i,j}$ and maximal mutual distance $\varepsilon_{\mathrm{gap}}$ based
upon $S_{i,j}$, i.e. the result of Algorithm 2.1 (initialise). As we assume
that $\varepsilon_{b}\ll\varepsilon_{\mathrm{gap}}$, we define the distance of
two triplets $x,y\in S_{i,j}$ as $\|x^{(i)}-y^{(i)}\|$ and the distance of $x$
to $\Gamma_{i,j}$ as the distance of $x^{(i,j)}$ to $\Gamma_{i,j}$. Following
a bottom-up approach, we start with discussing selected Building blocks
employed in fill.

######  Building block 2.6 (sort).

This Building block is to sort a set of triplets $S$ according to their
position along $\Gamma$. We omit the indices $i$ and $j$ for clarity. Let us
assume that $\Gamma$ is piecewise smooth and fulfills an inner cone condition
with angle $\beta_{\mathrm{angle}}$. We first search a triplet
$x_{\mathrm{start}}\in S$ closest to the boundary of $\Omega$ and assume
$x_{\mathrm{start}}$ to be the first triplet in the sorted set. We initially
set $\check{S}=\\{x_{\mathrm{start}}\\}$. Let now the triplets in
$\check{S}\subset S$ be already sorted with $x_{r}$ being the last of those,
$r>1$. We consider the $k_{\mathrm{sort}}$ nearest triplets
$y_{1},\ldots,y_{k_{\mathrm{sort}}}$ to $x_{r}$ in $S\setminus\check{S}$,
sorted by increasing distance to $x_{r}$. If for $s=1$,
$\angle(x_{r}-x_{r-1},x_{i}-y_{s})<\beta_{\mathrm{angle}}$, we set
$x_{r+1}=y_{s}$ and add $x_{r+1}$ to $\check{S}$; otherwise, we repeat with
$s+1$. If $s>k_{\mathrm{sort}}$, we store $\check{S}$, set
$S=S\setminus\check{S}$, and repeat the sorting procedure until $S=\emptyset$
ending up with a finite number of disjoint ordered subsets. At the end, we
combine all sorted subsets based upon the Euclidean distance between first and
last points of the subsets reversing the order of a subset if necessary.

Figure 6: Sorting due to Allasia may fail (left; sorting indicated in blue),
whereas our method succeeds in the present situation (right). We obtain two
ordered subsets (direction of sorting is indicated by arrows which after
combination yields the correct ordering. The dashed line represents a
connection to the nearest neighbour rejected due to angle in our approach.

######  Remark 2.7.

Allasia et al. [1] present a simpler sorting method than ours as they do not
enforce the condition
$\angle(x_{r}-x_{r-1},x_{i}-y_{s})<\beta_{\mathrm{angle}}$. However, it may
fail if $x_{\mathrm{start}}$ is not the true starting point and if
additionally the points are unevenly distributed along $\Gamma$ (Fig. 6) in
contrast to ours. Of course, our sorting method can fail as well, but due to
our experience, it is more reliable than Allasia’s method and works
sufficiently well. There are many more sophisticated sorting methods based
e.g. upon graph theory [2, 3, 13, 24], which are more reliable than our
approach, but more time-consuming and much harder to implement.

For describing Building block inipairs, which is part of fill, we introduce
another two Building blocks, which will be used in adapt as well.

######  Building block 2.8 (estcurv).

Let be $S_{\mathrm{loc}}=\\{x_{1},\ldots,x_{r}\\},r>2$, a set of ordered
points near $\Gamma$ up to $\varepsilon_{b}$. This Building block estimates
the curvature $c_{\ell}$ of $\Gamma$ in $x_{\ell},1<\ell<r$ by least-squares
fitting an approximation using Gaussian radial basis functions (RBFs) and then
$c_{\ell}$ by the curvature of that approximation in $x_{\ell}$. We hereby
assume that after shifting and suitable rotation, $\Gamma$ can be locally
represented as a graph of an unknown function. As the points in
$X_{\mathrm{loc}}$ are located on $\Gamma$ only up to $\varepsilon_{b}$, we
penalise the second derivative of the RBF approximation subject to a maximal
residual of $\varepsilon_{b}$. This coincides with the maximal deviation in
the value at an approximation point. We employ Tikhonov regularization with
parameter estimation using Morozov’s discrepancy principle.

If $\Gamma$ cannot be considered a graph of a function even after rotation, we
draw as a fallback a circle through the points with indices $\ell-1$, $\ell$
and $\ell+1$ and use the inverse of its radius for estimating $c_{\ell}$. We
estimate $c_{\ell}$ in the first or last point of $S_{\mathrm{loc}}$ by
drawing a circle through the first or last three points in $S_{\mathrm{loc}}$.

######  Building block 2.9 (esterror).

This Building block estimates the maximal deviation $\delta$ of a smooth curve
from a straight line between two points on the curve with distance $d$. A
straightforward calculation reveals that

$\delta=0.25cd^{2}+1/16c^{3}d^{4}+{\cal{O}}(d^{6}),$ (1)

where $c$ denotes the maximal curvature of the curve between the two points.
For $S_{\mathrm{loc}}$ as in Building block 2.8, we estimate the maximal
deviation $\delta$ of $\Gamma$ from the straight line between consecutive
points $x_{\ell}$ and $x_{\ell+1}$ by replacing $c$ with
$\max\\{c_{\ell},c_{\ell+1}\\}$ from Building block 2.8 (estcurv). Hereby, we
rely on (1) and neglect higher order terms.

Now we are prepared to discuss fill.

######  Algorithm 2.10 (fill).

Building block 2.6 (sort) sorts all triplets in $S_{i,j}$ according to their
position near $\Gamma_{i,j}$. We detect gaps in the representation of
$\Gamma_{i,j}$ by $S_{i,j}$ by considering subsequent triplets $x_{\ell}$ and
$x_{\ell+1}$. If the distance $d_{\ell}$ of $x_{\ell}$ to $x_{\ell+1}$ is
larger than a user-prescribed threshold $\varepsilon_{\mathrm{gap}}$, we
consider this a gap and aim to equidistantly add $R=\lceil
d_{\ell}/\gamma\rceil$ triplets near $\Gamma_{i,j}$ between $x_{\ell}$ and
$x_{\ell+1}$. To do so, Building block inipairs places new points
$z_{\ell,r},1\leq r\leq R$, equidistantly on
$\overline{x_{\ell}^{(i,j)}x_{\ell+1}^{(i,j)}}$ and computes from these
initial point pairs

$x^{+}_{\ell,r}=z_{\ell,r}+\alpha n,\quad x_{\ell,r}^{-}=z_{\ell,r}-\alpha
n,\quad 1\leq r\leq R.$

Here, $n$ denotes the (estimated) outer normal unit vector of $\Omega_{i}$
near $x_{\ell,r}$. Applying Building block 2.9 (esterror) to
$S_{\mathrm{loc}}=\\{x_{\ell-2}^{(i,j)},\ldots x_{\ell+2}^{(i,j)}\\}$ and some
safeguarding leads to

$\alpha=\min\\{\varepsilon_{\mathrm{safemax}}d_{\ell},\max\\{\delta,\varepsilon_{\mathrm{safemin}}\varepsilon_{b}\\}\\}$
(2)

with user-prescribed safety factors $\varepsilon_{\mathrm{safemax}}$ and
$\varepsilon_{\mathrm{safemin}}$.

######  Remark 2.11.

Having a local RBF approximation of $\Gamma$ at hand when computing $c$ in
estcurv, it seems to be straightforward for efficiency reasons to choose
points $z_{\ell,r}$ on that RBF approximation instead of just subdividing a
straight line. Numerical experiments did not show any significant advantage of
that approach compared to ours for fill. This could be related to the uneven
distribution of points on $\Gamma$ near gaps to fill, which may decrease the
quality of approximation. Therefore, we stick to the easier approach presented
here. However, estimating the curvature using Building block 2.8 (estcurv) is
sufficiently reliable for efficiently computing starting pairs.

However, the pairs of starting points (aka starting pairs) obtained from
inipairs are not necessarily valid. We call a starting pair for approximating
$\Gamma_{i,j}$ valid, if one of its points belongs to $\Omega_{i}$ and the
other one to $\Omega_{j}$. Therefore, we introduce Building block 2.12
(startpairs).

######  Building block 2.12 (startpairs).

This Building block obtains a valid starting pair from a pair of points
$(x_{\ell,r}^{+},x_{\ell,r}^{-})$. If the starting pair is already valid, we
return it as the result. If $x_{\ell,r}^{+}$ or $x_{\ell,r}^{-}$ belongs to a
third class, we stop without result. If both points belong to the same class,
we reflect $x_{\ell,r}^{+}$ on $z_{\ell,r}$ obtaining $x_{\ell,r}^{\prime}$
(Fig. 7). If both $(x_{\ell,r}^{+},x^{\prime}_{\ell,r})$ and
$(x_{\ell,r}^{-},x^{\prime}_{\ell,r})$ still belong to the same class, we
repeat this process with changing roles and escalating distances at most
$k_{\mathrm{rep}}$ times (typically $k_{\mathrm{rep}}=3$) and stop, if any of
the resulting point pairs is valid or one of the points belongs to a third
class.

######  Remark 2.13.

  1. 1.

We iterate the process of filling gaps in fill, as the arclength of $\Gamma$
between two subsequent triplets may considerably exceed the length of the
straight line between them, such that after a first pass, the mutual distance
of subsequent triplets may still exceed $\varepsilon_{\mathrm{gap}}$ in some
cases.

  2. 2.

If only two or even less triplets are known on $\Gamma_{i,j}$, estimating
curvature is impossible with estcurv, and computing $\alpha$ by (2) fails. For
this case, we implemented fallbacks.

Using bisection, fill creates new triplets, yielding new sets
$\overline{S}_{i,j}$.

Figure 7: Creating valid starting pairs in fill with startpairs: While
$(x_{\ell,1}^{+},x_{\ell,1}^{-})$ is valid, $(x_{\ell,2}^{+},x_{\ell,2}^{-})$
is not, as both points belong to the same class. However,
$(x_{\ell,2}^{-},x^{\prime}_{\ell,2})$ is valid. Figure 8: Algorithm 2.10
(fill) for Test problem 2.4. Previously existing points are greyed out in
$\overline{S}_{i,j}^{(i)}$. Otherwise, we stick to the colouring scheme of
Fig. 4.

######  Example 2.14.

For Test problem 2.3, we ordered the sets $S_{i,j}$ according to Building
block 2.6 with $\beta_{\mathrm{angle}}=\arccos(-0.9)\approx 154^{\circ}$ and
$k_{\mathrm{sort}}=5$. By filling gaps with $\varepsilon_{\mathrm{gap}}=0.05$,
$\varepsilon_{\mathrm{safemin}}=0.95$ and
$\varepsilon_{\mathrm{safemax}}=0.25$, we obtain $|\overline{S}_{1,2}|=33$,
$|\overline{S}_{1,3}|=18$, $|\overline{S}_{2,3}|=4$ (Fig. 8). We will stick to
the values of the parameters given here for all subsequent numerical examples.

######  Remark 2.15.

If $\Omega_{i}$ is not simply connected, some $\Gamma_{i,j}$ may consist of
several components. We detect this using fill. If there are some significant
gaps which can not be filled, it indicates the presence of several components.
We then subdivide $\overline{S}_{i,j}$ correspondingly and proceed with every
subset separately.

Fig. 8 indicates that even with filling gaps, $\overline{S}_{i,j}$ may not
appropriately represent $\Gamma_{i,j}$, as parts of $\Gamma_{i,j}$ before the
first known triplet $x_{1}$ and after the last known may be neglected.
Algorithm expand is used to expand $\overline{S}_{i,j}$ to a representation of
the complete curve $\Gamma_{i,j}$. It relies on several building blocks, which
we discuss first.

######  Building block 2.16 (extrapolate).

For a given ordered set $S$ with $n\geq 2$ triplets close up to
$\varepsilon_{b}$ to a curve $\Gamma$ and with average distance
$d_{\mathrm{avg}}$, we fit a polynomial with degree $n-1$ in local
coordinates. We compute these coordinates by least-squares-fitting a line to
$S$. As points in $S$ are located on $\Gamma$ up to $\varepsilon_{b}$ only, we
do not interpolate, but penalise the second derivative in a least-squares
approximation. Following Morozov’s discrepancy principle, we regularise such
that the maximal residual is approximately $\varepsilon_{b}$.

######  Building block 2.17 (stepsize).

Provided that $\Gamma$ extends sufficiently far before $x_{1}\in S$, it seems
to be straightforward to seek for a new triplet with distance
$\varepsilon_{\mathrm{gap}}$ to $x_{1}$. However, we limit the step size for
extrapolation based upon the curvature $c$ of $\Gamma$ in the vicinity of
$x_{1}$. Extrapolating far is unreliable in case of large curvature $c$, and
adapt will insert additional points afterwards in that region anyway for
accuracy reasons, such that it is much more efficient to adjust the step
length to the local properties of $\Gamma$ beforehand. Let us assume that a
polygonal final approximation of $\Gamma$ may deviate at most by
$\varepsilon_{\mathrm{err}}$ from $\Gamma$. We compute the step size which
would lead to a deviation of $\varepsilon_{\mathrm{err}}$ from a straight line
segment between $x_{1}$ and the new triplet yet to compute. This is a natural
upper bound for the step size $l_{\mathrm{extra}}$ in extrapolation.

Rearranging (1) and neglecting higher order terms leads to

$l_{\max}=\frac{2}{c^{2}}\left(\sqrt{1+4c\varepsilon_{\mathrm{err}}}-1\right)$
(3)

for the maximal admissible step length $l_{\max}$. However, evaluating (3) is
numerically unstable if $c\varepsilon_{\mathrm{err}}$ is small. We set
$(c\varepsilon_{\mathrm{err}})^{2}=v$ and search for the roots $v_{\min}$ and
$v_{\max}$ of $v^{2}+4v-16cd$. According to Vieta,
$v_{\min}=-\left(2+\sqrt{4+16cd}\right)$ and $v_{\max}=-16cd/v_{\min}$.
Resubstituting $v$ yields $l_{\max}=4\sqrt{d/(-cv_{\min})}$. Some safeguarding
of this result leads to a step length of

$l_{\mathrm{extra}}=\min\\{\varepsilon_{\mathrm{gap}},\beta_{\mathrm{growth}}d_{\mathrm{avg}},l_{\max}\\}\>,$

where we estimate $c$ using Building block 2.8 (estcurv) applied to
$S_{\mathrm{loc}}=\\{x_{1}^{i,j)},\ldots,x_{k_{\text{extra}}}^{(i,j)}\\}$ with
a user-defined parameter $k_{\text{extra}}$. The term
$\beta_{\mathrm{growth}}d_{\mathrm{avg}}$ increases robustness, as
extrapolation is reliable only near the points to extrapolate, and it may
happen that $d_{\mathrm{avg}}\ll\varepsilon_{\mathrm{gap}}$.

######  Algorithm 2.18 (expand).

This algorithm aims at finding triplets near $\Gamma_{i,j}$ beyond the first
or last known in $\overline{S}_{i,j}$ until the start or end of $\Gamma_{i,j}$
is reached or $\Gamma_{i,j}$ turns out to be a closed curve. For the sake of
simplicity, we refer in what follows to finding triplets before the first one
in $\overline{S}_{i,j}$. Finding triplets near $\Gamma_{i,j}$ beyond the last
one in $\overline{S}_{i,j}$ works analogously. We add a new triplet before
$x_{1}$ by extrapolating an approximation $\gamma$ of $\Gamma_{i,j}$ with
Building block 2.16 (extrapolate) setting
$S=\left\\{x_{1}^{(i,j)},\ldots,x_{k_{\mathrm{extra}}}^{(i,j)}\right\\}$ with
$x_{\ell}\in\overline{S}_{i,j}$ and choose the step size according to Building
block 2.17 (stepsize). This way, we obtain an extrapolating curve $\gamma$ and
some $s_{0}$ such that $\|x_{1}-\gamma(s_{0})\|\approx l_{\mathrm{extra}}$. We
then create a point pair $(x_{s_{0}}^{+},x_{s_{0}}^{-})$ based on
$\gamma(s_{0})$ in a similar way as in fill. Computing a valid starting pair
with Building block 2.12 (startpairs) and subsequent bisection yields a new
triplet in $\overline{S}_{i,j}$. We repeat this process until we reach the
true starting point of $\Gamma_{i,j}$. As heuristic criterion for
$\gamma(s_{0})$ exceeding this starting point, we consider

$x_{s_{0}}^{+}\not\in(\Omega_{i}\cup\Omega_{j})\vee
x_{s_{0}}^{-}\not\in(\Omega_{i}\cup\Omega_{j})$ (4)

(Fig. 9). In this case we employ Building block 2.19 (reducestepsize) for
obtaining a valid pair of points $(x_{\tilde{s}}^{+},x_{\tilde{s}}^{-})$,
which represents the starting point of $\Gamma_{i,j}$. After adding it to
$\overline{S}_{i,j}$, expand terminates. If $\Gamma_{ij}$ is closed, expanding
$\overline{S}_{i,j}$ as described above would lead to an endless loop.
Therefore, we start expanding $\overline{S}_{i,j}$, but detect after every
addition of a triplet, if $\Gamma_{i,j}$ is closed. If so, we stop expanding
and resort. We skip the details to keep the presentation uncluttered.
Algorithm expand yields approximating sets $\breve{S}_{i,j}$.

######  Building block 2.19 (reducestepsize).

Using the parameter value $s_{1}$ corresponding to $x_{1}^{(i,j)}$ as lower
bound and $s_{0}$ as upper bound, we obtain $\tilde{s}$, which leads to a
valid starting pair $(x_{\tilde{s}}^{+},x_{\tilde{s}}^{-})$ based upon
$\gamma(\tilde{s})$ and fulfils
$\|\gamma(\tilde{s})-\gamma(s^{\prime})\|<\varepsilon_{b}$ by bisection with
respect to condition (4). Here, $s^{\prime}$ denotes the second to last
parameter in the bisection process (compare Fig. 9).

Figure 9: Scheme of expanding $S_{1,2}$ until its end. Starting from
$x_{1}^{(1,2)},\ldots,x_{4}^{(1,2)}$ (displayed as blue dots), we construct
$\gamma$ by extrapolate. As $(x_{s_{0}}^{+},x_{s_{0}}^{-})$, displayed as
cyan-coloured dots, fulfills (4), we apply bisection with respect to the
parameter $s$ until a valid starting pair based upon $\gamma(\tilde{s})$ can
be constructed (displayed in magenta), from which we compute the final triplet
in $S_{1,2}$ by bisection.

Figure 10: Sets $\breve{S}_{i,j}^{(i)}$ (left) and $\hat{S}_{i,j}^{(i)}$
(right) for Test problem 2.4. We stick to the colouring scheme of Fig. 4.

######  Algorithm 2.20 (adapt).

Based on esterror, this algorithm inserts a triplet (approximately) halfway
between consecutive triplets $x_{\ell}$ and $x_{\ell+1}$, if esterror
indicates an error larger than $\varepsilon_{\mathrm{err}}$ and removes a
triplet, if the estimated error of both line segments the triplet belongs to
is smaller than $\varepsilon_{\mathrm{coarse}}$. In contrast to fill, we
employ the local RBF approximation from estcurv for computing an initial point
$x_{\mathrm{new}}$ between $x_{\ell}$ and $x_{\ell+1}$ when refining. With
$x_{\mathrm{new}}^{\pm}=x_{\mathrm{new}}\pm\alpha^{\prime}n$, we then proceed
as in fill. We compute $\alpha^{\prime}$ according to (2), but replace
$\delta$ by $\delta^{\prime}=1/16c^{3}d^{4}$, as the error due to (1) refers
to an approximation by line segments and is overly pessimistic for an RBF
approximation. For robustness, we never delete consecutive triplets in one
pass of the adaptive loop. After at most $k_{\mathrm{adap}}$ refinement and
coarsening sweeps, we end up with final sets $\hat{S}_{i,j}$.

######  Example 2.21.

For Test problem 2.3, we continue our calculations from Example 2.14. We set
$k_{\mathrm{extra}}=4,\varepsilon_{\mathrm{err}}=0.001$,
$\varepsilon_{\mathrm{coarse}}=0.0001$ and obtain $|\breve{S}_{1,2}|=43$,
$|\breve{S}_{1,3}|=28$, and $|\breve{S}_{2,3}|=18$ (Fig. 10, left). Algorithm
adapt with $k_{\text{adap}}=4$ yields $|\hat{S}_{1,2}|=77$,
$|\hat{S}_{1,3}|=26$, and $|\hat{S}_{2,3}|=22$ (Fig. 10, right). Computing
these sets requires 1088 classifications.

Figure 11: Reconstructed subdivision for Test problem 2.4.

###  2.2 Approximating $\Gamma_{i,j}$ in 3D

For approximating $\Gamma_{i,j}$ in three dimensions, we stick to the general
procedure for approximating these sets in two dimensions (Fig. 1). While
initialise remains unchanged, fill, expand and adapt differ from their 2D
counterparts, as these exploit ordering of the points on a decision curve or
fault line. There is however no straightforward ordering of points on a
surface. Because fill and adapt rely on esterror as in two dimensions, we
discuss error estimation first.

######  Building block 2.22 (esterror).

This Building block aims at estimating the maximal error $e$ of a linear
triangulation of $\Gamma$ based upon a finite set of points $S$ near $\Gamma$
up to $\varepsilon_{b}$. Let $\mathcal{T}$ be a Delaunay triangulation of
$S_{\text{loc}}\subset S$ and let us assume that $\Gamma$ can be represented
on the support of $\mathcal{T}$ by an unknown function $g$ after appropriate
change of coordinates. For some triangle $T\in\mathcal{T}$, it holds according
to [29, Theorem 4.1]

$\|g-I_{T}g\|_{T,\infty}\leq\frac{1}{2}\left(R^{2}-d^{2}\right)|g|_{2,\infty,T}.$

Here, $R$ describes the radius of the circumcircle of $T$, $d$ the distance
from its center to $T$ and $I_{T}g$ the Lagrange interpolant of $g$ on $T$;
with $|g|_{2,\infty,T}$, we denote the $L_{\infty}$-norm of the second
derivative of $g$ on $T$. It remains to estimate $|g|_{2,\infty,T}$. As $g$ is
unknown, we employ an RBF approximation $\varphi$ as in Building block 2.8
(estcurv) instead and approximate $|g|_{2,\infty,T}$ by evaluating its second
derivative in the vertices of $T$ and its center. The maximum of these four
values yields the desired approximation $\phi$ of $|g|_{2,\infty,T}$ and
therefore

$e\approx\frac{1}{2}\left(R^{2}-d^{2}\right)\phi$ (5)

######  Algorithm 2.23 (fill).

For any triplet $x$ in $S_{i,j}$, we search the $k_{\mathrm{near}}$ nearest
neighbours $x_{1},\ldots,x_{k_{\text{near}}}$ and switch to a local 2D
coordinate system by computing the optimal fitted plane in the sense that the
sum of the squared distances between the points and the plane is minimal (see
[27]). We then compute in local coordinates a Delaunay triangulation of the
patch $S_{\text{loc}}=\\{x,x_{1},\ldots,x_{k_{\text{near}}}\\}$. If the
maximal edge length of a triangle in this triangulation exceeds
$\varepsilon_{\text{gap}}$, we mark the centre of that triangle as starting
point for constructing a valid starting pair similar as in Building inipairs,
as long as the triangle is not too anisotropic, i.e. does not contain angles
approaching $0^{\circ}$ or $180^{\circ}$. However, large gaps in $S_{i,j}$ are
not covered by this local approach (Fig. 12). If there is a gap in $S_{i,j}$
and if the current point is at the boundary of this gap, it will be at the
boundary of the current patch. We detect this by exploiting the local
coordinate system. We assume that a patch in the vicinity of a large gap is at
least slightly elongated in tangential direction to the patch boundary.
Therefore, we consider the first local coordinate. If the first local
coordinate of the current point is almost the minimum or maximum of all
respective first coordinates of the patch, then it is at its boundary, and we
assume a large gap to fill. In this case, we construct $x_{\text{new}}$ on the
elongated line between the centre of gravity of the patch and $x$, enrich
$S_{\text{local}}$ by $x$ and continue as described above. Looping over all
triplets in $S_{i,j}$, removing duplicate starting points and starting points
which are extremely close enables us to compute triplets using Building block
bisection similar as in Algorithm 2.10 (fill). This yields an enriched
representation of $\Gamma_{i,j}$. Repeating this procedure with the enriched
set until no gaps are detected anymore leads to $\overline{S}_{i,j}$.

Figure 12: Left: Gradually filling large gaps in $S_{i,j}$ (displayed as grey
dots). The triplet $x$ is represented in golden colour; the triplets in
$S_{\text{local}}$, displayed in red, define a local coordinate system
sketched with black arrows. The centre of gravity of $S_{\text{local}}$ is
shown as a blue dot. Right: Expansion of $\overline{S}_{1,2}$ to the inner
boundary $\Gamma_{1,2}\cap\Gamma_{1,3}$. Triplets in
$\overline{S}_{1,2}\setminus\overline{S}_{1,2}^{I}$ are displayed as grey
dots, triplets in $\overline{S}_{1,2}^{I}$ as blue dots. 2D-Expanding
$\\{x,x^{\prime}\\}$ on $\mathcal{E}\cap\Gamma_{1,2}$ yields a new triplet
near $\Gamma_{1,2}\cap\Gamma_{1,3}$, displayed as red dot.

######  Test problem 2.24.

For $\Omega=[0,1]^{3}$, we consider the following partition: Let
$x_{M}=(1,0.5,0.5)^{\top}$ and
$\Omega_{3}:=\\{\|x-x_{M}\|_{6}<0.002\\}\cap\Omega$,
$\Omega_{1}^{\prime}:=\\{y+0.1z<0.7+0.1\sin(10x^{1.5})+0.05\sin(5z^{1.5})\\}\cap\Omega$
and $\Omega_{2}^{\prime}:=\Omega\setminus\Omega_{1}^{\prime}$. Then, we set
$\Omega_{1}:=\Omega_{1}^{\prime}\setminus\Omega_{3}$ and
$\Omega_{2}:=\Omega_{2}^{\prime}\setminus\Omega_{3}$ and study the partition
$\Omega=\Omega_{1}\cup\Omega_{2}\cup\Omega_{3}$ (Fig. 13).

Figure 13: Left: Partition of $[0,1]^{3}$ for Test problem 2.24 displayed by
$\Gamma_{1,2}$ (red), $\Gamma_{1,3}$ (blue) and $\Gamma_{2,3}$ (grey;
partially obscured by $\Gamma_{1,3}$); Right: approximating sets
$S_{i,j}^{(i)}$ after removing clusters.

######  Algorithm 2.25 (expand).

When expanding a representation $\overline{S}_{i,j}$ of $\Gamma_{i,j}$, we
distinguish between expansion towards inner boundaries
$\Gamma_{i,j}\cap\Gamma_{k,\ell}$ and expansion towards outer boundaries
$\Gamma_{i,j}\cap\partial\Omega$. We consider expanding to inner boundaries
first. Let $\overline{S}_{i,j}^{I}\subseteq\overline{S}_{i,j}$ be all triplets
closer than $0.75\varepsilon_{\text{gap}}$ to a triplet in another
$S_{k,\ell}$. For some $x\in\overline{S}_{i,j}^{I}$, let $x^{*}$ be the
triplet in $S_{k,\ell}$ closest to $x$. We estimate the normal vector
$\mathfrak{n}$ of $\Gamma_{i,j}$ in $x$ and project $\overline{xx^{*}}$ on the
corresponding tangential plane at $x$, yielding $\mathfrak{t}$. Let
$x^{\prime}\in\overline{S}_{i,j}$ be the nearest triplet to $x$ fulfilling
$\angle(-\mathfrak{t},\overline{xx^{\prime}})<\alpha_{\text{expand}}$ and
$\mathcal{E}$ the plane which contains $x$ and is spanned by
$\overline{xx^{\prime}}$ and $\mathfrak{n}$. We reduce the expansion to the
two-dimensional case by expanding the curve $\Gamma_{i,j}\cap\mathcal{E}$ with
Algorithm 2.18 (expand) from Section 2.1, using $x$ and $x^{\prime}$ as
initial set of triplets (compare Fig. 12, right).

For expanding to outer boundaries, we construct $\overline{S}_{i,j}^{B}$ as
follows: For each coordinate direction $i$, we select $n_{\text{expand},i}$
points in $S_{i,j}$ which are minimal or maximal with respect to this
coordinate and set
$\overline{S}_{i,j}^{B}=\overline{S}_{i,j}^{B}\setminus\overline{S}_{i,j}^{I}$
in order to avoid duplicate triplets. For determining $n_{\text{expand},i}$,
we rely on the axis-parallel bounding box of $\overline{S}_{i,j}$ with sizes
$b_{1},b_{2},b_{3}$. Then, $n_{\text{expand},i}=\lceil\max\\{b_{i+1\mod
3},b_{i\mod 3+1}\\}/\varepsilon_{\text{gap}}\rceil$. From these triplets, we
select the ones which are either closer than $\varepsilon_{\text{gap}}$ to one
of the boundary facets or fulfil

$\angle(\mathfrak{n},\mathfrak{n}_{\text{outer}})>\alpha_{\text{expbound}},$
(6)

where $\mathfrak{n}_{\text{outer}}$ stands for the outer normal vector of the
assigned boundary facet. Figure 15 illustrates the purpose of condition (6).
For any triplet in $\overline{S}_{i,j}^{B}$, we proceed as for expanding to
inner boundaries, but with $\mathfrak{t}=\mathfrak{n}_{\text{outer}}$. All new
triplets in the same facet $F$ of $\partial\Omega$ constitute the set
$\overline{S}_{i,j}^{F}$. As this set represents the curve $\Gamma_{i,j}\cap
F$ in the plane $F$, we are confronted with approximating a decision curve or
fault line in two dimensions. Therefore, we apply Algorithms 2.10 (fill) and
2.18 (expand) to $\overline{S}_{i,j}^{F}$. After adding the new triplets on
the boundary of $\Gamma$, we apply Algorithm 2.23 again and end up with an
enlarged set $\breve{S}_{i,j}$.

Figure 14: Left: Sets $\overline{S}_{1,2}^{(1)}$, $\overline{S}_{1,3}^{(1)}$
obtained applying Algorithm 2.23 (fill) for Test problem 2.24. Right:
$\breve{S}_{1,2}^{(1)}$, $\breve{S}_{1,3}^{(1)}$ after Algorithm 2.25 (expand)
for the same Test problem. Previously existing points are displayed in light
grey; we omit $\overline{S}_{2,3}^{(2)}$ and $\breve{S}_{2,3}^{(2)}$ in the
visualisation for the sake of clarity.

Figure 15: Suitable (left) and unsuitable (right) triplets for expanding to a
domain boundary.

######  Algorithm 2.26 (adapt).

In contrast to the two-dimensional case, we do not adaptively coarse the set
of triplets. While this is possible and does not harm accuracy, it complicates
computing a surface mesh from the set of triplets representing $\Gamma_{i,j}$.

For adaptive refinement, we do not rely on a global triangulation of
$\breve{S}_{i,j}$. Instead, for a given triplet $x\in\breve{S}_{i,j}$, let
$\breve{S}_{\text{loc}}$ consist of the $k_{\text{near}}$ nearest triplets in
$\breve{S}_{i,j}$ to $x$. We least-squares fit a plane $\mathcal{E}$ to
$\breve{S}_{i,j}^{\text{loc}}$ as in Algorithm 2.23 (fill) and create a
Delaunay triangulation $\mathcal{T}$ of $\breve{S}_{i,j}^{\text{loc}}$
projected to $\mathcal{E}$. For each non-degenerated triangle in
$\mathcal{T}$, we estimate the error applying esterror to
$\breve{S}_{\text{loc}}^{(i,j)}$. If the estimated error exceeds
$\varepsilon_{\text{err}}$, we employ the center of the triangle as a starting
point for adding a new triplet. Looping over all $x\in\breve{S}_{i,j}$
naturally leads to many duplicate or very close starting points which need to
eliminated before adding triples via bisection. However, these duplicates do
not harm the efficiency of our method, as no function evaluations are required
before computing triplets from starting points. We enrich $\breve{S}_{i,j}$
with the new triples created and repeat the refinement procedure at most
$k_{\text{adap}}$ times or until no more starting points for computing
triplets have been created.

######  Example 2.27.

We compute Test problem 2.24 starting with $200$ Halton-distributed points and
$k_{\text{near}}=10$. All other parameters are set as in the computation of
Test problem 2.4 in Section 2.1. For the corresponding numbers of triplets and
the number of function evaluations, we refer to Table 1. Figs. 13, 14 and 16
provide visualisations of the sets of triplets.

number of triplets | up to $\mathcal{M}^{2}$ | iniapprox | fill | expand | adapt
---|---|---|---|---|---
approx. of $\Gamma_{1,2}$ |  | 102 | 1459 | 1704 | 2460
approx. of $\Gamma_{1,3}$ |  | 83 | 1685 | 2058 | 2774
approx. of $\Gamma_{2,3}$ |  | 36 | 745 | 981 | 1510
function evaluations | 616 | 2157 | 25510 | 14232 | 5218
Table 1: Number of triplets per $\Gamma_{i,j}$ and number of function
evaluations for Example 2.27. Figure 16: Final sets $\hat{S}_{1,2}^{(1)}$ and
$\hat{S}_{1,3}^{(1)}$ for Test problem 2.24.

##  3 Direct and inverse acoustic scattering problem

In this section, we briefly explain the direct and inverse scattering for the
acoustic transmission problem. Most of the description can be found in [4] and
we refer the reader to [4] for more details.

Let the scatterer $D\subset\mathbb{R}^{3}$ be a given bounded domain with
boundary $\partial D$ assumed to be of class $C^{2,\alpha}$. The normal unit
vector $\nu$ points into the exterior $E=\mathbb{R}^{3}\backslash\overline{D}$
of the scatterer. The exterior $E$ assumed to be simply-connected is an
infinite homogeneous isotropic non-absorbing acoustic medium which is
characterised by the mass density $\varrho_{e}$, mean compressibility
$\kappa_{e}$, and sound speed $c_{e}=1/\sqrt{\kappa_{e}\varrho_{e}}$.
Likewise, the interior of $D$ is characterised by $\varrho_{i}$, $\kappa_{i}$,
and $c_{i}=1/\sqrt{\kappa_{i}\varrho_{i}}$. The given scatterer is excited by
a time-harmonic acoustic incident plane wave of the form

$\displaystyle
u^{\mathrm{inc}}(x;\widehat{d})=\mathrm{e}^{\mathrm{i}k_{e}x\cdotp\widehat{d}}\,,\qquad
x\in\mathbb{R}^{3}\,,$ (7)

where $k_{e}=\omega/c_{e}$ is the wave number of the acoustic wave in the host
medium, $\omega>0$ the angular frequency, and $\widehat{d}\in\mathbb{S}^{2}$
the direction of incidence with
$\mathbb{S}^{2}=\\{x\in\mathbb{R}^{3}:\|x\|=1\\}$ the unit sphere, where
$\|\cdotp\|$ denotes the standard Euclidean norm in $\mathbb{R}^{3}$. Note
that the incident field also depends on $k_{e}$, but this dependence is
suppressed.

The incident wave (7) interferes with the penetrable scatterer and creates two
waves. The first wave is the scattered field $u^{\mathrm{sca}}(x;\widehat{d})$
defined for $x\in E$ propagating outward and the second wave is the
transmitted field $u^{\mathrm{int}}(x;\widehat{d})$ defined for $x\in D$. The
total field in $E$ denoted by $u^{\mathrm{ext}}(x;\widehat{d})$ is the
superposition of $u^{\mathrm{int}}(x;\widehat{d})$ and
$u^{\mathrm{sca}}(x;\widehat{d})$ each of which satisfies the Helmholtz
equation (the reduced wave equation) in $E$ with wave number $k_{e}$.
Likewise, the transmitted field satisfies the Helmholtz equation in $D$ with
wave number $k_{i}$. Precisely, we have

$\displaystyle\Delta
u^{\mathrm{int}}(x;\widehat{d})+k_{i}^{2}u^{\mathrm{int}}(x;\widehat{d})=0\,,$
$\displaystyle x\in D\,,$ $\displaystyle\Delta
u^{\mathrm{ext}}(x;\widehat{d})+k_{e}^{2}u^{\mathrm{ext}}(x;\widehat{d})=0\,,$
$\displaystyle x\in E\,.$

Due to the continuity of the acoustic pressure and the normal component of the
particle velocity across $\partial D$ yields the transmission boundary
conditions

$u^{\mathrm{int}}(x;\widehat{d})=u^{\mathrm{ext}}(x;\widehat{d})\quad\text{
and
}\quad\partial_{\nu}u^{\mathrm{int}}(x;\widehat{d})=\tau\partial_{\nu}u^{\mathrm{ext}}(x;\widehat{d})\,,\qquad
x\in\partial D\,,$

where $\tau=\varrho_{i}/\varrho_{e}>0$ is the mass density ratio of the two
media. To ensure a well-posed boundary value problem, the scattered field
$u^{\mathrm{sca}}(x;\widehat{d})$ needs to satisfy the Sommerfeld radiation
condition

$\lim_{r\rightarrow\infty}r\left(\partial_{r}u^{\mathrm{sca}}(x;\widehat{d})-\mathrm{i}k_{e}u^{\mathrm{sca}}(x;\widehat{d})\right)=0$

with $r=\|x\|$. The classical acoustic transmission problem reads: find the
functions $u^{\mathrm{int}}(x;\widehat{d})\in C^{2}(D)\cap
C^{1}(\overline{D})$ and $u^{\mathrm{sca}}(x;\widehat{d}))\in C^{2}(E)\cap
C^{1}(E)$ satisfying

$\displaystyle\Delta
u^{\mathrm{int}}(x;\widehat{d})+k_{i}^{2}u^{\mathrm{int}}(x;\widehat{d})=0\,,$
$\displaystyle x\in D\,,$ (8) $\displaystyle\Delta
u^{\mathrm{ext}}(x;\widehat{d})+k_{e}^{2}u^{\mathrm{ext}}(x;\widehat{d})=0\,,$
$\displaystyle x\in E\,,$ (9) $\displaystyle
u^{\mathrm{int}}(x;\widehat{d})-u^{\mathrm{sca}}(x;\widehat{d})=u^{\mathrm{inc}}(x;\widehat{d})\,,$
$\displaystyle x\in\partial D\,,$ (10)
$\displaystyle\frac{1}{\tau}\partial_{\nu}u^{\mathrm{int}}(x;\widehat{d})-\partial_{\nu}u^{\mathrm{sca}}(x;\widehat{d})=\partial_{\nu}u^{\mathrm{inc}}(x;\widehat{d})\,,$
$\displaystyle x\in\partial D\,,$ (11)
$\displaystyle\lim_{r\rightarrow\infty}r\left(\partial_{r}u^{\mathrm{sca}}(x;\widehat{d})-\mathrm{i}k_{e}u^{\mathrm{sca}}(x;\widehat{d})\right)=0\,,$
$\displaystyle r=\|x\|$ (12)

###  3.1 The direct acoustic transmission problem

Given the incident field (i.e. the direction of incidence $\widehat{d}$ and
the wave number $k_{e}$), the scatterer $D$ (hence also its boundary $\partial
D$), the wave number $k_{i}$, and the parameter $\tau$, one has to solve (8) –
(12) for $u^{\mathrm{int}}(x;\widehat{d})$ and
$u^{\mathrm{sca}}(x;\widehat{d})$. In the direct problem, one is only
interested in the far-field $u^{\infty}(\widehat{x};\widehat{d})$ of
$u^{\mathrm{sca}}(x;\widehat{d})$ which is given by

$u^{\mathrm{sca}}(x;\widehat{d})=\frac{\mathrm{e}^{\mathrm{i}k_{e}\|x\|}}{\|x\|}u^{\infty}(\widehat{x};\widehat{d})+\mathcal{O}\left(\|x\|^{-2}\right)\,,\qquad\|x\|\rightarrow\infty$

uniformly in all directions $\widehat{x}\in\mathbb{S}^{2}$. The far-field can
be found by evaluating an integral equation over $\partial D$ given two
density functions determined by first solving a $2\times 2$ system of integral
equation of the second kind over $\partial D$ (see [4, Section 4.1]). Of
course, the integral equations at hand cannot be solved analytically and have
to be solved numerically for example by the boundary element collocation
method (see [21, Chapter 5] for more details).

To sum up, in the direct acoustic transmission problem one is interested in
$u^{\infty}(\widehat{x};\widehat{d})$ for $\widehat{x}\in\mathbb{S}^{2}$ given
the scatterer’s boundary $\partial D$ and the direction of incidence
$\widehat{d}\in\mathbb{S}^{2}$. The parameters $k_{e}$, $k_{i}$, and $\tau$
are given.

###  3.2 The inverse acoustic transmission problem

The parameters $k_{e}$, $k_{i}$, and $\tau$ are given. In the inverse acoustic
transmission problem one tries to find/reconstruct the domain’s boundary from
the knowledge of the far-field patterns $u^{\infty}(\widehat{x};\widehat{d})$
for all $\widehat{x},\widehat{d}\in\mathbb{S}^{2}$. This can be achieved with
the factorization method originally invented by Kirsch (see [20]). The
theoretical justification of the factorization method for the acoustic
transmission problem is given in [4, Chapter 3] and shown to work practically
in [4, Chapter 4]. We briefly outline the algorithm: Assume that the far-field
data are given for $\hat{x}_{i}$ and $\hat{d}_{j}$ with
$i,j\in\\{1,\ldots,m\\}$ stored in the matrix $A\in\mathbb{C}^{m\times m}$.
First, compute a singular decomposition of $A=U\Lambda V^{\ast}$ with
$\Lambda=\mathrm{diag}(\lambda_{1},\ldots,\lambda_{m})$. For a given point
$z\in\mathbb{R}^{3}$ compute the expansion coefficient of

$r_{z}=\left(\exp\left(-\mathrm{i}k_{e}z\cdotp\hat{d}_{j}\right)\right)_{j=1,\ldots,m}\in\mathbb{C}^{m}$

with respect to $V$ by

$\varrho_{\ell}^{(z)}=\sum_{j=1}^{M}V_{j,\ell}\mathrm{e}^{-\mathrm{i}k_{e}z\cdotp
d_{j}}\,,\quad\ell=1,\ldots,M\,,$

which is a matrix-vector multiplication $\varrho^{(z)}=V^{\top}r_{z}$.
Finally, we compute

$W(z)=\left[\sum_{\ell=1}^{M}\frac{|\varrho_{\ell}^{(z)}|^{2}}{|\lambda_{\ell}|}\right]^{-1}$

and plot the isosurfaces of $z\mapsto W(z)$. The values of $W(z)$ should be
much smaller for $z\notin D$ than those lying within $D$. The threshold value
can be approximately determined by a scatterer such as the unit sphere and
then reused for other scatterers as well. Note that until now, an equidistant
set of points $N$ within a predefined box have been used to find the values of
$W(z)$ leading to an amount of $N\times N\times N$ function evaluations for
such a “sampling” method (see also [8] for other sampling methods). This can
be considerably reduced as shown in the next section.

##  4 Applications

###  4.1 Test cases related to the detection of faults

In this section, we consider test cases which have been defined in the context
of the detection of fault lines or fault surfaces.

######  Test problem 4.1.

Following Allasia et al. [1] and Gutzmer et at. [16], we set
$x_{M}=(0.5,0.5)^{\top}$, $x=(x_{1},x_{2})^{\top}$, $\Omega=[0,1]^{2}$ and
consider the function

$f(x_{1},x_{2})=\left\\{\begin{array}[]{ccl}1+2\lfloor
3.5\|(x_{1},x_{2})\|_{2}\rfloor&,&\|(x-x_{M})\|_{2}<0.4\\\
0&,&\text{otherwise}\end{array}\right.$ (13)

This function is piecewise constant with smooth fault lines, it holds

$\displaystyle\Gamma=\bigcup\Gamma_{i,j}=$ $\displaystyle\hskip
9.95863pt\left\\{x\in\mathbb{R}^{2}\>:\>\|x-x_{m}\|_{2}<0.4\right\\}$
$\displaystyle\cup\left\\{\|x\|^{2}_{2}=4/7\>:\>\|x-x_{m}\|_{2}<0.4\>\right\\}$
$\displaystyle\cup\left\\{\|x\|^{2}_{2}=6/7\>:\|x-x_{m}\|_{2}<0.4\>\right\\}\>.$

The set $\Gamma_{1,3}$ consists of two separate components (Fig. 17, left).

######  Test problem 4.2.

We set

$f(x_{1},x_{2})=\left\\{\begin{array}[]{rcl}3&,&x_{1}>0.6\\\ 1&,&x_{1}<0.5\\\
2&,&\text{otherwise}\end{array}\right.\>.$ (14)

The set $\Gamma$ consists of two straight lines and coincides with the two
fault lines from [1], Example 4 (Fig. 17, middle).

######  Test problem 4.3.

We set

$f(x_{1},x_{2})=\left\\{\begin{array}[]{rcl}2&,&x_{1}>0.4\wedge
x_{2}>0.4\wedge x_{2}<0.2+x_{1}\\\ 1&,&\text{otherwise}\end{array}\right.$
(15)

The set $\Gamma_{1,2}$ coincides with the fault line from [1], Example 5 (Fig.
17, right).

Figure 17: Reconstructed subdivisions for Test problems 4.1, 4.2, and 4.3.

Figure 18: $S^{(i)}_{i,j}$ (left), $\breve{S}_{i,j}^{(i)}$ (middle), and
$\hat{S}_{i,j}^{(i)}$ (right) for Test problem 4.1. As for the computations of
Example 2.4 in Section 2.1, previously existing points are greyed out.

Figure 19: $S^{(i)}_{i,j}$ (left), $\breve{S}_{i,j}^{(i)}$ (middle), and
$\hat{S}_{i,j}^{(i)}$ (right) for Test problem 4.2.

Figure 20: $S^{(1)}_{1,2}$ (left), $\breve{S}_{1,2}^{(1)}$ (middle), and
$\hat{S}_{1,2}^{(1)}$ (right) for Test problem 4.3.

We compute all three test problems starting with $X$ from Example 2.4 and the
same parameters used there and in the subsequent examples in Section 2.1. Our
results (Figs. 18, 19, and 20) demonstrate a successful approximation of all
$\Gamma_{i,j}$. We provide details in Tables 2 and 3 combined with the results
for Test problem 2.3 from the preceding section. Algorithm iniapprox is the
most demanding Building block in terms of function evaluations in all examples
considered. Averaged over Test problems 2.3, 4.2, and 4.3, it takes on average
4.0 classifications for adding one triplet in fill. However, it takes 12.2
classifications per triplet added for Test problem 4.1. This is due to the
fact that $\Gamma_{1,3}$ consists of two components, which is detected by a
failed attempt of filling the large gap between these two components (Remark
2.15). This process adds classifications, but no triplets. In average over the
four examples, expand requires $5.5$ classifications per triplet added. A
significant part of the classifications is required for the very first and
very last step, as finding them requires bisection on the extrapolation curve
$\gamma$, which involves at least one classification per iteration step.
Considering the same ratio for adapt would be misleading, as triplets are
added and removed.

Allasia et al. present numbers of function evaluations for their method of
surface reconstruction applied to these test problems in [1]. They report
$5185$ evaluations of $f$ for Test problem 4.1 (our method: $2213$), $3479$
for Test Problem 4.2 (our method: $871$), and $2099$ for Test problem 4.3 (our
method: $450$). It turns out that the number of function evaluations is
significantly higher than for our algorithm. However, this comparison is
limited by several factors. In contrast to us, Allasia et al. need to employ
an explicit classification algorithm, as they consider more general functions
$f$ than we do, which may infer additional function evaluations. The number of
function evaluations required heavily depends on the accuracy of the
approximation and resolution of the fault line, such that comparing the number
of function evaluations given a prescribed tolerance for both algorithms would
provide a more sound basis for comparing the efficiency of the two algorithms.
Detailed information on accuracy is however lacking in Test problem 4.1.

Test P. | up to $\mathcal{M}^{2}$ | iniapprox | fill | expand | adapt
---|---|---|---|---|---
2.3 | 146 | 458 | 94 | 185 | 205
4.1 | 254 | 1017 | 465 | 350 | 127
4.2 | 176 | 595 | 26 | 74 | 0
4.3 | 87 | 202 | 34 | 44 | 83
Table 2: Number of function evaluations per Building block for reconstructing the subdomains. Test P. | iniapprox | fill | expand | adapt
---|---|---|---|---
2.3 | 40 | 55 | 89 | 125
4.1 | 66 | 98 | 138 | 171
4.2 | 32 | 45 | 60 | 6
4.3 | 17 | 26 | 36 | 18
Table 3: Total number of triplets after the respective Building block.

###  4.2 Decision analysis and modeling

MCDA methods depend on several parameters called “input factors”. Almost all
common MCDA models require to set up a _performance matrix_ $P=(p_{i,j})$,
where $p_{i,j}$ encodes the objective benefit or cost of the $i$-th
alternative with respect to the $j$-th criterion. Both finding suitable
criteria for modeling the decision process and setting up $P$ requires
expertise (e.g. [28] among many others) and is beyond the scope of this work.
Instead, we assume that $P$ is given and exactly known, although this
assumption may be questioned in many practical applications. In addition to
$P$, the decision maker needs to provide non-negative _weights_
$w=(w_{1},\ldots w_{c})^{\top}$ which reflect the importance of the criteria
from his point of view. Although originally intended as a decision support
tool, MCDA methods have been used for some time for decision analysis, e.g. to
predict decisions of actors under changed framework conditions, modelled by a
changed performance matrix. In the context of decision analysis, $w$ is in
many practical applications not exactly known and hard to obtain, as e.g.
surveys are time-consuming and prone to bias due to socially accepted answers
and other effects. Moreover, the restriction to fixed weights ignores possible
diversity within actors. All of the above motivates a robustness analysis of
the decision predicted, which in many cases focuses on examining the
robustness of the decision with respect to perturbations in $w$.

In what follows, we apply our algorithm for computing $\Gamma_{i,j}$ to such a
robustness analysis and consider as example an application from the scientific
monitoring of the mobility turnaround in Germany. Ball et al. [7] investigate
car users’ attitudes towards the purchase of hybrid (HEV) and electric
vehicles (BEV) versus conventional cars with internal combustion engine (ICE).
For this purpose, they identify 13 criteria and divide them into 5 categories
(“Ecological”,“Economic”, “Social”, “Comfort”, “Other”). The authors weight
both the criteria within the categories and the categories themselves. Taking
the former weights as given, we obtain the performance matrix in Table 4 from
the data in [7]. In contrast to SAW (Simple Additive Weighting) [11] as in
[7], we use the more complex MCDA method SIR-TOPSIS [30] here, which includes
the very widespread MCDA methods Promethee II [19] and SAW as special cases.
We omit all details on how TOPSIS works for the sake of brevity and refer
instead to [30].

| Ecological | Economic | Social | Comfort | Other
---|---|---|---|---|---
BEV | 0.5025 | 0.2792 | 0.6250 | 0.1497 | 0.1342
ICE | 0.1256 | 0.4167 | 0.1250 | 0.4300 | 0.6710
HEV | 0.3719 | 0.3042 | 0.2500 | 0.4202 | 0.1948
weightings | 2 | 7 | 0.1 | 8 | 0
Table 4: Performance matrix $\tilde{P}$ with corresponding criteria generated
from [7], and weightings from the car users’ point of view prior to
normalisation (last row).

SIR-TOPSIS requires as many other MCDA methods that the non-negative weights
are normalized, i.e. $\sum_{i=1}^{n}w_{i}=1$. Therefore, the set of normalized
admissible weights is the standard simplex in $\mathbb{R}^{n}$. For
visualisation, we consider $m=3$ or $m=4$ of the weights to be variable, and
the rest to be fixed. Let be $w=w_{v}+w_{f}$, where $w_{v}$ consists of the
variable weights and zeroes elsewhere, and $w_{f}$ of the fixed ones,
correspondingly. Therefore, normalisation implies $\sum w_{i,v}=1-\sum
w_{i,f}:=c_{f}$ such that the set of variable weights corresponds to a
downscaled standard simplex in $\mathbb{R}^{m}$, which we embed in
$\mathbb{R}^{m-1}$ by appropriate translation and rotation. This yields the
equilateral triangle (for $m=3$) and the regular tetrahedron (for $m=4$) shown
in Fig. 21.

For a 2D-visualisation of the decision space, we consider the weights
corresponding to “Ecological”, “Economic” and “Comfort” variable (Fig. 21) and
colour all weights leading to a decision in favour of ICEs black, for HEVs
blue and for BEVs green. We divide the weights $\tilde{w}$ proposed in [7] as
a representation of car user’s mindset in 2020 in a fixed and variable part,
setting $\tilde{w}=\tilde{w}_{v}+\tilde{w}_{f}$ and display $\tilde{w}_{v}$ as
bright red dot in Fig. 21. In contrast to [7], SIR-TOPSIS seems to predict a
shift to HEVs under today’s conditions. For measuring robustness, we consider
the largest sphere around $\tilde{w}_{v}$ still fully contained in the set of
weightings leading to HEVs $W_{\text{HEV}}$ and propose its radius $\rho$ as
simple measure of robustness. As we have a polygonal approximation of
$W_{\mathrm{HEV}}$ at hand, iteratively approximating $\rho$ boils down to an
intersection test of polygons, if we approximate the circle by a sufficiently
fine polygon. We end up with $\rho\approx 0.06$, which reconciles the findings
of Ball et al. and ours: As $\rho$ is rather small, the decision in favour of
HEVs is not very robust, as even small perturbations of the weightings may
lead to a decision in favour of ICEs. In [7], the decision towards ICEs was
found not be very robust. However, both ours and Ball et al.’s results
indicate that a broad shift towards BEVs is unlikely to happen under current
circumstances.

Figure 21: Visualisation of the decision space for $3$ (left) and $4$ (right)
variable criteria. The car users’ weightings according to [7] are displayed as
red dot along with the circumsphere with radius $\rho=0.06$.

For $m=3$ and $m=4$, we have $c_{f}=0.9942$. The downscaled standard simplex
is rotated and translated to the equilateral triangle with vertices
$c_{f}(0.4082,-0.7071)^{\top}$, $c_{f}(0.4082,0.7071)$, and
$c_{f}(-0.8165,0)$. Therefore, we set $\Omega=c_{f}[-0.9,0.4082]\times
c_{f}[-0.9,0.9]$ and employ an initial point set $X$ consisting of $100$
Halton-distributed points, where points far away from the triangle have been
discarded, as they cannot aid approximating $\Gamma_{i,j}$ (Fig. 22, left). We
choose the values of all algorithm-related parameters as for the examples in
Sections 2.1 and 4.1, except of $\varepsilon_{\mathrm{gap}}=c_{f}\cdot 0.05$.
For $m=4$, we start with $X$ consisting of $500$ Halton-distributed points in
the vicinity of the tetrahedron and take all values for the parameters of the
algorithm from Section 2.2, except of $k_{\text{adap}}=3$. For the final sets
$\hat{S}_{i,j}$, we refer to Fig. 22 (right and middle) and display the number
of function evaluations in Table 5. These sets have been used for
approximating $\rho$; Fig. 21 was generated with the proposed algorithm,
albeit using a finer initial point set and modified parameter settings, as
analytical descriptions of the decision curves and surfaces, resp., are
unknown.

Figure 22: Initial point set $X$ (left) and final sets $\hat{S}_{i,j}$ (middle) for $m=3$. Right: Selected final sets $\hat{S}_{i,j}$ for $m=4$. | up to $\mathcal{M}^{2}$ | iniapprox | fill | expand | adapt
---|---|---|---|---|---
evaluations, $m=3$ | 397 | 1325 | 198 | 283 | 44
evaluations, $m=4$ | 1518 | 4553 | 13352 | 14752 | 24344
no. of triplets, $m=3$ |  | 100 | 145 | 174 | 45
no. of triplets, $m=4$ |  | 335 | 1575 | 2401 | 6999
Table 5: Number of function evaluations per Building block for reconstructing
the decision boundaries (see Section 4.2).

Due to recent geopolitical events, car users today may attach a different
importance to issues of security of supply, e.g. with fuel, than in 2021 when
[7] was written. That kind of considerations are subsumed in the category
“Social” which motivates to additionally consider the weights of that category
to be variable. It turns out that with stronger emphasis on the category
“Social” car users tend to prefer a BEV (Fig. 21, right). One reason for this
may be that the dependence on oil imports makes the purchase of an ICE or even
an HEV seem less attractive.

###  4.3 Surface Reconstruction in 3D from scattering

In our tests, we will use different scatterers to apply our algorithm to such
as the unit sphere, the ellipsoid, the peanut, the acorn, the cushion, the
round short cylinder, and the round long cylinder. Their surfaces are given
parametrically in spherical coordinates $x=r_{1}\sin(\phi)\cos(\theta)$,
$y=r_{2}\sin(\phi)\cos(\theta)$, and $z=r_{3}\cos(\phi)$ with
$\theta\in[0,2\pi)$, $\theta\in[0,\pi]$ as $r_{1}=r_{2}=r_{3}=1$ for the unit
sphere, $r_{1}=r_{2}=1$ and $r_{3}=6/5$ for the ellipsoid,
$r_{1}^{2}=r_{2}^{2}=r_{3}^{2}=9(\cos^{2}(\phi)+\sin^{2}(\phi)/4)/4$ for the
peanut, $r_{1}^{2}=r_{2}^{2}=r_{3}^{2}=9(17/4+2\cos(3\phi))/25$ for the acorn,
$r_{1}=r_{2}=r_{3}=1-\cos(2\phi)/2$ for the cushion,
$r_{1}^{10}=r_{2}^{10}=r_{3}^{10}=1/((2\sin(\phi)/3)^{10}+\cos^{10}(\phi))$
for the round short cylinder, and
$r_{1}^{10}=r_{2}^{10}=r_{3}^{10}=1/((2\cos(\phi)/3)^{10}+\sin^{10}(\phi))$
for the round long cylinder, respectively. We will use $m=1026$ number of
incident and observation directions for the construction of the far-field data
with the parameters $k_{e}=2$, $k_{i}=1$, and $\tau=1/2$ as also used in [4,
p. 18]. Therefore, the factorization algorithm appears as classification
function $f$.

For our experiments, we set all algorithm-related parameters as in Section 2.2
except of $\varepsilon_{\mathrm{err}}=0.01$ and
$\varepsilon_{\text{gap}}=0.25$. The initial set $X$ consists of $200$ Halton-
distributed points in $\Omega=[-1.5,1.5]^{3}$ apart of the long and the short
cylinder, where we set $\Omega=[-2,2]^{3}$. The results (Figs. 23 and 24)
indicate successful reconstructions. For the number of triplets after each
part of the algorithm and the number of function evaluations, we refer to
Tables 6 and 7. As no triplet in these tests fulfilled condition (6), no
expansion took place, such that we omit the corresponding column in our
tables. As all scatters feature one closed surface, this behaviour of our
algorithm is as desired. In order to obtain a visually appealing
reconstruction of the scatterers, Anagnostopoulos et al. [4] create a tensor
product set of $55^{3}=166,375$ points and use the corresponding
classifications to compute isosurfaces based upon these data for visualisation
purposes. Our approach, on the other hand, requires only a fraction of these
function evaluations. Since about half of the total computation time was used
by Kirsch’s factorization method in the surface reconstruction with a Matlab
implementation of our method, our algorithm enables a significant speedup
compared to [4].

Moreover, in contrast to the level set approach, our method provides a set of
points near the scatterer, which can be used for computing a further refined
surface representation like triangulation or higher order interpolation. There
are two major sources of inaccuracy in surface reconstruction: the inaccuracy
of the factorization method itself and the error induced by the representation
of the reconstructed surface, be it an isosurface or a set of points. As our
algorithm controls the latter one, comparing with the ground truth aka the
analytical description of the true surface allows for analysing the former
error source far more easily than having only implicit surfaces available.

######  Remark 4.4.

In our tests, we do not commit to inverse crime. The far field data used for
reconstruction later one have been produced using boundary element collocation
with high accuracy. The inverse problem is solved using Kirsch’s factorization
method combined with the algorithm proposed in this work.

Figure 23: Reconstructed scatterers acorn, bumpy sphere, cushion and
ellipsoid.

Figure 24: Reconstructed Scatterers long cylinder, peanut, short cylinder, and sphere. scatterer | up to $\mathcal{M}^{2}$ | iniapprox | fill | adapt
---|---|---|---|---
acorn | 463 | 1891 | 7242 | 8272
bumpy sphere | 384 | 1326 | 4844 | 8911
cushion | 473 | 1956 | 6237 | 7881
ellipsoid | 383 | 1329 | 4694 | 5458
long cyl. | 378 | 1377 | 7119 | 9816
peanut | 394 | 1381 | 5149 | 6922
short cyl. | 465 | 2033 | 9021 | 9667
sphere | 384 | 1314 | 3954 | 5102
Table 6: Number of function evaluations per Building block for reconstructing the scatterers. scatterer | iniapprox | fill | adapt
---|---|---|---
acorn | 222 | 1291 | 3839
bumpy sphere | 158 | 827 | 3103
cushion | 232 | 1125 | 3613
ellipsoid | 157 | 843 | 2305
long cyl. | 149 | 1188 | 4040
peanut | 164 | 878 | 2864
short cyl. | 217 | 1594 | 4543
sphere | 159 | 715 | 2120
Table 7: Total number of triplets after the respective Building block for
reconstructing the scatterers.

##  5 Conclusions and Outlook

In this article, we presented a method for approximating manifolds of
discontinuity of a function $f$ in 2D and 3D and demonstrated successful
applications of our method to Multi-criteria Decision Aid, to generic test
cases connected with the detection of faults and to an inverse acoustic
scattering problem. In all cases, our method requires significantly fewer
evaluations of $f$ than previously existing algorithms we compared our method
to. At least for the inverse acoustic scattering problem, this leads to a
significant acceleration of the overall reconstruction, as computing the
factorization method dominates the computational time even in our
computations, where the number of such computations could be reduced by a
factor of approx. 7 to 12.

Our algorithm could be easily enhanced with classification algorithms as in
[1] in order to consider more general functions than we did, provided the
classification is certain. This enables for tackling fault detection problems
and for applying our method for interpolation of piecewise smooth functions
assuming that $f$ is smooth on $\Omega_{1},\ldots,\Omega_{n}$ with
$\Omega=\overline{\Omega_{1}}\cup\ldots\cup\overline{\Omega_{n}}$, but
globally discontinuous. In [23], the authors propose interpolation with radial
basis functions on each $\Omega_{i}$ in this setting. This however requires
knowledge about the boundaries of each $\Omega_{i}$ which could be obtained
using our method. For that purpose, [23] provides a fault detection algorithm
based on local approximation properties. Combining these two approaches is
subject of our current research.

Acknowledgements: The authors want to thank Rodin Eybesh and Luis Hasenauer
for porting a significant part of our Matlab implementation to python.

## Declarations

  * •

Conflict of Interest: The authors declare no competing interests.

  * •

Code Availability: The codes used for producing the results presented are
available at https://github.com/mgrajewski/faultapprox-matlab and
https://github.com/mgrajewski/faultapprox-python

## References

  * [1] Allasia, G., Besenghi, R., Cavoretto, R., and De Rossi, A. Efficient approximation algorithms. Part I: approximation of unknown fault lines from scattered data.  Dolomites Research Notes on Approximation 3 (2010), 7–38.
  * [2] Althaus, E., and Mehlhorn, K. Traveling salesman-based curve reconstruction in polynomial time.  SIAM Journal on Computing 31, 1 (2001), 27–66.
  * [3] Amenta, N., Bern, M., and Eppstein, D. The crust and the -skeleton: Combinatorial curve reconstruction.  Graphical Models and Image Processing 60, 2 (1998), 125–135.
  * [4] Anagnostopoulos, K. A., Charalambopoulos, A., and Kleefeld, A. The factorization method for the acoustic transmission problem.  Inverse Problems 29, 11 (2013), 115015.
  * [5] Arge, E., and Floater, M. Approximating scattered data with discontinuities.  Numerical Algorithms 8 (1994), 149–166.
  * [6] Audibert, L., and Haddar, H. A generalized formulation of the linear sampling method with exact characterization of targets in terms of farfield measurements.  Inverse Problems 30, 3 (2014), 035011.
  * [7] Ball, C. S., Vögele, S., Grajewski, M., and Kuckshinrichs, W. E-mobility from a multi-actor point of view: Uncertainties and their impacts.  Technological Forecasting and Social Change 170 (2021), 120925\.
  * [8] Bazán, F. S. V., Kleefeld, A., Leem, K. H., and Pelekanos, G. Sampling method based projection approach for the reconstruction of 3D acoustically penetrable scatterers.  Linear Algebra and its Applications 495 (2016), 289–323.
  * [9] Bishop, C. M. Pattern Recognition and Machine Learning.  Springer, New York, 2006.
  * [10] Bozzini, M., and Rossini, M. The detection and recovery of discontinuity curves from scattered data.  Journal of Computational and Applied Mathematics, 240 (2013), 148–162.
  * [11] Churchman, C. W., and Ackoff, R. L. An approximate measure of value.  Journal of the Operations Research Society of America 2, 2 (1954), 172–187.
  * [12] Colton, D., and Kirsch, A. A simple method for solving inverse scattering problems in the resonance region.  Inverse Problems 12, 4 (1996), 383–393.
  * [13] Dey, T. K., Mehlhorn, K., and Ramos, E. A. Curve reconstruction: Connecting dots with good reason.  Computational Geometry 15, 4 (2000), 229–244.
  * [14] Figueira, J. Multiple Criteria Decision Analysis: State of the Art Surveys, vol. 78 of SpringerLink Bücher.  Springer New York, 2005.
  * [15] Gout, C., Le Guyader, C., Romani, L., and Saint-Guirons, A.-G. Approximation of surfaces with fault(s) and/or rapidly varying data, using a segmentation process, ${D}^{m}$-splines and the finite element method.  Numerical Algorithms 48 (2008), 67–92.
  * [16] Gutzmer, T., and Iske, A. Detection of discontinuities in scattered data approximation.  Numerical Algorithms 16 (1997), 155–170.
  * [17] Harris, I., and Kleefeld, A. Analysis of new direct sampling indicators for far-field measurements.  Inverse Problems 35, 5 (2019), 054002.
  * [18] Ikehata, M. The probe method and its applications.  In Inverse Problems and Related Topics (London, 2000), G. Nakamura, S. Saitoh, J. Seo, and M. Yamamoto, Eds., vol. 419 of Research Notes in Mathematics, CRC Press.
  * [19] J.P. Brans, Ph. Vincke, and B. Mareschal.  How to select and how to rank projects: The PROMETHEE method.  European Journal of Operational Research 24 (1986), 228–238.
  * [20] Kirsch, A., and Grinberg, N. The Factorization Method for Inverse Problems.  Oxford University Press, Oxford, 2008.
  * [21] Kleefeld, A. The transmission problem for the Helmholtz equation in $\mathbb{R}^{3}$.  Comput. Methods Appl. Math. 12, 3 (2012), 330–350.
  * [22] Kleefeld, A., and Lin, T.-C. The nonlinear Landweber method applied to an inverse scattering problem for sound-soft obstacles in 3D.  Computer Physics Communications 182, 12 (2011), 2550–2560.
  * [23] Lenarduzzi, L., and Schaback, R. Kernel-based adaptive approximation of functions with discontinuities.  Applied Mathematics and Computation 307 (2017), 113–123.
  * [24] Ohrhallinger, S., and Mudur, S. An efficient algorithm for determining an aesthetic shape connecting unorganized 2D points.  Computer Graphics Forum 32, 8 (2013), 72–88.
  * [25] Papathanasiou, J., and Nikolaos, P. Multiple Criteria Decision Aid: Methods, Examples and Python Implementations.  Springer, 2019.
  * [26] Potthast, R. A survey on sampling and probe methods for inverse problems.  Inverse Problems 22, 2 (2006), R1.
  * [27] Shakarji, C. M. Least-squares fitting algorithms of the NIST algorithm testing system.  Journal of Research of the National Institute of Standards and Technology 103, 6 (1998), 633–641.
  * [28] Vögele, S., Ball, C., and Kuckshinrichs, W. Multi-criteria approaches to ancillary effects: The example of E-mobility.  In Ancillary Benefits of Climate Policy: New Theoretical Developments and Empirical Findings (Cham, 2020), W. Buchholz, A. Markandya, D. Rübbelke, and S. Vögele, Eds., Springer International Publishing, pp. 157–178.
  * [29] Waldron, S. The error in linear interpolation at the vertices of a simplex.  SIAM Journal on Numerical Analysis 35, 3 (1998), 1191–1200.
  * [30] Xu, X. The SIR method: A superiority and inferiority ranking method for multiple criteria decision making.  European Journal of Operational Research 131, 3 (2001), 587–602.
  * [31] Zeng, F., Suarez, P., and Sun, J. A decomposition method for an interior inverse scattering problem.  Inverse Problems and Imaging 7, 1 (2013), 291–303.