# Degree of Convexity and Expected Distances in Polygons

Mikkel Abrahamsen∗ Viktor Fredslund-Hansen∗

(October 2022)

###### Abstract

We present an algorithm for computing the so-called _Beer-index_ of a polygon
$P$ in $O(n^{2})$ time, where $n$ is the number of corners. The polygon $P$
may have holes. The Beer-index is the probability that two points chosen
independently and uniformly at random in $P$ can see each other. Given a
finite set $M$ of $m$ points in a simple polygon $P$, we also show how the
number of pairs in $M$ that see each other can be computed in $O(n\log
n+m^{4/3}\log^{\alpha}m\log n)$ time, where $\alpha<1.78$ is a constant.

We likewise study the problem of computing the expected geodesic distance
between two points chosen independently and uniformly at random in a simple
polygon $P$. We show how the expected $L_{1}$-distance can be computed in
optimal $O(n)$ time by a conceptually very simple algorithm. We then describe
an algorithm that outputs a closed-form expression for the expected
$L_{2}$-distance in $O(n^{2})$ time.

†† ∗ Department of Computer Science, University of Copenhagen,
<EMAIL_ADDRESS>Mikkel Abrahamsen is supported by Starting Grant
1054-00032B from the Independent Research Fund Denmark under the Sapere Aude
research career programme. Viktor Fredslund-Hansen is supported by Starting
Grant 7027-00050B from the Independent Research Fund Denmark under the Sapere
Aude research career programme. Both authors are part of Basic Algorithms
Research Copenhagen (BARC), supported by the VILLUM Foundation grant 16582.

##  1 Introduction

In this paper, we consider a function defined for pairs of points in a polygon
$P$, and we want to compute the expected value of the function when choosing
two points uniformly and independently at random from some set $M\subset P$
(in the sequel just called “two random points in $M$,” for brevity). Here, $M$
can either be a given finite set of points in $P$ or all points in $P$, i.e.,
$M=P$. For a pair of points $p,q\in P$, the functions we will consider are:

  * •

the indicator $[pq\subset P]\in\\{0,1\\}$, i.e., whether $p$ and $q$ can see
each other in $P$, and

  * •

the length of the geodesic shortest path from $p$ to $q$ in $P$ in the
$L_{1}$\- or $L_{2}$-metric.

Table 1 summarizes our results.

#### Beer-index.

With $M=P$ and the first function $[pq\subset P]$, we obtain the probability
$B(P)$ that two random points $p,q$ in $P$ see each other,

$B(P):=\frac{1}{|P|^{2}}\iint_{p\in P}\iint_{q\in P}[pq\subset
P]\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}p,$

where $|\cdot|$ denotes the area. Throughout the paper, our integrals are
defined with respect to Lebesgue measure. The measure $B(\cdot)$ is known as
the _Beer-index_ after Beer [4].

One of the well-known characterizations of a polygon $P$ being convex is that
all pairs of points in $P$ can see each other. Thus, the number $B(P)\in[0,1]$
quantifies the degree to which this holds, so it is a natural measure for the
degree of convexity of $P$. The problem of partitioning a polygon into
components that are close to convex is useful in many practical settings, as
such partitions provide similar benefits as convex partitions, while the
number of components can be significantly smaller when the components are
allowed to be slightly non-convex [16, 23]. This motivates ways to quantify
the degree of convexity and algorithms for computing these measures. Stern
[28] also described the measure $B(P)$, using the equivalent form

$B(P)=\frac{1}{|P|^{2}}\iint_{p\in
P}\left|\textrm{vis}(p)\right|\mathop{}\\!\mathrm{d}p,$

where $\textrm{vis}(p)=\textrm{vis}_{P}(p)$ is the _visibility polygon_ of the
point $p$, i.e., the region of points in $P$ that $p$ can see. Stern suggested
that $B(P)$ be estimated by measuring $|\textrm{vis}(p)|$ for every grid point
$p$ in $P$ from a rectangular grid. Rote [25] showed that $B(P)$ can be
expressed as a sum of $O(n^{9})$ closed-form expressions by an algorithm that
uses the same amount of time. Here, $P$ is a polygon which may have holes, and
$n$ is the number of corners. Using a bound from [21], Rote [26] later
observed that the number is in fact bounded by $O(n^{7})$.

Buchin, Kostitsyna, Löffler and Silveira [9] described another algorithm that
outputs $B(P)$, again for a given polygon $P$ which may have holes. They
claimed the running time of the algorithm to be $O(n^{2})$. Unfortunately, as
the algorithm is described in the journal paper [9], it may output $B(P)$ as a
sum of $\Omega(n^{4})$ closed-form expressions and likewise have a running
time of $\Omega(n^{4})$. In Appendix A, we show that even when $P$ is a simple
polygon, i.e., when $P$ has no holes, the running time of the algorithm may be
$\Omega(n^{4})$. We have made the authors aware of this mistake, and they have
informally described a different algorithm to us in private communication
which seems to have a running time of $O(n^{2}\log n)$. However, at the time
of writing, the algorithm has apparently not been written down in detail.

For a given polygon $P$ with $n$ corners and possibly with holes, we describe
an algorithm outputting $B(P)$ as a sum of $O(n^{2})$ closed-form expressions
using $O(n^{2})$ time. The algorithm performs a sort of rotational sweep
around each edge $e$. In a usual rotational sweep, a single ray is rotating
around a point, but in our case, two parallel rays anchored at the endpoints
of $e$ are rotating, and we keep track of everything “seen” from $e$ in the
strip between these rays. Computing the Beer-index $B(P)$ then boils down to
computing the contribution from pairs of points in certain trapezoids with a
parallel pair of rotating edges, and there are $O(n^{2})$ of these trapezoids.
The algorithm is arguably simpler than the algorithm from [9], which relies on
a transformation using geometric duality and the computation of overlays of
polygons in the dual plane, whereas our algorithm deals exclusively with
primitive objects in the usual primal space.

######  Theorem 1.

Given a polygon $P$ which may have holes, there is an algorithm that returns
the Beer-index $B(P)$ as a sum of $O(n^{2})$ closed-form expressions, each of
size $O(1)$, where $n$ is the number of corners of $P$. The algorithm uses
$O(n^{2})$ time and space.

Problem | Time | Space | Reference
---|---|---|---
Beer-index (pol. with holes) | $O(n^{2})$ | $O(n^{2})$ | Section 2
No. of visible pairs (points in simple pol.) | $O(n\log n+m^{4/3}\log^{\alpha}m\log n)$ | $O(n+m^{4/3})$ | Section 3
Total $L_{1}$-distance (points in $\mathbb{R}^{d}$) | $O(d\cdot m\log m)$ | $O(m)$ | Section 4.1
Expected $L_{1}$-distance (simple pol.) | $O(n)$ | $O(n)$ | Section 4.2
Expected $L_{2}$-distance (simple pol.) | $O(n^{2})$ | $O(n^{2})$ | Section 5
Table 1: Our results. Here, $n$ denotes the number of corners of a polygon,
$m$ denotes the number of discrete given points and $\alpha<1.78$ is a
constant.

#### Counting visible pairs of points.

Very recently, Buchin, Custers, van der Hoog, Löffler, Popov, Roeloffzen and
Staals [8] studied the problem of computing the number of visible pairs among
$m$ points in a simple polygon $P$ with $n$ corners. If $m$ is small compared
to $n$, they suggest to use a data structure by Hershberger and Suri [19] to
test if each pair is visible in $O(\log n)$ time. This results in an algorithm
with running time $O(n+m^{2}\log n)$. For the case that $m$ is large, they
describe an algorithm with running time $O(n+m^{3/2+\varepsilon}\log
n\log(nm))$ for any constant $\varepsilon>0$. We present an algorithm for
counting the number of visible pairs among $m$ points in time $O(n\log
n+m^{4/3}\log^{\alpha}m\log n)$. Our algorithm is therefore polynomially
faster than the algorithms from [8] except if $m$ is small, i.e.,
$m=o(n^{2/3})$. We furthermore believe that our algorithm serves as a simple
and instructive demonstration of two important concepts in computational
geometry: (i) the use of divide-and-conquer by splitting a simple polygon
along a diagonal into two parts with nearly equally many corners, and (ii) the
transformation of one problem to an apparently more tractable problem using
geometric duality.

######  Theorem 2.

Given a simple polygon $P$ with $n$ corners and a set $M$ of $m$ points in
$P$, the number of pairs of points in $M$ that can see each other can be
computed in $O(n\log n+m^{4/3}\log^{\alpha}m\log n)$ time, where $\alpha<1.78$
is a constant.

Finding the edges of the visibility graph of a given polygon is one of the
classical problems in computational geometry. Lee [22] described a simple and
well-known algorithm with running time $O(n^{2}\log n)$ already in 1979, where
$n$ is the number of corners. The algorithm works by performing a rotational
sweep around all corners. The same algorithm was described by Sharir and
Schorr [27], and it also appears in the book [5]. Asano, Asano, Guibas,
Hershberger and Imai [2] and Welzl [30] gave algorithms with running time
$O(n^{2})$. Hershberger [18] gave an output-sensitive algorithm using $O(k)$
time to compute the visibility graph of a triangulated simple polygon, where
$k$ is the number of edges in the visibility graph. Together with Chazelle’s
algorithm [12] for triangulating a simple polygon in $O(n)$ time, this yields
an optimal algorithm with time $O(k)$ (note that $k=\Omega(n)$). When the set
$M$ of points in $P$ are the corners of $P$, our algorithm from Theorem 2
returns the _size_ of the visibility graph of a simple polygon in
$O(n^{4/3}\log^{\alpha+1}n)$ time. The number $k/\binom{n}{2}\in[0,1]$ can be
considered a discrete variant of the Beer-index. Note that using the
algorithms from [8] result in running time $\omega(n^{3/2})$.

######  Corollary 3.

Given a simple polygon $P$, the number $k$ of edges in the visibility graph of
$P$ can be computed in $O(n^{4/3}\log^{\alpha+1}n)$ time.

#### Expected distances.

The problem of determining the expected distance between two random points
from a given domain has a long history. Czuber’s book from 1884 [13] contains
calculations of the values for equilateral triangles, squares and rectangles.
Bäsel [10] recently derived formulas for the expected distance, as well as
higher moments, in regular $n$-gons for $n=3,4,5,6,8,10,12$. The paper
likewise contains more historical information about these problems. In another
recent paper, Bonnet, Gusakova, Thäle and Zaporozhets [6] proved that the
expected distance between two points in a convex body in the plane with
perimeter $1$ is between $7/60$ and $1/6$, and that these bounds are tight.
They also provide bounds for higher dimensional convex bodies. The problem of
computing expected distances often appears in distribution management,
transportation system analysis and airport terminal operations analysis and
modeling [20, 3, 29].

As a warmup for computing expected distances in simple polygons, we consider
the conceptually simpler problem of computing the sum of $L_{1}$-distances
between all pairs of points from a finite set $M\subset\mathbb{R}^{d}$, i.e.,
with no polygon involved. We give a very simple algorithm computing the sum in
$O(d\cdot m\log m)$ time, where $m=|M|$, which nicely demonstrates the
technique underlying our algorithm for computing the expected $L_{1}$-distance
in a simple polygon.

######  Theorem 4.

Given a set $M$ of $m$ points in $\mathbb{R}^{d}$, the sum of pairwise
distances in the $L_{1}$-metric between points in $M$ can be computed in
$O(d\cdot m\log m)$ time.

Hsu [20] studied the algorithmic problems of computing the expected distance
between two random points in a given polygon. He gave a $O(n^{2})$-time
algorithm for the expected geodesic $L_{1}$-distance in a simple polygon and a
$O(n^{3})$-time algorithm for the expected $L_{2}$-distance in a _convex_
polygon. In this paper, we describe a very simple algorithm for computing the
expected geodesic $L_{1}$-distance in optimal $O(n)$ time and another
algorithm outputting an expression for the $L_{2}$-distance in $O(n^{2})$
time, both in a simple polygon.

######  Theorem 5.

Given a simple polygon $P$ with $n$ corners, there is an algorithm for
computing the expected geodesic $L_{1}$-distance between two random points in
$P$ using $O(n)$ time and space.

In the case of the expected $L_{2}$-metric, we give a computer-assisted proof
that a closed-form expression of the expected distance can be computed in
$O(n^{2})$ time. The computer-assistance consists of finding closed-form
expressions for certain integrals using Maple [24]. Appendix C shows an
example of a Maple document for calculating a closed-form expression for an
integral (in this case, the integral expresses the contribution of one
trapezoid to the Beer-index, but expressions for the integrals arising when
computing the expected $L_{2}$-distance are obtained by slight variations of
this code). To the best of our knowledge, we describe the first algorithm that
computes the expected $L_{2}$-distance between two points in a given simple
polygon that may not be convex.

######  Theorem 6.

Given a simple polygon $P$ with $n$ corners, there is an algorithm that as an
output defines $O(n^{2})$ constants using closed-form expressions of size
$O(1)$, where the names of previously defined constants appear in the
expressions for subsequent ones. The algorithm then outputs the expected
geodesic $L_{2}$-distance in $P$ as a sum of a subset of these constants. The
algorithm uses $O(n^{2})$ time and space.

Let us remark that if we want to output the expected distance as one single
closed-form expression (i.e., not relying on predefined constants), the length
of the expression would instead be $O(n^{3})$.

The algorithm works by separately computing the contribution from (i) pairs of
points in $P$ that see each other, and (ii) pairs where the shortest path
contains a corner of $P$. For group (i), an algorithm identical to the one for
the Beer-index is used, except that we are summing over a different integral
for each trapezoid. For group (ii), we compute for each corner $u$ of $P$ the
contribution from pairs of points $p,q\in P$ where the first corner of $P$ on
the shortest path from $p$ to $q$ is $u$. We show that this boils down to
computing certain areas and integrals of distances from $u$ to all points on
one side of a diagonal $uv$, as well as two types of integrals over points in
triangles.

More prior attention has been given to the problem of computing the average
distance between two vertices in a graph $G$. This known as the _Wiener-index_
of $G$, and it is a fundamental measure with important applications in
mathematical chemistry and appears in thousands of publications. Note that the
Wiener-index is equivalent to the sum of pairwise distances. For more
information on the problem, see the papers [11, 15].

##  2 Beer-index

We now present an algorithm for computing the Beer-index $B(P)$ in a polygon
$P$ which may have holes. Assume without loss of generality that $P$ has area
$1$ and let $n$ be the number of corners of $P$. We furthermore assume that no
three corners of $P$ are on a line. This assumption can be avoided using a
symbolic perturbation of the points [14].

Figure 1: Left: The pair $p_{0},q_{0}$ contributes to $e_{0}$, since the
extension $p^{\prime}_{0}$ from $p_{0}$ is on $e_{0}$. The pair $p_{1},q_{1}$
contributes to $e_{1}$. Middle: The polygon $\textrm{vis}(e,\varphi)$ shown in
green with a trapezoid $\tau$ in darker green. The left and right pivots of
$\tau$ are $p_{\ell}$ and $p_{r}$, and the top edge is $e_{t}$. Right: The
same trapezoid $\tau$ as in the middle, which is created at the angle
$\varphi_{0}$ and removed at $\varphi_{1}$. The pair $p,q_{0}$ contributes to
$\tau$, but the pair $p,q_{1}$ does not.

The algorithm iterates over the edges of $P$ and computes for each edge $e$ a
contribution from pairs of points $p,q\in P$ that see each other and so that
the extension $pp^{\prime}$ of $qp$ from $p$ has the endpoint $p^{\prime}$ at
$e$; see Figure 1 (left). Here, the _extension_ of a segment $ab\subset P$
from $b$ is the maximal segment $bb^{\prime}$ on the same line as $ab$ such
that $bb^{\prime}\subset P$ and $ab\cap bb^{\prime}=\\{b\\}$. The extension
from $a$ is defined analogously. We assume without loss of generality that $e$
is horizontal with the interior of $P$ above $e$.

To compute the contribution corresponding to $e$, we perform a rotational
sweep around the endpoints of $e$ with two parallel rays. Informally, we keep
track of the region seen from $e$ in the strip between the rays and in the
direction of the rays. We say that a point $q\in P$ is _seen_ from $e$ at an
angle of $\varphi\in[0,\pi]$ if there exists $p\in e$ with $pq\subset P$ and
so that the counterclockwise angle from $e$ to $pq$ is $\varphi$. Define
$\textrm{vis}(e,\varphi)=\textrm{vis}_{P}(e,\varphi)$, $\varphi\in[0,\pi]$, to
be the region of points in $P$ seen from $e$ at an angle of $\varphi$; see
Figure 1 (middle). We keep track of the region $\textrm{vis}(e,\varphi)$ as
$\varphi$ increases from $0$ to $\pi$.

Let $d_{\varphi}=(\cos\varphi,\sin\varphi)$ be the direction corresponding to
the angle $\varphi$. For each corner $c$ of $P$ contained in
$\textrm{vis}(e,\varphi)$, let us draw a line through $c$ with direction
$d_{\varphi}$. These lines partition $\textrm{vis}(e,\varphi)$ into a set of
_rotating trapezoids_ (sometimes just called _trapezoids_ for brevity) which
we denote as $T(e,\varphi)$.

Each rotating trapezoid $\tau\in T(e,\varphi)$ is in fact a map from an
interval of angles $[\varphi_{0},\varphi_{1}]$ to trapezoids contained in $P$,
so that $\tau(\varphi)$ is a trapezoid in $\textrm{vis}(e,\varphi)$ bounded by
$e$, a pair of parallel edges with direction $d_{\varphi}$, and another edge
of $P$, called the _top_ edge of $\tau$. As the angle $\varphi$ increases, the
parallel edges are thus rotating in a parallel fashion, and each is rotating
around a corner of $P$. These corners are called the _left_ and _right pivot_
of $\tau$ and the associated edges of $\tau$ are the _left_ and _right_ edges
of $\tau$, respectively.

Whenever a corner of $P$ starts or stops being seen from $e$, trapezoids are
created in or removed from $T(e,\varphi)$, or both. Consider a trapezoid
$\tau$ that is created at some angle $\varphi_{0}$ and removed again at some
angle $\varphi_{1}>\varphi_{0}$; see Figure 1 (right). Once we get to the
angle $\varphi_{1}$ and observe that $\tau$ is removed, we compute a
contribution to $B(P)$ corresponding to $\tau$. That is, we integrate over all
pairs $p,q\in P$ where there exists $\varphi\in[\varphi_{0},\varphi_{1}]$ such
that $pq\in\tau(\varphi)$ and the oriented segment $pq$ has direction
$d_{\varphi}$.

#### Events.

The events of the rotational sweep are the angles $\varphi$ at which the
combinatorial structure of the trapezoids of $\textrm{vis}(e,\varphi)$ change.
This can happen because a corner becomes visible or invisible from $e$.
Whenever a change happens, we will compute a contribution to the Beer-index
from the trapezoids that stop existing at the event. It then follows that we
compute the contribution to the Beer-index $B(P)$ of all pairs of points in
$P$ that see each other, since every pair appears in a trapezoid that will
eventually stop existing. Furthermore, the pairs contributing to different
trapezoids are disjoint, except for the shared boundary between the
neighboring trapezoids, which is a set of measure $0$ in $P\times P$, so we do
not over-count any contribution.

Figure 2: An example with all the events that happen when scanning around the
edge $e$. The blue segments and points indicate the diagonal (and its
extensions) causing the event.

When a corner of $P$ becomes visible from $e$, an old trapezoid will be
replaced by two new ones, and we call it an _appearance_ event. Likewise, when
a corner of $P$ becomes invisible from $e$, two old trapezoids will be
replaced by one new, and we call it a _disappearance_ event. The angles
$\varphi=0$ and $\varphi=\pi$ are exceptional in that for $\varphi=0$, the
first trapezoid is created and for $\varphi=\pi$, the last trapezoid is
removed. See Figure 2 for an example of all events occurring in one scan.

Note that when a corner $c$ becomes visible, it is because an edge $f$ stops
blocking $c$ from being visible. Likewise, when a corner $c$ stops being
visible, it is because an edge $f$ starts blocking $c$ from being visible.
Therefore, appearance and disappearance events happen when $c$ and an endpoint
$c^{\prime}$ of $f$ are both visible at an angle of $\varphi$ from the same
point on $e$. We conclude that each event corresponds to a diagonal of $P$
whose extension intersects $e$. (Recall that the edges of $P$ are also
diagonals by definition.) Such a diagonal will either share an endpoint with
$e$ or the extension will have an endpoint at an interior point of $e$.

We first describe a version of the algorithm that handles the events in order
of angle. Here it seems unavoidable to use a priority queue to keep track of
the forthcoming events. As there are $O(n^{2})$ diagonals defining the events
and the queue operations take $O(\log n)$ time, we end up with a running time
of $O(n^{2}\log n)$. In the end of the section, we will explain how to get
down to $O(n^{2})$ time. This comes at the cost that we will generally not
handle the events in that order, so the faster algorithm may be conceptually
more challenging to understand, although it is strictly speaking a simpler
algorithm.

#### Computing diagonals.

In order to keep track of the events, we compute for each corner $c$ the set
$\mathcal{D}_{c}$ of all diagonals of $P$ with an endpoint at $c$ in order
sorted by angle. We can compute the diagonals by running a simple rotating
sweep-line algorithm centered at each corner [5], which takes $O(n^{2}\log n)$
time in total.

Figure 3: The full segments are diagonals from the corner $c$ and the dashed
segments are their extensions. The successors of $d_{0},d_{1},d_{2},d_{3}$ are
$d_{3},d_{2},d_{0},d_{1}$, respectively.

Let $c$ be a corner of $P$ and consider a segment $cd\subset P$ and the full
line $l(cd)$ containing $cd$. The _successor_ of $cd$ around $c$ is the
diagonal $cd_{0}$ we first meet when rotating $l(cd)$ counterclockwise around
$c$; see Figure 3. We will precompute the successor of each diagonal
$cd\in\mathcal{D}_{c}$ as follows. After the diagonals $\mathcal{D}_{c}$ have
been computed, we can perform a rotational sweep with a full line around $c$,
and every time we meet a diagonal, it is the successor of the previously
encountered diagonal. We can therefore precompute all successors of diagonals
around $c$ in $O(|\mathcal{D}_{c}|)$ time.

In order to find the top edges of trapezoids, we define for each oriented
diagonal $ab$ of $P$ an edge $\eta(ab)$ as follows; see Figure 4 for an
illustration. The edge $\eta(ab)$ has the informal description as the edge of
$P$ we will see when looking “between” $a$ and $b$ from an angle slightly
larger than that of $ab$. Let $c$ be the corner following $b$ on the boundary
$\partial P$, so that the interior of $P$ is to the left of $bc$. If $c$ is to
the left of the line containing $ab$, oriented from $a$ to $b$, or $c=a$, then
we define $\eta(ab)=bc$. Otherwise, let $bb^{\prime}$ be the extension of $ab$
from $b$. Then $\eta(ab)$ is the edge of $P$ containing $b^{\prime}$. By our
general position assumption, this edge is unique. We can easily define the
$\eta$-edges while finding the diagonals with a rotational sweep around each
corner.

Figure 4: The $\eta$-edge for three different diagonals in the same polygon
showing the different cases of the definition. Left: $c$ is to the left of
$ab$. Middle: $c=a$. Right: $c$ is to the right of $ab$.

#### Representing rotating trapezoids.

We maintain each existing trapezoid $\tau$ as a collection of a few pointers
to objects in $P$, as follows. We use a pointer to the left and right pivot
and the top edge. For each of the two pivots, we will furthermore keep track
of the successor of the current direction $d_{\varphi}$. The trapezoid $\tau$
will then be removed once we reach the event corresponding to the first of
these two successors. We store the trapezoids in a priority queue $Q$ ordered
by the angle of removal, so that the next trapezoid to leave $Q$ is the next
to be removed.

Figure 5: The three ways in which a trapezoid can be removed, corresponding to
cases 1–3.

#### Handling events.

Consider a trapezoid $\tau$ that we dequeue from $Q$. We will then compute the
contribution to $B(P)$ corresponding to $\tau$, which will be explained later.
The trapezoid $\tau$ is removed because the left edge $e_{\ell}$ hits another
corner than the left pivot $p_{\ell}$ or the right edge $e_{r}$ hits another
corner than the right pivot $p_{r}$. As the edges sweep in a counterclockwise
direction, this can only happen because of one of these three cases; see
Figure 5:

  1. 1.

The part of $e_{\ell}$ above $p_{\ell}$ hits a corner $c$ of $P$ that becomes
visible.

  2. 2.

The part of $e_{r}$ below $p_{r}$ hits a corner $c$ of $P$ and $p_{r}$ becomes
invisible.

  3. 3.

The part of $e_{\ell}$ below $p_{\ell}$ hits $p_{r}$ or the part of $e_{r}$
above $p_{r}$ hits $p_{\ell}$. Here, $p_{\ell}$ becomes invisible.

We handle the cases in the following way.

  1. 1.

We create two new trapezoids $\tau_{0}$ and $\tau_{1}$. Here, $\tau_{0}$ is
similar to $\tau$ but with $c$ as left pivot instead of $p_{\ell}$. The
trapezoid $\tau_{1}$ has $c$ as right pivot, $p_{\ell}$ as left pivot and
$\eta(p_{\ell}c)$ as top edge.

  2. 2.

We create a new trapezoid $\tau_{0}$ which is similar to $\tau$ but with $c$
as right pivot instead of $p_{r}$.

  3. 3.

We don’t create any new trapezoids.

When a new trapezoid $\tau^{\prime}$ is created, we know the diagonal
$d^{\prime}$ that caused $\tau^{\prime}$ to be created. Recall that we need to
find the successor of $d^{\prime}$ around each of the two pivots of
$\tau^{\prime}$. Let us consider an appearance event as in case 1. For the
trapezoid $\tau_{0}$, the successor around the right pivot $p_{r}$ can be
inherited from the removed trapezoid $\tau$, and the successor around the left
pivot $c$ is the precomputed successor of the diagonal $p_{\ell}c$. The
successors around the pivots of the trapezoid $\tau_{1}$ are the precomputed
successors of the diagonal $p_{\ell}c$. The new trapezoid created in case 2
can be handled in a similar way. We conclude that the successors can be found
in $O(1)$ time.

In order to prove that the algorithm works, we note that each trapezoid
$\tau$, except the very first one, has a unique predecessor $\tau^{\prime}$,
which is the one whose removal causes $\tau$ to be created. It then follows
that the algorithm will create and remove each trapezoid exactly once.

#### Contribution of one rotating trapezoid.

The perhaps most involved part of our algorithm is to compute the contribution
of a single trapezoid. Consider a trapezoid $\tau$ existing in the angle
interval $[\varphi_{0},\varphi_{1}]$. We evaluate the contribution of $\tau$
by integrating over quadruples of coordinates $(x_{0},y_{0},x_{1},y_{1})\in
P\times P$, such that the segment from $(x_{0},y_{0})$ to $(x_{1},y_{1})$ is
contained in $\tau(\varphi)$ and has direction $d_{\varphi}$, for some
$\varphi\in[\varphi_{0},\varphi_{1}]$.

We will apply a sequence of transformations to $\tau$ in order to obtain an
expression for the contribution of $\tau$ which is as simple as possible; see
Figure 6. In the previous description of the algorithm, we assumed without
further explanation that the edge $e$ is horizontal and that the interior of
$P$ is above $e$. If $e$ does not have this property, we first apply the
unique linear angle-preserving transformation $A$ that maps the endpoints $u$
and $v$ of $e$ to $(0,0)$ and $(1,0)$, respectively. The transformation $A$
scales all distances by a factor of $s_{0}=1/\|e\|_{2}$. Since we integrate
over a set of points in $P\times P\subset\mathbb{R}^{4}$, this causes the
integral to be scaled by $s_{0}^{4}$, so we should scale the final result by
$1/s_{0}^{4}$.

Figure 6: Left: Interpolation between the defining vectors $V_{0}$ and
$V_{1}$. Right: For all points $(h,0)$ on the segment $uv$ in the trapezoid,
we integrate over all pairs of points in the direction of $V_{\rho}$.

Denote by $x(\cdot)$ and $y(\cdot)$ the $x$\- and $y$-coordinates of a point
or a vector, respectively. Let the left and right pivots of $\tau$ be
$p_{\ell}$ and $p_{r}$, respectively. We use a horizontal translation to
obtain that $x(p_{\ell})=0$. Suppose that $x(p_{r})>0$ and define
$s_{1}=1/x(p_{r})$. Scaling the $x$-coordinates of all objects by $s_{1}$, we
can assume that $x(p_{r})=1$. Since we integrate over a set of points in
$P\times P\subset\mathbb{R}^{4}$ where the first and third coordinates are
$x$-coordinates, this causes the integral to be scaled by $s_{1}^{2}$, so we
should scale the final result by $1/s_{1}^{2}$. If instead $x(p_{r})<0$, we
can obtain $x(p_{r})=-1$. In the exceptional case that $x(p_{r})=0$, we
already have a situation with simple coordinates and can proceed in an
analogous way.

Let $f$ and $g$ be the endpoints of the top edge $e_{t}$ of $\tau$. We are
given $f$ and $g$ as particular points, but when computing the contribution of
$\tau$, we are free to choose other points on the same line containing the
edge $e_{t}$, i.e., we can “slide” $f$ and $g$ along the line containing
$e_{t}$ and choose some points with particularly simple coordinates on the
line. If the edge $fg$ is not vertical, we can therefore assume that $x(f)=1$
and $x(g)=0$. In the exceptional case that $fg$ is vertical, we can instead
assume that $y(f)=0$ and $y(g)=1$. In the following, we will focus on the
general case where we have $x(f)=1$ and $x(g)=0$, and in the other case, we
could proceed in an analogous way. Let $s_{2}=1/y(f)$ and note that by scaling
the $y$-coordinates of all objects by $s_{2}$, we obtain that $y(f)=1$. This
again causes the result of the integral to be scaled by $s_{2}^{2}$, so we
should multiply the final result by $1/s_{2}^{2}$. We then have $f=(1,1)$.

The start and end angles $\varphi_{0},\varphi_{1}$ are defined by vectors
$V_{0}$ and $V_{1}$, respectively. Note that $y(V_{0})\geq 0$ and
$y(V_{1})\geq 0$, as each vector points from a corner above $e$ to another
corner further above $e$. In particular, the counterclockwise angle from
$V_{0}$ to $V_{1}$ is also at most $\pi$. If the angle is at least $\pi/4$, we
introduce one or two auxiliary vectors that lie in between the angular span
from $V_{0}$ and $V_{1}$ so that we “split” the trapezoid $\tau$ into two or
three trapezoids. These have the same pivots and top edge as $\tau$, but the
range of rotation of each of them is only a subset of the range of the
original trapezoid $\tau$. In the following, we therefore assume that the
angle from $V_{0}$ to $V_{1}$ is less than $\pi/2$.

In the general case that $y(V_{0})>0$ and $y(V_{1})>0$, we can scale $V_{0}$
and $V_{1}$ to obtain $y(V_{0})=y(V_{1})=1$. Otherwise, note that it is
impossible that $y(V_{0})=0$ and $y(V_{1})=0$, since this would mean that a
third corner of $P$ is on the line through $uv$, which is excluded by our
general position assumption. If $y(V_{0})=0$ and $x(V_{0})>0$, we know that
$x(V_{1})>0$ and $y(V_{1})>0$ since the angle from $V_{0}$ to $V_{1}$ is less
than $\pi/2$, so we can scale to obtain $V_{0}=(1,0)$ and $x(V_{1})=1$.
Otherwise, we have that $y(V_{1})=0$, $x(V_{0})<0$ and $x(V_{1})<0$, and we
scale to obtain $V_{0}=(-1,0)$ and $x(V_{1})=-1$. In the following, we will
focus on the general case where we have $y(V_{0})=y(V_{1})=1$. Expressions for
the contributions in the other cases can be computed in a similar way.

For $\rho\in[0,1]$, we define a corresponding interpolated vector
$V_{\rho}=V_{0}\cdot(1-\rho)+V_{1}\cdot\rho$, which corresponds to the time
where the left and right edges of $\tau$ are parallel to $V_{\rho}$. In order
to evaluate the contribution of $\tau$, we will use a transformation of
coordinates. Our new coordinates are quadruples $(\rho,h,p,q)$, where $h$
denotes where we “look” from the $x$-axis in direction $V_{\rho}$, and $p$ and
$q$ are scalars that can be multiplied with $V_{\rho}$ to obtain the two
points in $\tau$ we are integrating over. The usual coordinates cooresponding
to the quadruple $(\rho,h,p,q)$ are given by the transformation

$\displaystyle F(\rho,h,p,q)$ $\displaystyle=(h+p\cdot x(V_{\rho}),p\cdot
y(V_{\rho}),h+q\cdot x(V_{\rho}),q\cdot y(V_{\rho}))$ $\displaystyle=(h+p\cdot
x(V_{\rho}),p,h+q\cdot x(V_{\rho}),q).$

We thus have the Jacobian

$J_{F}=\begin{bmatrix}p\cdot x(V_{1}-V_{0})&1&x(V_{\rho})&0\\\ 0&0&1&0\\\
q\cdot x(V_{1}-V_{0})&1&0&x(V_{\rho})\\\ 0&0&0&1\end{bmatrix},$

with the determinant

$\det J_{F}=(q-p)\cdot x(V_{0})+(p-q)\cdot x(V_{1}).$

We now describe the bounds of the integral that we evaluate in order to
compute the contribution to the Beer-index of the trapezoid $\tau$. Let
$h_{\ell}(\rho)$ be the $x$-coordinate of the intersection of $e$ with the
line through the left pivot $p_{\ell}$ with direction $V_{\rho}$ and define
$h_{r}(\rho)$ analogously for the right pivot $p_{r}$. For a value
$h\in[h_{\ell}(\rho),h_{r}(\rho)]$, let $T(\rho,h)\geq 0$ be the
$y$-coordinate of the intersection of the line through $(h,0)$ with direction
$V_{\rho}$ and the line containing $e_{t}$. The contribution of $\tau$ can
then be expressed as

$\displaystyle\int_{0}^{1}\int_{h_{\ell}(\rho)}^{h_{r}(\rho)}\int_{0}^{T(\rho,h)}\int_{0}^{q}\left|\det
J_{F}\right|\mathop{}\\!\mathrm{d}p\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}h\mathop{}\\!\mathrm{d}\rho.$
(1)

We now show that there is a closed-form expression equal to this integral.
Note that as $q\geq p$ and $x(V_{0})>x(V_{1})$, we have $\det J_{F}>0$, so the
absolute value can be ignored. It is straightforward to check that

$\int_{0}^{q}\det
J_{F}\mathop{}\\!\mathrm{d}p=\frac{q^{2}(x(V_{0})-x(V_{1}))}{2}.$

One can verify that

$T(\rho,h)=\frac{(h-1)\,y(g)-h}{(x(V_{0})\,(\rho-1)-x(V_{1})\,\rho)\,y(g)+x(V_{0})\,(1-\rho)+x(V_{1})\,\rho-1}.$

Since we have the antiderivative

$\int\frac{q^{2}(x(V_{0})-x(V_{1}))}{2}\mathop{}\\!\mathrm{d}q=\frac{q^{3}\,(x(V_{0})-x(V_{1}))}{6},$

we get

$\displaystyle\int_{0}^{T(\rho,h)}\int_{0}^{q}\left|\det
J_{F}\right|\mathop{}\\!\mathrm{d}p\mathop{}\\!\mathrm{d}q$
$\displaystyle=\left[\frac{q^{3}\,(x(V_{0})-x(V_{1}))}{6}\right]_{0}^{T(\rho,h)}$
$\displaystyle=\frac{((h-1)\,y(g)-h)^{3}\,(x(V_{0})-x(V_{1}))}{6\,((x(V_{0})\,(\rho-1)-x(V_{1})\,\rho)\,y(g)+x(V_{0})\,(1-\rho)+x(V_{1})\,\rho-1)^{3}}.$
(2)

We note that (2) is a cubic polynomial in $h$, so it is straightforward to
compute the antiderivative

$\int\int_{0}^{T(\rho,h)}\int_{0}^{q}\left|\det
J_{F}\right|\mathop{}\\!\mathrm{d}p\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}h,$

which will be a quartic polynomial in $h$. Since we have

$\displaystyle h_{\ell}(\rho)$
$\displaystyle=y(p_{\ell})\,(x(V_{0})\,(\rho-1)-x(V_{1})\,\rho)\quad\quad\quad\textrm{and}$
$\displaystyle h_{r}(\rho)$
$\displaystyle=y(p_{r})\,(x(V_{0})\,(\rho-1)-x(V_{1})\,\rho)+1,$

we then get that the integral

$\displaystyle\int_{h_{\ell}(\rho)}^{h_{r}(\rho)}\int_{0}^{T(\rho,h)}\int_{0}^{q}\left|\det
J_{F}\right|\mathop{}\\!\mathrm{d}p\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}h$
(3)

is a rational function in $\rho$ (note that $\rho$ appears in the denominator
in (2), so we don’t get a polynomial in $\rho$).

Now, as (3) is a rational function, it is well-known [7, ch. 2] that the
antiderivative

$\int\int_{h_{\ell}(\rho)}^{h_{r}(\rho)}\int_{0}^{T(\rho,h)}\int_{0}^{q}\left|\det
J_{F}\right|\mathop{}\\!\mathrm{d}p\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}h\mathop{}\\!\mathrm{d}\rho$

can be written as a function of $\rho$ involving only the usual arithmetic
operations, logarithms, and arctangents. This finishes the proof that a
closed-form expression exists.

Using Maple [24], we have computed an explicit formula for (1) involving the
five parameters $y(p_{\ell})$, $y(p_{r})$, $y(g)$, $x(V_{0})$, $x(V_{1})$. The
formula does involve logarithms, but not arctangents. Unfortunately, the
formula is too long to be written here, but it can be seen in the output from
Maple in Appendix C.

#### Running time.

We use $O(n^{2}\log n)$ time to find all diagonals in sorted order around each
corner of $P$.

It is not hard to construct a polygon $P$ (with holes) where we create and
remove $\Omega(n^{2})$ trapezoids when sweeping around a single edge $e$.
However, we claim that when scanning over all edges $e$, the total number of
handled trapezoids will also be $O(n^{2})$.

Figure 7: Three cases showing when the diagonal $ab$ causes an event to
happen. The gray regions indicate the exterior of $P$.

To this end we observe that there are at most two edges such that when
sweeping around these, a diagonal $ab$ can cause an event to happen. At most
two trapezoids can be created in each event, so the claim will follow. Recall
that the diagonal $ab$ will only correspond to an event when sweeping around
an edge $e$ if the extension $aa^{\prime}$ from $a$ intersects $e$. Let $c$
and $d$ be the neighboring corners of $a$ on the boundary $\partial P$ such
that the interior of $P$ is to the left of $ca$ and $ad$; see Figure 7. If $c$
and $d$ are on different sides of $ab$, then it must be the case that $c$ is
to the left and $d$ is to the right, since otherwise, $ab$ would not be
contained in $P$. We conclude that $a=a^{\prime}$ and that $ab$ only causes an
event when we sweep around $ca$ and $ad$.

Suppose on the other hand that $c$ and $d$ are on the same side of $ab$. Then
$ab$ causes an event when sweeping around the edge containing $a^{\prime}$,
which is unique by the general position assumption. If $c$ and $d$ are to the
left of $ab$, then $ab$ also causes an event when sweeping around $ca$, and
otherwise when sweeping around $ad$.

We use $O(\log n)$ time to handle each trapezoid, which is the time it takes
to dequeue from $Q$ and insert one or two trapezoids in $Q$. Since each of the
$O(n^{2})$ diagonals can lead to at most two trapezoids and it takes $O(\log
n)$ time to handle each trapezoid, we use $O(n^{2}\log n)$ time in total.

#### Getting to $O(n^{2})$ time.

Using the algorithm from [2], we can compute the visibility polygon
$\textrm{vis}(c)$ of each corner $c$ in total time $O(n^{2})$. We then know
all diagonals and can likewise compute the successors and $\eta$-edges in
$O(n^{2})$ total time.

We need to avoid using a priority queue, as it takes $\Omega(\log n)$ time to
insert or dequeue. We can instead use an ordinary queue (or even a stack). We
will then not handle the trapezoids in order of when they are removed, but we
will still handle each trapezoid exactly once. We then obtain Theorem 1.

##  3 Counting visible pairs

In this section, we describe an algorithm for computing the number of visible
pairs among a set $M$ of $m$ points contained in a simple polygon $P$ with $n$
corners. It is a divide-and-conquer algorithm, where we choose a diagonal $uv$
of $P$ that separates $P$ into two sub-polygons $P_{1}$ and $P_{2}$ with
nearly equally many corners. We then count the number $k$ of pairs of one
point in $P_{1}$ and the other in $P_{2}$ that see each other across $uv$ and
then recurse on $P_{1}$ and $P_{2}$. We are able to compute $k$ in time
$O(n+m^{4/3}\log^{\alpha}m)$, where $\alpha<1.78$ is a constant. Computing the
number of pairs of points in $M$ that see each other in a single triangle
trivially takes $O(1)$ time once we know how many points from $M$ the triangle
contains. Let $\triangle=n-2$ be the number of triangles in a triangulation of
$P$. The recursive formula for the running time $T(\triangle,m)$ is then

$\displaystyle T(1,m)$ $\displaystyle=O(1),$ $\displaystyle T(\triangle,m)$
$\displaystyle=T(\triangle_{1},m_{1})+T(\triangle_{2},m_{2})+O(\triangle+m^{4/3}\log^{\alpha}m),\quad\triangle>1,$

where $\triangle_{1},\triangle_{2}\leq 2\triangle/3$,
$\triangle_{1}+\triangle_{2}=\triangle$ and $m_{1}+m_{2}=m$. It follows that
the algorithm has running time
$T(\triangle,m)=O(\triangle\log\triangle+m^{4/3}\log^{\alpha}m\log\triangle)=O(n\log
n+m^{4/3}\log^{\alpha}m\log n)$.

In order to find a suitable diagonal, we first triangulate $P$, for instance
using the simple algorithm from [5] running in $O(n\log n)$ time. It is well
known that there exists a diagonal $uv$ in the triangulation that splits $P$
into two sub-polygons $P_{1}$ and $P_{2}$, each consisting of at most $\lceil
2\triangle/3\rceil$ triangles. By partitioning $P_{1}$ and $P_{2}$
recursively, we obtain a balanced hierarchical decomposition of $P$
represented as a balanced binary tree $\mathcal{T}$. Here, the root of
$\mathcal{T}$ represents the diagonal $uv$, the left and right subtrees
represent $P_{1}$ and $P_{2}$, respectively, and the leafs represent the
individual triangles.

The triangulation of $P$ induces a dual graph $G$ which is a tree with maximum
degree $3$, and the diagonal $uv$ can be chosen as incident to the triangle
corresponding to the centroid of $G$. We can therefore find $uv$ in $O(n)$
time by a simple traversal of $G$, which yields an algorithm for constructing
the tree $\mathcal{T}$ in $O(n\log n)$ time. By a more refined approach,
$\mathcal{T}$ can also be computed in $O(n)$ time [17].

Figure 8: Left: We find the visible parts of the diagonal $uv$ from the three
points $p,q,r$. Right: We map the visibility segments to the dual space and
identify the segment $d=pr$ as the intersection $d^{*}$.

Our algorithm relies on the principle of geometric duality [5]. Here, a point
$p=(a,b)$ is mapped to a line $p^{*}:y=ax-b$ and a line $\ell:y=cx+d$ is
mapped to a point $\ell^{*}=(c,-d)$. Assume without loss of generality that
the splitting diagonal $uv$ is on the $y$-axis. For a segment $s$, denote by
$l(s)$ the full line containing $s$. For each point $p\in M$ and every point
$q\in uv$ that $p$ sees, we draw the point $l(pq)^{*}$ in the dual; see Figure
8. Since $p$ sees an interval (or nothing) of $uv$ and all the lines
$l(pq)^{*}$ pass through the same point $p$, these dual points form a segment
$s_{p}$.

Our algorithm relies on the following observation: If $p_{1},p_{2}\in M$ are
on different sides of $uv$ and the segments $s_{p_{1}}$ and $s_{p_{2}}$
intersect, then $p_{1}$ and $p_{2}$ see each other across $uv$. To see this,
note that the intersection point $r^{*}$ of $s_{p_{1}}$ and $s_{p_{2}}$
corresponds to a line $r$ in primal space that contains visibility segments
from both $p_{1}$ and $p_{2}$ to $uv$. Therefore, $p_{1}$ and $p_{2}$ see each
other across $uv$.

When $p\in M$ is in $P_{1}$, we color the segment $s_{p}$ blue. Otherwise,
when $p$ is in $P_{2}$, we color $s_{p}$ red. Then the number of pairs that
see each other across $uv$ is exactly the number of intersections between a
red and a blue segment in the dual space. Agarwal [1] gave an algorithm for
counting such bichromatic intersections in $O(m^{4/3}\log^{\alpha}m)$ time,
where $\alpha<1.78$ is a constant. The algorithm uses $O(m^{4/3})$ space.

In order to construct the segments $s_{p}$ efficiently, we use a data
structure described by Guibas, Hershberger, Leven, Sharir and Tarjan [17].
Given the segment $uv$ in $P$, we can preprocess $P$ in $O(n)$ time so that
given a query point $p$, we can compute the part of $uv$ that is visible from
$p$ in $O(\log n)$ time. The authors assumed that $P$ was given together with
a triangulation, presumably because the paper appeared before it was known
that a triangulation can be found in $O(n)$ time [12].

In order to handle the recursive calls, we need to find the sets $M_{1}=M\cap
P_{1}$ and $M_{2}=M\cap P_{2}$. This can be done by creating a data structure
supporting ray shooting queries in $P$ in $O(\log n)$ time, as described by
Hershberger and Suri [19]. For each point $m\in M$, we shoot a ray vertically
up and check if the point where we hit the boundary is in $P_{1}$ or $P_{2}$.
The data structure takes $O(n)$ time to construct. We therefore use $O(n+m\log
n+m^{4/3}\log^{\alpha}m)=O(n+m^{4/3}\log^{\alpha}m)$ time on preprocessing,
constructing the segments $s_{p}$ for all $p\in M$, counting the intersections
and finding the sets $M_{1}$ and $M_{2}$. Taking the recursive calls into
account, this results in an algorithm with total running time $O(n\log
n+m^{4/3}\log^{\alpha}m\log n)$, proving Theorem 2.

##  4 Expected distances in the $L_{1}$-metric

###  4.1 Finite set of points in $\mathbb{R}^{d}$

Consider $m$ points
$M=\left\\{p_{1},p_{2},\ldots,p_{m}\right\\}\subseteq\mathbb{R}^{d}$, where we
want to compute the sum of the $L_{1}$-distance between all pairs. Our
proposed method demonstrates a technique that we will also use in Section 4.2
when computing the expected $L_{1}$-distance between a pair of random points
in a polygon, so this section serves as a warmup before that algorithm is
presented. The overall idea is to compute the contribution to the total
distance from each dimension separately.

Figure 9: The width of the slab $S$ is part of the difference between the
$x$-coordinates of each pair of one point among $p_{1},\ldots,p_{\pi(i)}$ and
another among $p_{\pi(i+1)},\ldots,p_{n}$, and the number of such pairs is
$i\cdot(m-i)$.

We first explain how the algorithm works in the planar case $d=2$ and then
observe that it generalizes to an arbitrary dimension $d$. When $d=2$, the
algorithm works by computing the total differences of $x$\- and
$y$-coordinates of the points separately and then adding these numbers. Let
$\pi$ be a permutation of $\left\langle 1,2,\ldots,m\right\rangle$ that sorts
the points by $x$-coordinates, i.e., $x(p_{\pi(i)})\leq x(p_{\pi(i+1)})$ for
all $1\leq i<m$. Let us consider the slab $S$ between the vertical lines
through $x(p_{\pi(i)})$ and $x(p_{\pi(i+1)})$ and count how many pairs of
points are separated by $S$; see Figure 9. As there are $i$ points to the left
of $S$ and $m-i$ points to the right, there are $i\cdot(m-i)$ such pairs, and
the width of $S$ is part of the difference of $x$-coordinates for all these
pairs. Thus, the contribution to the total difference of $x$-coordinates of
$S$ is $i\cdot(m-i)\cdot(x(p_{\pi(i+1)})-x(p_{\pi(i)}))$. Summing over all
values of $i$, we get the total difference of $x$-coordinates over all pairs
of points, i.e.,

$\displaystyle\sum_{i=1}^{m}\sum_{j=1}^{i-1}(x(p_{\pi(i)})-x(p_{\pi(j)}))=\sum_{i=1}^{m-1}i\cdot(m-i)\cdot(x(p_{\pi(i+1)})-x(p_{\pi(i)})).$
(4)

Note that the sum (4) can be evaluated in $O(m)$ time once the permutation
$\pi$ has been found. We can evaluate the total difference of the
$y$-coordinates in an analogous way. Hence, we can compute the total
$L_{1}$-distance in $O(m\log m)$ time.

For points of arbitrary dimension $d$, consider a coordinate
$k\in\\{1,\ldots,d\\}$, and let $\pi_{k}$ be a permutation of $\left\langle
1,2,\ldots,m\right\rangle$ sorting the points with respect to this coordinate,
i.e., we have $p_{\pi_{k}(i),k}\leq p_{\pi_{k}(i+1),k}$ for all $1\leq i<m$,
where $p_{i,k}$ denotes the $k$th coordinate of $p_{i}$. Using a formula
analogous to (4) for each coordinate separately, we get

$\displaystyle\sum_{i=1}^{m}\sum_{j=1}^{i-1}||p_{i}-p_{j}||_{1}$
$\displaystyle=\sum_{i=1}^{m}\sum_{j=1}^{i-1}\sum_{k=1}^{d}|p_{i,k}-p_{j,k}|$
$\displaystyle=\sum_{k=1}^{d}\sum_{i=1}^{m}\sum_{j=1}^{i-1}|p_{i,k}-p_{j,k}|$
$\displaystyle=\sum_{k=1}^{d}\sum_{i=1}^{m-1}i\cdot(m-i)\cdot(p_{\pi_{k}(i+1),k}-p_{\pi_{k}(i),k}).$

The above implies that the sum of all pairwise distances may be computed in
time $O(d\cdot m\log m)$ by sorting the coordinates according to their $k$th
coordinate and evaluating the inner sum for each $k$, which establishes
Theorem 4.

###  4.2 Expected $L_{1}$-distance in a polygon

We are given a simple polygon $P$ with $n$ corners. Let us assume for
simplicity that no two corners of $P$ have the same $x$-coordinates or the
same $y$-coordinates. Note that a shortest path in $P$ in the $L_{1}$-metric
can be realized as a union of horizontal and vertical segments. For points
$a,b\in P$, denote by $a\,\includegraphics[page=16,scale={0.7}]{figs.pdf}\,b$
such a shortest path from $a$ to $b$. More generally, when $a$ or $b$ is a
compact set in $P$, we also use
$a\,\includegraphics[page=16,scale={0.7}]{figs.pdf}\,b$ to denote a shortest
rectilinear path $P$ from a point in $a$ to a point in $b$. Let
$|a\,\includegraphics[page=16,scale={0.7}]{figs.pdf}\,b|$ be the
$L_{1}$-length of the path
$a\,\includegraphics[page=16,scale={0.7}]{figs.pdf}\,b$. As in our warmup-case
of points in $\mathbb{R}^{d}$ from Section 4.1, we find the expected distance
by computing the contribution from horizontal and vertical segments
separately.

Let $\delta_{1}(p)$ (resp. $\delta_{2}(p)$) be the segments of a path $p$ that
are horizontal (resp. vertical) and define
$\Delta_{i}(a,b)=\sum_{s\in\delta_{i}(a\,\includegraphics[page=16,scale={0.7}]{figs.pdf}\,b)}|s|$.
Define for $A,B\subset P$ the notation

$D_{i}(A,B)=\iint_{p\in A}\iint_{q\in
B}\Delta_{i}(p,q)\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}p.$

We now observe

$\displaystyle\iint_{p\in P}\iint_{q\in
P}\left|p\,\includegraphics[page=16,scale={0.7}]{figs.pdf}\,q\right|\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}p$
$\displaystyle=\iint_{p\in P}\iint_{q\in
P}\left(\Delta_{1}(p,q)+\Delta_{2}(p,q)\right)\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}p$
$\displaystyle=D_{1}(P,P)+D_{2}(P,P).$

We explain how to compute $D_{1}(P,P)$ in $O(n)$ time, and $D_{2}(P,P)$ can be
computed in an analogous way, leading to an algorithm for computing the
expected $L_{1}$-distance with linear running time. We first construct a
vertical trapezoidation of $P$: For each corner $c$ of $P$, we add the maximal
vertical segment contained in $P$ and containing $c$; see Figure 10. These
vertical segments partition $P$ into trapezoids, each of which has a pair of
vertical edges. Some trapezoids degenerate into triangles. The trapezoids
induce a tree $T$, where the vertices are the trapezoids and two vertices are
connected by an edge if the trapezoids share a vertical edge. We choose an
arbitrary root $r$ of $T$, which induce parent-child relationships among
neighbouring pairs of trapezoids in $T$. Define $p(t)$ to be the parent of $t$
and $c(t)$ to be the set of children of $t$.

For each trapezoid $t$ which is not the root, let $l(t)$ be vertical edge of
$t$ that separates $t$ from the parent $p(t)$. Furthermore, let $P[t]\subset
P$ be the region consisting of $t$ and all decendants of $t$. Define

$L_{i}(t)=\iint_{p\in P[t]}\Delta_{i}(p,l(t))\mathop{}\\!\mathrm{d}p.$

We now show how we can compute the numbers $D_{1}(P[t],P[t])$, $L_{1}(t)$ and
$|P[t]|$ for all trapezoids $t$ in $O(n)$ time in total. We have then in
particular computed $D_{1}(P,P)=D_{1}(P[r],P[r])$.

It is trivial to evaluate each of the numbers $D_{1}(P[t],P[t])$, $L_{1}(t)$
and $|P[t]|=|t|$ for a leaf $t$ of $T$ in $O(1)$ time, since $t$ has
complexity $O(1)$, so suppose now that $t$ is not a leaf.

By our general position assumption, a trapezoid $t$ has at most four children.
Assuming that the values $D_{1}(P[t^{\prime}],P[t^{\prime}])$,
$L_{1}(t^{\prime})$ and $|P[t^{\prime}]|$ have been computed for each of the
children $t^{\prime}$ of $t$, we show how $D_{1}(P[t],P[t])$, $L_{1}(t)$ and
$|P[t]|$ can be computed in $O(1)$ time. The claim that we can compute the
numbers for all $t$ in $O(n)$ time then follows.

The area $|P[t]|$ can be simply computed as

$|P[t]|=|t|+\sum_{t^{\prime}\in c(t)}|P[t^{\prime}]|.$

We now explain how to exaluate $L_{1}(t)$, so suppose that $t\neq r$. We have

$L_{1}(t)=\iint_{p\in
t}\Delta_{1}(p,l(t))\mathop{}\\!\mathrm{d}p+\sum_{t^{\prime}\in
c(t)}\iint_{p\in P[t^{\prime}]}\Delta_{1}(p,l(t))\mathop{}\\!\mathrm{d}p.$

Clearly, the first term can be evaluated in $O(1)$ time, as the trapezoid $t$
has complexity $O(1)$. There are two cases when evaluating the integral in the
sum in the second term: (i) $t^{\prime}$ and $p(t)$ are on the same side of
$t$, and (ii) $t^{\prime}$ and $p(t)$ are on different sides of $t$; see
Figure 10. In case (i), the shortest path from a point $p\in P[t^{\prime}]$ to
$l(t)$ first follows a shortest path to $l(t^{\prime})$ and then follows a
vertical segment to $l(t)$, so we have

$\iint_{p\in
P[t^{\prime}]}\Delta_{1}(p,l(t))\mathop{}\\!\mathrm{d}p=\iint_{p\in
P[t^{\prime}]}\Delta_{1}(p,l(t^{\prime}))\mathop{}\\!\mathrm{d}p=L_{1}(t^{\prime}).$

Figure 10: A vertical trapezoidation of a simple polygon $P$ and the tree it
induces. Case (i) and (ii) in the computation of $L_{1}(t)$ are shown as
shortest paths from points $p^{\prime}\in P[t^{\prime}]$ and
$p^{\prime\prime}\in P[t^{\prime\prime}]$ for children $t^{\prime}$ and
$t^{\prime\prime}$ of $t$, respectively, to the segment $l(t)$.

In case (ii), the shortest path first follows a shortest path from $p$ to
$l(t^{\prime})$ and then crosses $t$ to get to $l(t)$. We then have

$\iint_{p\in
P[t^{\prime}]}\Delta_{1}(p,l(t))\mathop{}\\!\mathrm{d}p=\iint_{p\in
P[t^{\prime}]}(\Delta_{1}(p,l(t^{\prime}))+w(t))\mathop{}\\!\mathrm{d}p=L_{1}(t^{\prime})+\left|P[t^{\prime}]\right|w(t),$

where $w(t)$ is the width of $t$.

In order to evaluate $D_{1}(P[t],P[t])$, we observe

$D_{1}(P[t],P[t])=D_{1}(t,t)+2\sum_{t^{\prime}\in
c(t)}D_{1}(P[t^{\prime}],t)+\sum_{t^{\prime},t^{\prime\prime}\in
c(t)}D_{1}(P[t^{\prime}],P[t^{\prime\prime}]).$

Again, the first term takes $O(1)$ to evaluate as $t$ has complexity $O(1)$.
To evaluate the sum in the second term, we observe

$\displaystyle D_{1}(P[t^{\prime}],t)$ $\displaystyle=\iint_{p^{\prime}\in
P[t^{\prime}]}\iint_{p\in
t}\Delta_{1}(p^{\prime},p)\mathop{}\\!\mathrm{d}p\mathop{}\\!\mathrm{d}p^{\prime}$
$\displaystyle=\iint_{p^{\prime}\in P[t^{\prime}]}\iint_{p\in
t}(\Delta_{1}(p^{\prime},l(t^{\prime}))+\Delta_{1}(l(t^{\prime}),p))\mathop{}\\!\mathrm{d}p\mathop{}\\!\mathrm{d}p^{\prime}$
$\displaystyle=\left|t\right|L_{1}(t^{\prime})+\left|P[t^{\prime}]\right|\iint_{p\in
t}\Delta_{1}(l(t^{\prime}),p)\mathop{}\\!\mathrm{d}p.$

Again, the integral in the final expression can be evaluated in $O(1)$ time,
so the same holds for $D_{1}(P[t^{\prime}],t)$.

Consider now a term of the form $D_{1}(P[t^{\prime}],P[t^{\prime\prime}])$. If
$t^{\prime}=t^{\prime\prime}$, the number
$D_{1}(P[t^{\prime}],P[t^{\prime\prime}])$ has already been computed, so
suppose $t^{\prime}\neq t^{\prime\prime}$. We have two cases: (i) $t^{\prime}$
and $t^{\prime\prime}$ are on the same side of $t$ and (ii) $t^{\prime}$ and
$t^{\prime\prime}$ are on different sides. In case (i), we get

$\displaystyle D_{1}(P[t^{\prime}],P[t^{\prime\prime}])$
$\displaystyle=\iint_{p^{\prime}\in P[t^{\prime}]}\iint_{p^{\prime\prime}\in
P[t^{\prime\prime}]}\Delta_{1}(p^{\prime},p^{\prime\prime})\mathop{}\\!\mathrm{d}p^{\prime\prime}\mathop{}\\!\mathrm{d}p^{\prime}$
$\displaystyle=\iint_{p^{\prime}\in P[t^{\prime}]}\iint_{p^{\prime\prime}\in
P[t^{\prime\prime}]}(\Delta_{1}(p^{\prime},l(t^{\prime}))+\Delta_{1}(l(t^{\prime\prime}),p^{\prime\prime}))\mathop{}\\!\mathrm{d}p^{\prime\prime}\mathop{}\\!\mathrm{d}p^{\prime}$
$\displaystyle=\left|P[t^{\prime}]\right|L_{1}(t^{\prime\prime})+\left|P[t^{\prime\prime}]\right|L_{1}(t^{\prime}).$

In case (ii), we get

$\displaystyle D_{1}(P[t^{\prime}],P[t^{\prime\prime}])$
$\displaystyle=\iint_{p^{\prime}\in P[t^{\prime}]}\iint_{p^{\prime\prime}\in
P[t^{\prime\prime}]}\Delta_{1}(p^{\prime},p^{\prime\prime})\mathop{}\\!\mathrm{d}p^{\prime\prime}\mathop{}\\!\mathrm{d}p^{\prime}$
$\displaystyle=\iint_{p^{\prime}\in P[t^{\prime}]}\iint_{p^{\prime\prime}\in
P[t^{\prime\prime}]}(\Delta_{1}(p^{\prime},l(t^{\prime}))+w(t)+\Delta_{1}(l(t^{\prime\prime}),p^{\prime\prime}))\mathop{}\\!\mathrm{d}p^{\prime\prime}\mathop{}\\!\mathrm{d}p^{\prime}$
$\displaystyle=\left|P[t^{\prime}]\right|L_{1}(t^{\prime\prime})+\left|P[t^{\prime\prime}]\right|L_{1}(t^{\prime})+\left|P[t^{\prime}]\right|\left|P[t^{\prime\prime}]\right|w(t).$

Since a trapezoidation of $P$ can be computed in $O(n)$ time using Chazelle’s
algorithm [12], we get Theorem 5. Instead of Chazelle’s algorithm, which is
known for being complicated, we can also use a simpler sweep-line algorithm
[5] and obtain an algorithm running in time $O(n\log n)$.

##  5 Expected distances in the $L_{2}$-metric

We now show how to compute the expected geodesic $L_{2}$-distance between two
random points in a simple polygon $P$. For $p,q\in P$, we denote by $p\leadsto
q$ the unique shortest path in $P$ from $p$ to $q$ in the $L_{2}$-metric, and
we denote by $\|p\leadsto q\|$ the length of this path. We thus want to
compute

$\mathbb{E}(\|p\leadsto q\|)=\frac{1}{|P|^{2}}\iint_{p\in P}\iint_{q\in
P}\left\|p\leadsto q\right\|\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}p.$

We split this integral into two parts; see Figure 11. Let $A\subset P\times P$
be the pairs $(p,q)$ such that $p\leadsto q=pq$, i.e., $p$ and $q$ see each
other in $P$. Let $B\subset P\times P$ be the pairs $(p,q)$ such that
$p\leadsto q$ contains one or more corners of $P$. Then $A\cup B=P\times P$.
Furthermore, the intersection $A\cap B$ contains the pairs $(p,q)$ such that
$pq\subset P$ and $pq$ contains a corner of $P$. This set has measure $0$ in
$P\times P$, so we have

$\iint_{p\in P}\iint_{q\in P}\left\|p\leadsto
q\right\|\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}p=\iiiint_{(p,q)\in
A}\left\|pq\right\|\mathop{}\\!\mathrm{d}(p,q)+\iiiint_{(p,q)\in
B}\left\|p\leadsto q\right\|\mathop{}\\!\mathrm{d}(p,q).$

In the following subsections, we describe how we compute the integral over $A$
and $B$, respectively. In each subsection, we explain how to output $O(n^{2})$
constants defined as closed-form expressions of size $O(1)$. As described in
Theorem 6, some of these will refer to previously defined constants, and in
the end, the integral over both $A$ and $B$ can be expressed as a sum of a
subset of these $O(n^{2})$ constants, and the theorem follows.

Figure 11: A polygon $P$ with $(p_{1},q_{1})\in A$ and $(p_{2},q_{2})\in B$.

###  5.1 Shortest paths that are segments

Let us first consider the computation of the integral over $A$. Note that when
computing the Beer-index as described in Section 2, we were also integrating
over the set $A$ of pairs of points that see each other. We will therefore use
a similar algorithm here that evaluates the contribution from each rotating
trapezoid. The only difference is that instead of computing the volume in
$\mathbb{R}^{4}$ corresponding to each trapezoid, we will now integrate the
distance function. In other words, instead of the integral (1), we end up with

$\displaystyle\int_{0}^{1}\int_{h_{\ell}(\rho)}^{h_{r}(\rho)}\int_{0}^{T(\rho,h)}\int_{0}^{q}\left\|q\cdot
V_{\rho}-p\cdot V_{\rho}\right\|\left|\det
J_{F}\right|\mathop{}\\!\mathrm{d}p\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}h\mathop{}\\!\mathrm{d}\rho=$
$\displaystyle\int_{0}^{1}\int_{h_{\ell}(\rho)}^{h_{r}(\rho)}\int_{0}^{T(\rho,h)}\int_{0}^{q}\left\|V_{\rho}\right\|\left(q-p\right)\left|\det
J_{F}\right|\mathop{}\\!\mathrm{d}p\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}h\mathop{}\\!\mathrm{d}\rho.$
(5)

We would again like to transform the trapezoid to make it simpler, but here we
need to be more careful than we were in the case of the Beer-index. In
particular, we cannot scale the $x$\- and $y$-coordinates independently, since
such a scaling does not just scale the Euclidean distance between two points.
We can however scale $x$\- and $y$-coordinates simultaneously by some factor
$s>0$. By a transformation of coordinates, we have for any set $X\subset
P\times P$ that

$\displaystyle\iiiint_{(p,q)\in
X}\left\|p-q\right\|\mathop{}\\!\mathrm{d}(p,q)$
$\displaystyle=\iiiint_{(p,q)\in
sX}\left\|p/s-q/s\right\|\left|\det(I_{4}/s)\right|\mathop{}\\!\mathrm{d}(p,q)$
$\displaystyle=\frac{1}{s^{5}}\iiiint_{(p,q)\in
sX}\left\|p-q\right\|\mathop{}\\!\mathrm{d}(p,q),$

where $I_{4}$ denotes the $4\times 4$ identity matrix. Using this, we can
obtain the following properties by scaling the final result suitably:

  1. 1.

By a rotation and scaling, the edge $e$ is contained in the $x$-axis and the
interior of $P$ is above $e$.

  2. 2.

By a horizontal translation, we have $x(p_{\ell})=0$, where $p_{\ell}$ is the
left pivot.

  3. 3.

By a scaling, we have $x(p_{r})=1$, where $p_{r}$ is the right pivot.

  4. 4.

By sliding the points $f$ and $g$ defining the top edge, we have $x(f)=1$ and
$x(g)=0$.

  5. 5.

By scaling the vectors $V_{0}$ and $V_{1}$, we have $y(V_{0})=y(V_{1})=1$.

Again, it is not always possible to obtain properties 3–5 due to some special
cases, but when it is not possible, we can get similar properties and proceed
in an analogous way.

We obtain a rotating trapezoid depending on the six parameters
$y(p_{\ell}),y(p_{r}),y(f),y(g),x(V_{0}),x(V_{1})$. Note that in the case of
the Beer-index, we had $y(f)=1$, but we do not see how to obtain that here
(while also maintaining the other properties). We can now rewrite the integral
(5) as

$\int_{0}^{1}\int_{h_{\ell}(\rho)}^{h_{r}(\rho)}\int_{0}^{T(\rho,h)}\int_{0}^{q}\sqrt{1+((1-\rho)\,x(V_{0})+\rho\,x(V_{1}))^{2}}\left(q-p\right)\left|\det
J_{F}\right|\mathop{}\\!\mathrm{d}p\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}h\mathop{}\\!\mathrm{d}\rho.$

Maple [24] is able to compute a closed-form expression of this integral,
involving only elementary functions. Unfortunately, when written as clear text
and using 3-letter names for the six parameters mentioned above, the formula
has a length of more than 12 million characters, so it is hardly of any
practical use. We conclude that using $O(n^{2})$ time, we can output
$O(n^{2})$ closed-form expressions of size $O(1)$ whose sum equals the
integral over $A$.

###  5.2 Shortest paths containing corners

Figure 12: Left: The regions $P_{-}[uv_{1}]$ and $P_{+}[uv_{2}]$ for chords
$uv_{1}$ and $uv_{2}$. The dashed arrows bound the wedge from $uv_{1}$
counterclockwise to $e_{-}$. Right: The sets $E_{-},E_{0},E_{+}$ for a corner
$u$.

For a corner $u$ of $P$, let $K_{u}\subset B$ be pairs of points $(p,q)$ such
that the first corner of $P$ on the shortest path $p\leadsto q$ from $p$ is
$u$. We then define

$M_{u}=\iiiint_{(p,q)\in K_{u}}\left\|p\leadsto
q\right\|\mathop{}\\!\mathrm{d}(p,q)$

and thus have

$\displaystyle\iiiint_{(p,q)\in B}\left\|p\leadsto
q\right\|\mathop{}\\!\mathrm{d}(p,q)=\sum_{u\in\mathcal{C}}M_{u},$ (6)

where $\mathcal{C}$ is the set of corners of $P$. We will compute the
contribution to the integral over $B$ for each set $K_{u}$. To this end, we
introduce some regions and integrals defined by diagonals of $P$. Consider a
chord $uv$ of $P$ where $u$ is a corner. Let $e_{-},e_{+}$ be the edges of $P$
preceding and succeeding $u$ in counterclockwise order, respectively; see
Figure 12 (left). We define $P_{-}[uv]\subset P$ to be the points $q$ in $P$
such that the shortest path from $u$ to $q$ starts with a segment in the
closed wedge from $uv$ counterclockwise to $e_{-}$. Similarly,
$P_{+}[uv]\subset P$ is the set of points $q$ in $P$ such that the shortest
path from $u$ to $q$ starts with a segment in the closed wedge from $uv$
clockwise to $e_{+}$. We also include $u$ in $P_{-}[uv]$ and $P_{+}[uv]$ to
get closed sets. We note that $P$ can be written as an interior-disjoint union
$P=P_{-}[uv]\cup P_{+}[uv]$. If $uv$ is a diagonal, we say that $P_{-}[uv]$
and $P_{+}[uv]$ are _diagonal regions_.

For $\pm\in\\{-,+\\}$, we now define

$L_{\pm}[uv]=\iint_{q\in P_{\pm}[uv]}\left\|u\leadsto
q\right\|\mathop{}\\!\mathrm{d}q.$

We will be explain in Section 5.3 how the areas $|P_{\pm}[uv]|$ and values
$L_{\pm}[uv]$ can be computed for all diagonals $uv$. Here we will assume that
it has already been done, and we will then see how the integral $M_{u}$ can be
computed. First, observe that if $u$ is a convex corner, then
$K_{u}=\textrm{vis}(u)\times\\{u\\}$, where $\textrm{vis}(u)$ is the
visibility polygon of $u$. This set has measure $0$ in $P\times P$, so
clearly, $M_{u}=0$ in this case. Assume therefore that $u$ is a concave
corner. Note that $K_{u}$ is a subset of $\textrm{vis}(u)\times P$. We now
partition $\textrm{vis}(u)$ into three parts as follows. Let $e^{\prime}_{-}$
and $e^{\prime}_{+}$ be the extensions of $e_{-}$ and $e_{+}$ from $u$,
respectively. Then $e^{\prime}_{-}$ and $e^{\prime}_{+}$ partition
$\textrm{vis}(u)$ into three parts $E_{-},E_{0},E_{+}$, such that $e_{-}$ is
on the boundary of $E_{-}$, $e_{+}$ is on the boundary of $E_{+}$, and $E_{0}$
is in between $e^{\prime}_{-}$ and $e^{\prime}_{+}$; see Figure 12 (right). We
can now partition $K_{u}$ into the three sets $K_{u}\cap(E_{-}\times P)$,
$K_{u}\cap(E_{0}\times P)$ and $K_{u}\cap(E_{+}\times P)$. We note that
$K_{u}\cap(E_{0}\times P)=E_{0}\times\\{u\\}$, which has measure $0$, so we
can express $M_{u}$ as

$M_{u}=\iiiint_{(p,q)\in K_{u}\cap(E_{-}\times P)}\left\|p\leadsto
q\right\|\mathop{}\\!\mathrm{d}(p,q)+\iiiint_{(p,q)\in K_{u}\cap(E_{+}\times
P)}\left\|p\leadsto q\right\|\mathop{}\\!\mathrm{d}(p,q).$

We describe how to evaluate the integral over $K_{u}\cap(E_{+}\times P)$, and
the one over $K_{u}\cap(E_{-}\times P)$ can be handled in an analogous way.

Figure 13: A triangle $\Lambda$ in $E_{+}$. The shortest paths from $p$ that
pass through $u$ have endpoints in
$\Delta(u\lambda^{\prime}_{2}u^{\prime}_{p})$,
$\Delta(uv^{\prime}\lambda^{\prime}_{2})$ or $P_{-}[uv]$.

The diagonals $uv$ and their extensions from $u$ partition $\textrm{vis}(u)$
into triangles; see Figure 13. For each triangle $\Lambda$ in $E_{+}$, we will
compute the contribution from paths from $p\in\Lambda$ to points $q$ such that
the first corner of $P$ on $p\leadsto q$ is $u$,

$M_{u}^{\Lambda}=\iiiint_{(p,q)\in K_{u}\cap(\Lambda\times P)}\left\|p\leadsto
q\right\|\mathop{}\\!\mathrm{d}(p,q).$

Taking the sum of this integral over all such triangles $\Lambda$ in $E_{+}$
will then give us the integral over $K_{u}\cap(E_{+}\times P)$.

For a point $p\in\textrm{vis}(u)$, let the extension of $pu$ from $u$ be
$uu^{\prime}_{p}$. We can then write

$M_{u}^{\Lambda}=\iint_{p\in\Lambda}\iint_{q\in
P_{-}[uu^{\prime}_{p}]}\left\|p\leadsto
q\right\|\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}p=\iint_{p\in\Lambda}\iint_{q\in
P_{-}[uu^{\prime}_{p}]}\left\|pu\right\|+\left\|u\leadsto
q\right\|\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}p.$

Let $uv$ be the diagonal which is the next after $uu^{\prime}_{p}$ in
counterclockwise direction around $u$ and let $vv^{\prime}$ be the extension
of $uv$ from $v$ (where we may have $v^{\prime}=v$). Note that this diagonal
$uv$ does not depend on the choice of $p\in\Lambda$, so we can partition
$P_{-}[uu^{\prime}_{p}]$ as
$P_{-}[uu^{\prime}_{p}]=P_{-}[uv]\cup\Delta(uv^{\prime}u^{\prime}_{p})$, and
we get

$\displaystyle M_{u}^{\Lambda}$
$\displaystyle=\iint_{p\in\Lambda}\left(\iint_{q\in
P_{-}[uv]}\left\|pu\right\|+\left\|u\leadsto
q\right\|\mathop{}\\!\mathrm{d}q+\iint_{q\in\Delta(uv^{\prime}u^{\prime}_{p})}\left\|pu\right\|+\left\|uq\right\|\mathop{}\\!\mathrm{d}q\right)\mathop{}\\!\mathrm{d}p$
$\displaystyle=\left|P_{-}[uv]\right|\iint_{p\in\Lambda}\left\|pu\right\|\mathop{}\\!\mathrm{d}p+\left|\Lambda\right|L_{-}[uv]+\iint_{p\in\Lambda}\iint_{q\in\Delta(uv^{\prime}u^{\prime}_{p})}\left\|pu\right\|+\left\|uq\right\|\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}p.$
(7)

Let the corners of $\Lambda$ be $u,\lambda_{1},\lambda_{2}$ in
counterclockwise order. Let $u\lambda^{\prime}_{1}$ and
$u\lambda^{\prime}_{2}$ be the extensions of $\lambda_{1}u$ and $\lambda_{2}u$
from $u$, so that we have an interior-disjoint union
$\Delta(uv^{\prime}u^{\prime}_{p})=\Delta(uv^{\prime}\lambda^{\prime}_{2})\cup\Delta(u\lambda^{\prime}_{2}u^{\prime}_{p})$.
We can rewrite the last integral as

$\iint_{p\in\Lambda}\iint_{q\in\Delta(uv^{\prime}\lambda^{\prime}_{2})}\left\|pu\right\|+\left\|uq\right\|\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}p+\iint_{p\in\Lambda}\iint_{q\in\Delta(u\lambda^{\prime}_{2}u^{\prime}_{p})}\left\|pu\right\|+\left\|uq\right\|\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}p.$

Here, the first part equals

$\displaystyle\left|\Delta(uv^{\prime}\lambda^{\prime}_{2})\right|\iint_{p\in\Lambda}\left\|pu\right\|\mathop{}\\!\mathrm{d}p+\left|\Lambda\right|\iint_{q\in\Delta(uv^{\prime}\lambda^{\prime}_{2})}\left\|uq\right\|\mathop{}\\!\mathrm{d}q,$
(8)

and the second equals

$\displaystyle\iint_{p\in\Lambda}\left\|pu\right\|\left|\Delta(u\lambda^{\prime}_{2}u^{\prime}_{p})\right|\mathop{}\\!\mathrm{d}p+\iint_{p\in\Lambda}\iint_{q\in\Delta(u\lambda^{\prime}_{2}u^{\prime}_{p})}\left\|uq\right\|\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}p.$
(9)

By changing order of integration, we obtain

$\displaystyle\iint_{p\in\Lambda}\iint_{q\in\Delta(u\lambda^{\prime}_{2}u^{\prime}_{p})}\left\|uq\right\|\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}p=$
$\displaystyle\iint_{q\in\Delta(u\lambda^{\prime}_{2}\lambda^{\prime}_{1})}\iint_{p\in\Delta(u\lambda_{1}u^{\prime}_{q})}\left\|uq\right\|\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}p=$
$\displaystyle\iint_{q\in\Delta(u\lambda^{\prime}_{2}\lambda^{\prime}_{1})}\left\|uq\right\|\left|\Delta(u\lambda_{1}u^{\prime}_{q})\right|\mathop{}\\!\mathrm{d}q.$
(10)

It therefore remains to explain how we compute each of the integrals

$\displaystyle\iint_{p\in\Lambda}\left\|pu\right\|\mathop{}\\!\mathrm{d}p,$
[in (7) and (8)] (11)
$\displaystyle\iint_{q\in\Delta(uv^{\prime}\lambda^{\prime}_{2})}\left\|uq\right\|\mathop{}\\!\mathrm{d}q,$
[in (8)] (12)
$\displaystyle\iint_{p\in\Lambda}\left\|pu\right\|\left|\Delta(uv^{\prime}u^{\prime}_{p})\right|\mathop{}\\!\mathrm{d}p,$
[in (9)] (13)
$\displaystyle\iint_{q\in\Delta(u\lambda^{\prime}_{2}\lambda^{\prime}_{1})}\left\|uq\right\|\left|\Delta(u\lambda_{1}u^{\prime}_{q})\right|\mathop{}\\!\mathrm{d}q.$
[in (10)] (14)

For three points $a,b,c$, define the integral $\Xi(abc)$ over the triangle
$\Delta(abc)$ as

$\Xi(abc)=\iint_{q\in\Delta(abc)}\left\|aq\right\|\mathop{}\\!\mathrm{d}q.$

We then note that (11) is just $\Xi(u\lambda_{1}\lambda_{2})$ and (12) is
$\Xi(uv^{\prime}\lambda^{\prime}_{2})$. Furthermore, note that (13) and (14)
have the same form, so that a closed-form expression for one automatically
gives one for the other as well. In Appendix B, we show how the integrals
(11)–(14) can be computed using Maple.

Given that we know all diagonals $uv$ from $u$, it takes $O(n)$ time to
compute the contribution $M_{u}$ by iterating over all triangles $\Lambda$ in
$E_{+}$ and $E_{-}$. We therefore use $O(n^{2})$ time in total on outputting
closed-form expressions, each of size $O(1)$, whose sum equals the integral
over $B$.

###  5.3 Computing the values $|P_{\pm}[uv]|$ and $L_{\pm}(uv)$

Figure 14: The figures show how a diagonal region $P_{-}[uv]$ can be described
as interior-disjoint unions of smaller diagonal regions and triangles. This
depends on which regions $P_{-}[uv]$ and $P_{+}[uv]$ contain the edges
incident at $v$. The region $P_{-}[uv]$ is colored and $P_{+}[uv]$ is white.
Left: One edge is in $P_{+}[uv]$. Middle: Both edges are in $P_{+}[uv]$.
Right: Both edges are in $P_{-}[uv]$.

We explain how the values $|P_{-}[uv]|$ and $L_{-}(uv)$ are computed; the
values $|P_{+}[uv]|$ and $L_{+}(uv)$ are computed in an analogous way. We
first observe that each diagonal region $P_{-}[uv]$ can be described as an
interior-disjoint union of at most two smaller diagonal regions and at most
two triangles. Let $vv^{\prime}$ be the extension of $uv$ from $v$. We now
have two cases.

  1. 1.

One or both of the edges of $P$ incident at $v$ are in $P_{+}[uv]$; see Figure
14 (left and middle). Let $uv_{0}$ be the diagonal after $uv$ in
counterclockwise order around $u$, and let $v_{0}v_{0}^{\prime}$ be the
extension of $uv_{0}$ from $v_{0}$. Then

$\displaystyle P_{-}[uv]=P_{-}[uv_{0}]\cup\Delta(uv^{\prime}_{0}v^{\prime}).$
(15)

  2. 2.

Both edges of $P$ incident at $v$ are in $P_{-}[uv]$; see Figure 14 (right).
Let $uv_{0}$ be the diagonal after $uv$ in counterclockwise order around $u$,
and let $v_{0}v_{0}^{\prime}$ be the extension of $uv_{0}$ from $v_{0}$. Let
$vw$ be the diagonal after $vv^{\prime}$ in counterclockwise order around $v$,
and let $ww^{\prime}$ be the extension of $vw$ from $w$. Then

$\displaystyle P_{-}[uv]=P_{-}[uv_{0}]\cup\Delta(uv^{\prime}_{0}v)\cup
P_{-}[vw]\cup\Delta(vw^{\prime}v^{\prime}).$ (16)

These formulas naturally lead to a recursive approach for computing the values
$|P_{-}[uv]|$. Recall that we can compute all diagonals of $P$ in $O(n^{2})$
time. For each diagonal $uv$ we then compute the area $|P_{-}[uv]|$ using (15)
or (16). When a value for a smaller region is needed which we have not yet
computed, we do it recursively. We use $O(1)$ time per diagonal $uv$ to
compute the area $|P_{-}[uv]|$.

Using the notation $\Xi(\cdot)$ introduced in Section 5.2 for integrals of
distances over triangles, we get the following two formulas corresponding to
(15) and (16).

$\displaystyle L_{-}[uv]$ $\displaystyle=L_{-}[uv_{0}]+\Xi(uv^{\prime}_{0}v),$
(17) $\displaystyle L_{-}[uv]$
$\displaystyle=L_{-}[uv_{0}]+\Xi(uv^{\prime}_{0}v)+L_{-}[vw]+\Xi(vw^{\prime}v^{\prime})+(\left|P_{-}[vw]\right|+\left|\Delta(vw^{\prime}v^{\prime})\right|)\left\|uv\right\|.$
(18)

We can therefore likewise compute all the values $L_{-}[uv]$ recursively. In
total, closed-form expressions for all the values $|P_{\pm}[uv]|$ and
$L_{\pm}[uv]$ can be returned in $O(n^{2})$ time.

Figure 15: A square polygon $P$ where the segment $c_{1}c_{2}$ is a
(degenerate) hole. The curve $\gamma$ is the locus of points $g$ where the
paths from $p$ to $g$ above $c_{1}$ and below $c_{2}$ are equally long. Each
such point is an intersection of circles $C_{1}$ and $C_{2}$ centered at
$c_{1}$ and $c_{2}$, where the radius of $C_{1}$ is $\|pc_{2}\|-\|pc_{1}\|$
larger than the radius of $C_{1}$. Therefore, $\gamma$ is an arc on a
hyperbola.

##  6 Concluding remarks

It is interesting whether the expected distance between two points in a
polygon with holes can also be expressed in closed form. This seems much
harder. Even the expected distance from a fixed point $p$ to another random
point is complicated to compute in the presence of holes; see Figure 15.
Whether the shortest path from $p$ to another point $q$ goes above or below a
single hole sometimes depends on whether $q$ is above or below an arc $\gamma$
on a certain hyperbola which depends on $p$.

Another interesting question is whether our $O(n^{2})$-time algorithms for the
Beer-index and the expected $L_{2}$-distance are conditionally optimal under
some reasonable hypothesis. We could not find faster algorithms, even for
simple polygons when computing the Beer-index and convex polygons when
computing the expected distance.

## References

  * [1] Pankaj K. Agarwal  “Parititoning Arrangements of Lines II: Applications”  In _Discret. Comput. Geom._ 5, 1990, pp. 533–573  DOI: 10.1007/BF02187809
  * [2] Takao Asano, Tetsuo Asano, Leonidas J. Guibas, John Hershberger and Hiroshi Imai  “Visibility of Disjoint Polygons”  In _Algorithmica_ 1.1, 1986, pp. 49–63  DOI: 10.1007/BF01840436
  * [3] S. Bandara and S.C. Wirasinghe  “Walking distance minimization for airport terminal configurations”  In _Transportation Research Part A: Policy and Practice_ 26.1, 1992, pp. 59–74  DOI: https://doi.org/10.1016/0965-8564(92)90045-9
  * [4] Gerald Beer  “The index of convexity and the visibility function”  In _Pacific Journal of Mathematics_ 44.1 Mathematical Sciences Publishers, 1973, pp. 59–67  DOI: 10.2140/pjm.1973.44.59
  * [5] Mark Berg, Otfried Cheong, Marc J. Kreveld and Mark H. Overmars  “Computational geometry: algorithms and applications, 3rd Edition”  Springer, 2008  URL: https://www.worldcat.org/oclc/227584184
  * [6] Gilles Bonnet, Anna Gusakova, Christoph Thäle and Dmitry Zaporozhets  “Sharp inequalities for the mean distance of random points in convex bodies”  In _Advances in Mathematics_ 386, 2021, pp. 1–27  DOI: https://doi.org/10.1016/j.aim.2021.107813
  * [7] Manuel Bronstein  “Symbolic Integration I: Transcendental Functions”, Algorithms and Computation in Mathematics (AACIM, volume 1), 2005  DOI: 10.1007/978-3-662-03386-9
  * [8] Kevin Buchin, Bram Custers, Ivor Hoog, Maarten Löffler, Aleksandr Popov, Marcel Roeloffzen and Frank Staals  “Segment Visibility Counting Queries in Polygons”  arXiv, 2022  DOI: 10.48550/ARXIV.2201.03490
  * [9] Kevin Buchin, Irina Kostitsyna, Maarten Löffler and Rodrigo I. Silveira  “Region-Based Approximation of Probability Distributions (for Visibility Between Imprecise Points Among Obstacles)”  In _Algorithmica_ 81.7, 2019, pp. 2682–2715  DOI: 10.1007/s00453-019-00551-2
  * [10] Uwe Bäsel  “The moments of the distance between two random points in a regular polygon”  arXiv, 2021  DOI: 10.48550/ARXIV.2101.03815
  * [11] Sergio Cabello  “Subquadratic Algorithms for the Diameter and the Sum of Pairwise Distances in Planar Graphs”  In _ACM Trans. Algorithms_ 15.2, 2019, pp. 21:1–21:38  DOI: 10.1145/3218821
  * [12] Bernard Chazelle  “Triangulating a Simple Polygon in Linear Time”  In _Discret. Comput. Geom._ 6, 1991, pp. 485–524  DOI: 10.1007/BF02574703
  * [13] Emanuel Czuber  “Geometrische Wahrscheinlichkeiten und Mittelwerte”, 1884  URL: https://books.google.de/books/about/Geometrische_Wahrscheinlichkeiten_und_Mi.html?id=2jQAAAAAQAAJ
  * [14] Herbert Edelsbrunner and Ernst Peter Mücke  “Simulation of Simplicity: A Technique to Cope with Degenerate Cases in Geometric Algorithms”  In _ACM Trans. Graph._ 9.1, 1990, pp. 66–104  DOI: 10.1145/77635.77639
  * [15] Pawel Gawrychowski, Haim Kaplan, Shay Mozes, Micha Sharir and Oren Weimann  “Voronoi Diagrams on Planar Graphs, and Computing the Diameter in Deterministic $\tilde{O}(n^{5/3})$ Time”  In _Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms (SODA 2018)_ , 2018, pp. 495–514  DOI: 10.1137/1.9781611975031.33
  * [16] Mukulika Ghosh, Nancy M. Amato, Yanyan Lu and Jyh-Ming Lien  “Fast approximate convex decomposition using relative concavity”  In _Comput. Aided Des._ 45.2, 2013, pp. 494–504  DOI: 10.1016/j.cad.2012.10.032
  * [17] Leonidas J. Guibas, John Hershberger, Daniel Leven, Micha Sharir and Robert Endre Tarjan  “Linear-Time Algorithms for Visibility and Shortest Path Problems Inside Triangulated Simple Polygons”  In _Algorithmica_ 2, 1987, pp. 209–233  DOI: 10.1007/BF01840360
  * [18] John Hershberger  “An Optimal Visibility Graph Algorithm for Triangulated Simple Polygons”  In _Algorithmica_ 4.1, 1989, pp. 141–155  DOI: 10.1007/BF01553883
  * [19] John Hershberger and Subhash Suri  “A Pedestrian Approach to Ray Shooting: Shoot a Ray, Take a Walk”  In _J. Algorithms_ 18.3, 1995, pp. 403–431  DOI: 10.1006/jagm.1995.1017
  * [20] Arthur C. Hsu  “The expected distance between two random points in a polygon”, 1990  URL: https://dspace.mit.edu/bitstream/handle/1721.1/14232/22504138-MIT.pdf
  * [21] Rolf Klein  “Algorithmische Geometrie”  Addison-Wesley-Longman, 1997
  * [22] Der-Tsai Lee  “Proximity and reachability in the plane”, 1978
  * [23] Jyh-Ming Lien and Nancy M. Amato  “Approximate convex decomposition of polygons”  In _Comput. Geom._ 35.1-2, 2006, pp. 100–123  DOI: 10.1016/j.comgeo.2005.10.005
  * [24] Maple 2019.2, Maplesoft, a division of Waterloo Maple Inc.,ΩWaterloo, Ontario. Maple is a trademark of Waterloo Maple Inc.
  * [25] Günter Rote  “The degree of convexity”  In _Abstracts of the 29th European Workshop on Computational Geometry (EuroCG 2013)_ , 2013, pp. 69–72  URL: http://page.mi.fu-berlin.de/rote/Papers/pdf/The+degree+of+convexity.pdf
  * [26] Günter Rote  “Note after publication of [25]”, 2017  URL: https://page.mi.fu-berlin.de/rote/Papers/abstract/The+degree+of+convexity.html
  * [27] Micha Sharir and Amir Schorr  “On Shortest Paths in Polyhedral Spaces”  In _SIAM J. Comput._ 15.1, 1986, pp. 193–215  DOI: 10.1137/0215014
  * [28] Helman Stern  “Polygonal entropy: A convexity measure”  In _Pattern Recognit. Lett._ 10.4, 1989, pp. 229–235  DOI: 10.1016/0167-8655(89)90093-7
  * [29] Vojin Tošić  “A review of airport passenger terminal operations analysis and modelling”  In _Transportation Research Part A: Policy and Practice_ 26.1, 1992, pp. 3–26  DOI: https://doi.org/10.1016/0965-8564(92)90041-5
  * [30] Emo Welzl  “Constructing the Visibility Graph for $n$-Line Segments in $O(n^{2})$ Time”  In _Inf. Process. Lett._ 20.4, 1985, pp. 167–171  DOI: 10.1016/0020-0190(85)90044-4

##  Appendix A Counterexample to a claim from [9]

The algorithm for computing the Beer-index of a polygon with holes from [9]
works as follows. We consider a convex polygon $P$ and a family
$H_{1},\ldots,H_{k}$ of obstacles, and we want to find the probability that
for two random points $p,q$ in $P\setminus\bigcup H_{i}$, the line segment
$pq$ is disjoint from all obstacles. Note that this is just a slight
reformulation of the problem of computing the Beer-index of a polygon with
holes.

We partition $P\setminus\bigcup H_{i}$ into trapezoids $T_{1},\ldots,T_{t}$.
We now compute the Beer-index as the sum of two entities: (i) the probability
that two random points $p,q$ are in the same trapezoid, and (ii) the
probability that two random points $p,q$ are in different trapezoids and the
segment $pq$ is disjoint from all obstacles. For the second case, we compute
the contribution for each pair $T_{i},T_{j}$ of trapezoids and sum over all
these pairs.

Let us consider a fixed pair $T_{i},T_{j}$ of trapezoids. For each pair
$s_{1}^{i},s_{2}^{i}$ of two segments of $T_{i}$ and each pair
$s_{1}^{j},s_{2}^{j}$ of two segments of $T_{j}$, we do as follows. Using
geometric duality, we map each line intersecting all four segments
$s_{1}^{i},s_{2}^{i},s_{1}^{j},s_{2}^{j}$ to a point in the dual space, and
these points together form a convex polygon $L^{*}$. Likewise, the dual points
of the set of lines intersecting an obstacle $H_{l}$ form an “hourglass-
shaped” unbounded region $H_{l}^{*}$ in the dual space. We now compute the
overlay of $L^{*}$ with all the regions $H_{l}^{*}$. This results in an
arrangement $\mathcal{S}^{*}$ of the plane. Each cell $C$ of $\mathcal{S}^{*}$
contained in $L^{*}$ corresponds to a set of lines in primal space that all
cross the four segments $s_{1}^{i},s_{2}^{i},s_{1}^{j},s_{2}^{j}$.
Furthermore, the visibility from $T_{i}$ to $T_{j}$ is either not blocked
along any of these lines or blocked for all the lines by one or more
obstacles. It is then described in the paper how to compute the contribution
to the Beer-index of all pairs of points $p\in T_{1}$, $q\in T_{2}$, where the
dual point $l(pq)^{*}$ is in $C$. The time it takes to compute this
contribution is proportional to the complexity of $C$.

Running over all cells $C$ contained in $L^{*}$, we obtain that computing the
contribution for all pairs of points $p\in T_{1}$, $q\in T_{2}$, where $l(pq)$
intersects $s_{1}^{i},s_{2}^{i},s_{1}^{j},s_{2}^{j}$, can be done in
$O(n^{2})$ time, where $n$ is the total complexity of $P$ and the obstacles.
Since there are $O(1)$ choices of the four segments
$s_{1}^{i},s_{2}^{i},s_{1}^{j},s_{2}^{j}$, the contribution from all pairs
$p\in T_{1}$, $q\in T_{2}$ can then also be computed in time $O(n^{2})$.

It is then erroneously claimed in Lemma 3 that the total size of the
arrangements $\mathcal{S}^{*}$, as we sum over all pairs of trapezoids
$T_{i},T_{j}$, is also $O(n^{2})$. This is not true. For $\Omega(n^{2})$ of
the pairs, the complexity may be $\Omega(n^{2})$, resulting in a total
complexity of $\Omega(n^{4})$. This can even be the case for simple polygons,
as the polygon in Figure 16 shows. The polygon is simply a rectangle with some
spikes removed along the top edge. The obstacles will therefore be the removed
spikes. Note that a line can intersect any consecutive interval of obstacles
and no other obstacle, so there are $\Omega(n^{2})$ subsets of obstacles that
can be intersected. The arrangement we get by overlaying the hour-glasses
$H_{l}^{*}$ therefore has complexity $\Omega(n^{2})$. It then follows that for
$\Omega(n^{2})$ pairs of trapezoids $T_{i},T_{j}$, the arrangement
$\mathcal{S}^{*}$ has complexity $\Omega(n^{2})$, resulting in a running time
of $\Omega(n^{4})$.

Figure 16: For polygons of this type, the algorithm from [9] has running time
$\Omega(n^{4})$.

##  Appendix B Triangle integrals

We now explain how to compute the integral

$\Xi(abc)=\iint_{q\in\Delta(abc)}\left\|aq\right\|\mathop{}\\!\mathrm{d}q.$

By rotation and scaling, we can assume that $a=(0,0)$ and $y(b)=y(c)=1$, where
$b$ is to the left of $c$. We make a transformation of coordinates so that the
usual coordinates corresponding to the pair $(\rho,q)$ is

$F(\rho,q)=q\cdot(b\cdot(1-\rho)+c\cdot\rho)=q\cdot(x(b)\cdot(1-\rho)+x(c)\cdot\rho,1).$

We then have the Jacobian

$J_{F}=\begin{bmatrix}q\cdot(-x(b)+x(c))&x(b)\cdot(1-\rho)+x(c)\cdot\rho\\\
0&1\end{bmatrix},$

with the determinant $\det J_{F}=q\cdot(x(c)-x(b))$. As $x(b)<x(c)$ and
$q\in[0,1]$, we have $|\det J_{F}|=q\cdot(x(c)-x(b))$.

We then have

$\Xi(abc)=\int_{0}^{1}\int_{0}^{1}q\,\sqrt{(x(b)\cdot(1-\rho)+x(c)\cdot\rho)^{2}+1}\left|\det
J_{F}\right|\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}\rho.$

Using Maple [24], we can find a find a formula for this integral.

Figure 17: The setting of integrals (13) and (14).

We now turn our attention to the integrals (13) and (14). We focus on (13)
since (14) is analogous; see Figure 17. The corners of $\Lambda$ are
$u,\lambda_{1},\lambda_{2}$, and by rotation, scaling and translation, we can
assume that $u=(0,0)$, $y(\lambda_{1})=y(\lambda_{2})=1$ and
$x(\lambda_{2})<x(\lambda_{1})$. Recall that $u\lambda_{1}^{\prime}$ and
$u\lambda_{2}^{\prime}$ are the extensions of $\lambda_{1}u$ and
$\lambda_{2}u$ from $u$. Let the edge containing $\lambda_{1}^{\prime}$ and
$\lambda_{2}^{\prime}$ be defined by corners $f$ and $g$. Unless $fg$ is
vertical, we can assume $x(f)=1$ and $x(g)=0$ by sliding the points on the
line containing the edge. We can then express the integral in terms of the
four parameters $x(\lambda_{1}),x(\lambda_{2}),y(f),y(g)$.

We use the same transformation of coordinates to pairs $(\rho,q)$ as used
above for the integral $\Xi(abc)$. We rewrite (13) as the integral

$\int_{0}^{1}\int_{0}^{1}q\,\sqrt{(x(\lambda_{1})\cdot(1-\rho)+x(\lambda_{1})\cdot\rho)^{2}+1}\left|\Delta(\rho)\right|\left|\det
J_{F}\right|\mathop{}\\!\mathrm{d}q\mathop{}\\!\mathrm{d}\rho,$

where $\Delta(\rho)$ is shorthand for the triangle corresponding to
$\Delta(u\lambda^{\prime}_{2}u^{\prime}_{p})$ (note that using the transformed
coordinates $(\rho,q)$, the triangle depends only on $\rho$).

As it turns out, Maple refuses to evaluate this integral directly, but we can
evaluate the inner integral and we observe that the result has the form

$\frac{(a\,\rho+b)\,\sqrt{c\,\rho^{2}+d\,\rho+e}}{f\,\rho+g},$

where $a,b,c,d,e,f$ are constants depending on the parameters
$x(\lambda_{1}),x(\lambda_{2}),y(f),y(g)$. Maple is able to find an anti-
derivative of this function. By inserting the expressions for $a,b,c,d,e,f$,
we can then evaluate the integral (13).

##  Appendix C Maple calculation

The following code shows a calculation of a closed-form expression for the
integral (1).