Department of Computer Science, University of Maryland, College Park, USA and
~octavo@umd.eduhttps://orcid.org/0000-0002-5704-312X Wellesley College,
Wellesley, MA and<EMAIL_ADDRESS>Montgomery Blair High School, Silver
Spring, MD and<EMAIL_ADDRESS>Indiana University Boolmington,
Bloomington, IN and<EMAIL_ADDRESS>Stony Brook University, Stony Brook, NY
and<EMAIL_ADDRESS>Department of Computer Science,
University of Maryland, College Park, USA and
https://www.cs.umd.edu/~mount/mount@umd.eduhttps://orcid.org/0000-0002-3290-8932
Auguste H. Gezalyan, Soo H. Kim, Carlos Lopez, Daniel Skora, Zofia
Stefankovic, and David M. Mount, [500]Theory of computation Computational
geometry

# Delaunay Triangulations in the Hilbert Metric

Auguste H. Gezalyan  Soo H. Kim  Carlos Lopez  Daniel Skora  Zofia Stefankovic
David M. Mount

###### Abstract

The Hilbert metric is a distance function defined for points lying within the
interior of a convex body. It arises in the analysis and processing of convex
bodies, machine learning, and quantum information theory. In this paper, we
show how to adapt the Euclidean Delaunay triangulation to the Hilbert geometry
defined by a convex polygon in the plane. We analyze the geometric properties
of the Hilbert Delaunay triangulation, which has some notable differences with
respect to the Euclidean case, including the fact that the triangulation does
not necessarily cover the convex hull of the point set. We also introduce the
notion of a Hilbert ball at infinity, which is a Hilbert metric ball centered
on the boundary of the convex polygon. We present a simple randomized
incremental algorithm that computes the Hilbert Delaunay triangulation for a
set of $n$ points in the Hilbert geometry defined by a convex $m$-gon. The
algorithm runs in $O(n(\log n+\log^{3}m))$ expected time. In addition we
introduce the notion of the Hilbert hull of a set of points, which we define
to be the region covered by their Hilbert Delaunay triangulation. We present
an algorithm for computing the Hilbert hull in time $O(nh\log^{2}m)$, where
$h$ is the number of points on the hull’s boundary.

###### keywords:

Delaunay Triangulations, Hilbert metric, convexity, randomized algorithms

##  1 Introduction

David Hilbert introduced the Hilbert metric in 1895 [15]. Given any convex
body $\Omega$ in $d$-dimensional space, the Hilbert metric defines a distance
between any pair of points in the interior of $\Omega$ (see Section 2 for
definitions). The Hilbert metric has a number of useful properties. It is
invariant under projective transformations, and straight lines are geodesics.
When $\Omega$ is a Euclidean ball, it realizes the Cayley-Klein model of
hyperbolic geometry. When $\Omega$ is a simplex, it provides a natural metric
over discrete probability distributions (see, Nielsen and Sun [20, 21]). An
excellent resource on Hilbert geometries is the handbook of Hilbert geometry
by Papadopoulos and Troyanov [22].

The Hilbert geometry provides new insights into classical questions from
convexity theory. Efficient approximation of convex bodies has a wide range of
applications, including approximate nearest neighbor searching both in
Euclidean space [6] and more general metrics [1], optimal construction of
$\varepsilon$-kernels [4], solving the closest vector problem approximately
[11, 24, 12, 18], and computing approximating polytopes with low combinatorial
complexity [5, 3]. These works all share one thing in common—they approximate
a convex body by covering it with elements that behave much like metric balls.
These covering elements go under various names: Macbeath regions, Macbeath
ellipsoids, Dikin ellipsoids, and $(2,\varepsilon)$-covers. Vernicos and Walsh
showed that these shapes are, up to a constant scaling factor, equivalent to
Hilbert balls [26, 2]. In addition, the Hilbert metric behaves nicely in the
context flag approximability of convex polytopes as studied by Vernicos and
Walsh [27].

Other applications of the Hilbert metric include machine learning [20],
quantum information theory [23], real analysis [17], and optimal mass
transport [9]. Despite its obvious appeals, only recently has there been any
work on developing classical computational geometry algorithms that operate in
the Hilbert metric. Nielsen and Shao characterized balls in the Hilbert metric
defined by a convex polygon with $m$ sides [19]. Hilbert balls are convex
polygons bounded by $2m$ sides. Nielsen and Shao showed that Hilbert balls can
be computed in $O(m)$ time, and they developed dynamic software for generating
them. Gezalyan and Mount presented an $O(mn\log n)$ time algorithm for
computing the Voronoi diagram of $n$ point sites in the Hilbert polygonal
metric [13] (see Figure 1(a) and (b)). They showed that the diagram has worst-
case combinatorial complexity of $O(mn)$. Bumpus et al. [8] further analyzed
the properties of balls in the Hilbert metric and presented software for
computing Hilbert Voronoi diagrams.

Figure 1: (a) A convex polygon $\Omega$ and sites, (b) the Voronoi diagram,
and (c) the Delaunay triangulation.

In this paper, we present an algorithm for computing the Delaunay
triangulation of a set $P$ of $n$ point sites in the Hilbert geometry defined
by an $m$-sided convex polygon $\Omega$. The Hilbert Delaunay triangulation is
defined in the standard manner as the dual of the Hilbert Voronoi diagram (see
Figure 1(c)). Our algorithm is randomized and runs in $O(n(\log n+\log^{3}m))$
expected time. Excluding the polylogarithmic term in $m$, this matches the
time of the well-known randomized algorithm for Euclidean Delaunay
triangulations [14]. It is significantly more efficient than the time to
compute the Hilbert Voronoi diagram, which is possible because of the smaller
output size. A central element of our algorithm is an $O(\log^{3}m)$ time
algorithm for computing Hilbert circumcircles.

Unlike the Euclidean case, the Delaunay triangulation does not necessarily
triangulate the convex hull of $P$. We also provide a characterization of when
three points admit a Hilbert ball whose boundary contains all three points,
and introduce the notion of the Hilbert hull, which we define to be the region
covered by the Hilbert Delaunay triangulation. We give an algorithm for the
Hilbert hull that runs in time $O(nh\log^{2}m)$, where $h$ is the number of
points on the Hilbert hull’s boundary.

##  2 Preliminaries

###  2.1 The Hilbert Metric and Hilbert Balls

The Hilbert metric is defined over the interior of a convex body $\Omega$ in
$\mathbb{R}^{d}$ (that is, a closed, bounded, full dimensional convex set).
Let $\partial\kern 1.0pt\Omega$ denote $\Omega$’s boundary. Unless otherwise
stated, we assume throughout that $\Omega$ is an $m$-sided convex polygon in
$\mathbb{R}^{2}$. Given two points $p$ and $q$ in $\Omega$, let
$\overline{pq}$ denote the _chord_ of $\Omega$ defined by the line through
these points.

Figure 2: (a) The Hilbert distance $d_{\Omega}(p,q)$, (b) spokes defined by
$p$, and (c) the Hilbert ball $B(p,\rho)$.

######  Definition 2.1 (Hilbert metric).

Given a convex body $\Omega$ in $\mathbb{R}^{d}$ and two distinct points
$p,q\in\operatorname{int}(\Omega)$, let $p^{\prime}$ and $q^{\prime}$ denote
endpoints of the chord $\overline{pq}$, so that the points are in the order
$\langle p^{\prime},p,q,q^{\prime}\rangle$. The _Hilbert distance_ between $p$
and $q$ is defined

$d_{\Omega}(p,q)\leavevmode\nobreak\ =\leavevmode\nobreak\
\frac{1}{2}\ln\left(\frac{\|q-p^{\prime}\|}{\|p-p^{\prime}\|}\frac{\|p-q^{\prime}\|}{\|q-q^{\prime}\|}\right),$

and define $d_{\Omega}(p,p)=0$ (see Figure 2(a)).

Note that the quantity in the logarithm is the cross ratio
$(q,p;p^{\prime},q^{\prime})$. Since cross ratios are preserved by projective
transformations, it follows that the Hilbert distances are invariant under
projective transformations. Straight lines are geodesics, but not all
geodesics are straight lines. The Hilbert distance satisfies all the axioms of
a metric, and in particular it is symmetric and the triangle inequality holds
[22]. When $\Omega$ is a probability simplex [7], this symmetry distinguishes
it from other common methods of calculating distances between probability
distributions, such as the Kullback-Leibler divergence [16].

Given $p\in\operatorname{int}(\Omega)$ and $\rho>0$, let $B(p,\rho)$ denote
the Hilbert ball of radius $\rho$ centered at $p$. Nielsen and Shao showed how
to compute Hilbert balls [19]. Consider the set of $m$ chords $\overline{pv}$
for each vertex $v$ of $\Omega$ (see Figure 2(b)). These are called the
_spokes_ defined by $p$. For each spoke $\overline{pv}$, consider the two
points at Hilbert distance $\rho$ on either side of $p$. $B(p,\rho)$ is the
(convex) polygon defined by these $2m$ points (see Figure 2(c)). Given any
line that intersects $\Omega$, a simple binary search makes it possible to
determine the two edges of $\Omega$’s boundary intersected by the line.
Applying this to the line passing through $p$ and $q$, it follows that Hilbert
distances can be computed in $O(\log m)$ time.

###  2.2 The Hilbert Voronoi Diagram

Given a set $P$ of $n$ point sites in $\operatorname{int}(\Omega)$, the
_Hilbert Voronoi diagram_ of $P$ is defined in the standard manner as a
subdivision of $\operatorname{int}(\Omega)$ into regions, called _Voronoi
cells_ , based on which site of $P$ is closest in the Hilbert distance. It was
shown by Gezalyan and Mount [13] that each Voronoi cell is star-shaped with
respect to its site. Given two sites $p,q\in P$, the _$(p,q)$ -bisector_ is
the set of points that are equidistant from both sites in the Hilbert metric
(see Figure 3(a)).

Figure 3: The $(p,q)$-bisector is a piecewise conic with joints on the spokes
of $p$ and $q$.

The distances from any point $x$ on the $(p,q)$-bisector to $p$ and $q$ are
determined by the edges of $\partial\kern 1.0pt\Omega$ incident to the chords
$\overline{px}$ and $\overline{qx}$. We can subdivide the interior $\Omega$
into equivalence classes based upon the identity of these edges. It is easy to
see that as $x$ crosses a spoke of either $p$ or $q$, it moves into a
different equivalence class. It has been shown that within each equivalence
class, the bisector is a conic [13, 8]. It follows that the $(p,q)$-bisector
is a piecewise conic, whose joints lie on sector boundaries (see Figure 3(b)).
It is also known that the worst-case combinatorial complexity of the Hilbert
Voronoi diagram is $\Theta(mn)$ [13].

###  2.3 The Hilbert Delaunay Triangulation

The _Hilbert Delaunay triangulation_ of $P$, denoted $\operatorname{DT}(P)$,
is defined as the dual of the Hilbert Voronoi diagram, where two sites $p$ and
$q$ are connected by an edge if their Voronoi cells are adjacent. Throughout,
we make the general-position assumption that no four sites of $P$ are Hilbert
equidistant from a single point in $\operatorname{int}(\Omega)$. It is easy to
see that familiar concepts from Euclidean Delaunay triangulations apply, but
using Hilbert balls rather than Euclidean balls. For example, two sites are
adjacent in the triangulation if and only if there exists a Hilbert ball whose
boundary contains both sites, and whose interior contains no sites. (The
center of this ball lies on the Voronoi edge between the sites.) Also, three
points define a triangle in $DT(P)$ if and only if there is a Hilbert ball
whose boundary contains all three points, but is otherwise empty.

Consider a triangle $\triangle pqr$ in $\Omega$’s interior. A Hilbert ball
whose boundary passes through all three points is said to _circumscribe_ this
triangle. As we shall see in Section 4, some triangles admit no circumscribing
Hilbert ball, but the following lemma shows that if it does exist, then it is
unique. We refer to the boundary of such a ball as a _Hilbert circumcircle_.

######  Lemma 2.2.

There is at most one Hilbert circumcircle for any triangle in $\Omega$’s
interior.

######  Proof 2.3.

Let $p$, $q$, and $r$ be the vertices of the triangle, which we assume are not
collinear. Suppose for the sake of a contradiction there did exist two
circumscribing Hilbert balls $B_{1}$ and $B_{2}$. Let $c_{1}$ and $c_{2}$
denote their respective centers, and let $E=\Omega\setminus(B_{1}\cup B_{2})$
denote the region exterior to both balls. Since Hilbert balls are convex, $p$,
$q$, and $r$ all lie on the boundary of $E$. It follows that for any point $w$
in $E$, there exist three nonintersecting curves joining $v$ to $p$, $q$, and
$r$. Since both balls are empty, points $c_{1}$ and $c_{2}$ are vertices in
the Voronoi diagram of $\\{p,q,r\\}$. Since Voronoi cells are star-shaped, the
six open line segments joining each center to each point do not intersect.
Further, these lie within $B_{1}\cup B_{2}$, and so they are disjoint from the
curves joining $w$ to these points. Therefore the complete bipartite graph
$\\{p,q,r\\}\times\\{c_{1},c_{2},w\\}$ is planar, but this is impossible since
it would imply a planar embedding of $K_{3,3}$ (the complete bipartite
$3\times 3$ graph).

Note that the non-collinear assumption is necessary. To see why, let $\Omega$
be an axis aligned rectangle, and let $\ell$ be a horizontal line through the
rectangle’s center. It is easy to construct two Hilbert balls, one above the
line and one below, whose boundaries contain a segment along this line.
Placing the three points within this segment would violate the lemma.

The following lemma shows that $\operatorname{DT}(P)$ is a connected, planar
graph. Its proof is similar to the Euclidean case.

######  Lemma 2.4.

Given any discrete set of points $P$, $\operatorname{DT}(P)$ is a planar
(straight-line) graph that spans $P$.

######  Proof 2.5.

Suppose to the contrary that there existed $a,b,c,d\in P$ such that the edges
$(a,b)$ and $(c,d)$, viewed as open line segments, intersected each other.
These edges are witnessed by two empty Hilbert balls, $B_{1}$ for $(a,b)$ and
$B_{2}$ for $(c,d)$ (see Figure 4(a)). By convexity, $B_{1}$ contains the
segment $ab$, and by emptiness, its interior does not contain $c$ or $d$.
Similarly, $B_{2}$ contains the segment $cd$, and its interior does not
contain $a$ or $b$. and $B_{2}$. It follows the boundaries of $B_{1}$ and
$B_{2}$ must intersect in at least four points. Applying Lemma 2.2 to any
three of these four intersection points yields a contradiction.

Figure 4: Proof of Lemma 2.4.

To prove the $\operatorname{DT}(P)$ spans $P$, suppose to the contrary that
this was not so. Consider the two sites $p,q\in P$ that are in different
components but closest in terms of their Hilbert distance. Let $c$ denote the
center of the smallest Hilbert ball $B$ containing $p$ and $q$. Because
straight lines are geodesics, $c$ lies on the line segment $pq$, and
$d_{\Omega}(p,q)=d_{\Omega}(p,c)+d_{\Omega}(c,q)$ (see Figure 4(b)). We assert
that this ball can contain no other site in its interior. If there were such a
site $r$, it lies in a different connected component from $p$ or $q$ (possibly
both). Let us assume without loss of generality that it is not in $p$’s
connected component. Because it is in the ball’s interior, we have
$d_{\Omega}(c,r)<d_{\Omega}(c,q)$. By the triangle inequality, we obtain

$d_{\Omega}(p,r)\leavevmode\nobreak\ \leq\leavevmode\nobreak\
d_{\Omega}(p,c)+d_{\Omega}(c,r)\leavevmode\nobreak\ <\leavevmode\nobreak\
d_{\Omega}(p,c)+d_{\Omega}(c,q)\leavevmode\nobreak\ =\leavevmode\nobreak\
d_{\Omega}(p,q),$

but this contradicts the hypothesis that $p$ and $q$ were the closest sites in
different components.

Another useful property of the Delaunay triangulation, viewed as a graph, is
its relationship to other geometric graph structures. The _Hilbert minimum
spanning tree_ of a point set $P$, denoted $\text{MST}_{\Omega}(P)$, is the
spanning tree over $P$ where edge weights are Hilbert distances. The _Hilbert
relative neighborhood graph_ for a point set $P$, denoted
$\text{RNG}_{\Omega}(P)$ is a graph over the vertex set $P$, where two points
$p,q\in P$ are joined by an edge if
$d_{\Omega}(p,q)\leq\min(d_{\Omega}(p,r),d_{\Omega}(q,r))$, for all $r\in
P\setminus\\{p,q\\}$. Toussaint proved that in the Euclidean metric, the
minimum spanning tree is contained in the relative neighborhood graph, which
is contained in the Delaunay graph [25]. The following shows that this holds
in the Hilbert metric as well.

######  Lemma 2.6.

Given any discrete set of points $P$,
$\text{MST}_{\Omega}(P)\subseteq\text{RNG}_{\Omega}(P)\subseteq\operatorname{DT}_{\Omega}(P)$.

######  Proof 2.7.

The fact that the minimum spanning tree is contained in the relative
neighborhood graph holds in any metric. If two points $p$ and $q$ are not
joined by an edge in the RNG, then there is a point $r$ that is closer to both
than they are to each other. The edge $pq$ cannot be in the minimum spanning
tree, because it is the highest weight edge in the cycle formed by $p$, $q$,
and $r$.

To prove that the
$\text{RNG}_{\Omega}(P)\subseteq\operatorname{DT}_{\Omega}(P)$, consider an
edge $pq$ in the Hilbert RNG. Let $B(p,q)$ denote the Hilbert ball centered at
$p$ with radius $d_{\Omega}(p,q)$ and define $B(q,p)$ analogously. The
condition for two points $p,q\in P$ to be connected in
$\text{RNG}_{\Omega}(P)$ is that there is no other point of $P$ lying in the
interior of “lune” $B(p,q)\cap B(q,p)$. Consider the ball $B$ centered at the
point $c$ that is midway between $p$ and $q$, as introduced in the proof of
Lemma 2.4. For any $r\in B$, we showed that $d_{\Omega}(p,r)<d_{\Omega}(p,q)$.
Symmetrically, we have $d_{\Omega}(q,r)<d_{\Omega}(q,p)$, which implies that
$r$ lies in the lune. Since the lune must be empty, $B$ must also be empty.
This empty ball witnesses the fact that $pq$ is an edge of
$\operatorname{DT}_{\Omega}(P)$.

##  3 Hilbert Bisectors

In this section we present some utilities for processing Hilbert bisectors,
which will be used later in our algorithms. We start by recalling a
characterization given by Gezalyan and Mount for when a point lies on the
Hilbert bisector [13]. Recall that three lines are said to be _concurrent_ in
projective geometry if they intersect in a common point or are parallel.

######  Lemma 3.1.

Given a convex body $\Omega$ and two points
$p,q\in\operatorname{int}(\Omega)$, consider any other point $x\in\Omega$. Let
$p^{\prime}$ and $p^{\prime\prime}$ denote the endpoints of the chord
$\overline{px}$ so the points appear in the order $\langle
p^{\prime},p,x,p^{\prime\prime}\rangle$. Define $q^{\prime}$ and
$q^{\prime\prime}$ analogously for $\overline{qx}$.

  1. $(i)$

If $x\in\operatorname{int}(\Omega)$, then it lies on the $(p,q)$-bisector if
and only if the lines $\overleftrightarrow{pq}$,
$\overleftrightarrow{p^{\prime}q^{\prime}}$, and
$\overleftrightarrow{p^{\prime\prime}q^{\prime\prime}}$ are concurrent (see
Figure 5(a)).

  2. $(ii)$

If $x\in\partial\kern 1.0pt\Omega$ (implying that
$x=p^{\prime\prime}=q^{\prime\prime}$), then $x$ is the limit point of the
$(p,q)$-bisector if and only if there is a supporting line through $x$
concurrent with $\overleftrightarrow{pq}$ and
$\overleftrightarrow{p^{\prime}q^{\prime}}$ (see Figure 5(b)).

Figure 5: Properties of points on the Hilbert bisector.

The following utility lemmas arise as consequences of this characterization.
Both involve variants of binary search.

######  Lemma 3.2.

Given an $m$-sided convex polygon $\Omega$, two points
$p,q\in\operatorname{int}(\Omega)$, and any ray emanating from $p$, the point
of intersection between the ray and the $(p,q)$-bisector (if it exists) can be
computed in $O(\log^{2}m)$ time.

######  Proof 3.3.

Through binary search on the boundary of $\Omega$, in $O(\log m)$ time we can
determine the endpoints $p^{\prime}$ and $p^{\prime\prime}$ of the chord
passing through $p$ along the ray (see Figure 6(a)). Let $\ell$ denote this
chord, and let $V$ denote the points along the chord where the spokes of $q$
intersect $\ell$ (see Figure 6(b)).

Figure 6: Computing where a ray through a site hits the bisector.

While the points of $V$ will not be computed explicitly, the search algorithm
will simulate a binary search among these points. Observe that for any
subsegment $[x_{0},x_{1}]\subseteq\ell$, the two chords $\overline{qx_{0}}$
and $\overline{qx_{1}}$ subdivide the boundary of $\Omega$ into two polygonal
arcs so that every chord through $q$ that intersects this subsegment has its
endpoints in these two arcs. We can effectively simulate a binary search along
$[x_{0},x_{1}]$ by sampling the median vertex $v$ from whichever arc contains
the greater number of vertices, and then taking the point
$x=\overline{qv}\cap\ell$ (see Figure 6(c)). Each such probe will succeed in
eliminating at least one quarter of the remaining points of $V$.

Given a probe point $x$ along $\ell$, we first determine in $O(1)$ time
whether it is on the proper side of $p$ relative to the ray’s direction. If
not, we eliminate the portion lying behind $p$. If so, we compute the opposite
endpoint $v^{\prime}$ of $\overline{qv}$ in $O(\log m)$ time. Next, we compute
the distances $d_{\Omega}(p,x)$ and $d_{\Omega}(q,x)$ (which can be done in
$O(1)$ time since both chords are known). If the distance to $q$ is smaller,
we recurse on $[x_{0},x]$, and otherwise we recurse on $[x,x_{1}]$. In either
case, a constant fraction of remaining points have been eliminated.

The search eventually terminates either with an interval along $\ell$ between
two consecutive spokes of $q$, or it falls off the end of the ray. In the
former case, we can determine the exact location of the bisector within the
interval in $O(1)$ time, and report this answer. In the latter case, we report
that the ray does not intersect the bisector.

Each probe takes $O(\log m)$ time (to determine the opposite endpoint of each
ray through $q$). Since $V$ contains $O(m)$ points, and we eliminate a
constant fraction with each probe, the overall time is $O(\log^{2}m)$, as
desired.

######  Lemma 3.4.

Given an $m$-sided convex polygon $\Omega$ and any two points
$p,q\in\operatorname{int}(\Omega)$, the endpoints of the $(p,q)$-bisector on
$\partial\kern 1.0pt\Omega$ can be computed in $O(\log^{2}m)$ time.

######  Proof 3.5.

For the sake of illustration, let us assume that the line through $p$ and $q$,
denoted $\ell$, is horizontal (see Figure 7(a)). In $O(\log m)$ time we can
compute the two points, $z_{0}$ and $z_{1}$, where $\ell$ intersects
$\partial\kern 1.0pt\Omega$, with $z_{0}$ on the left side of $\Omega$ and
$z_{1}$ on the right. These points implicitly split the boundary of $\Omega$
into two convex chains, which we call the _upper_ and _lower chains_. We will
show how to compute the endpoint of the $(p,q)$-bisector that lies on the
upper chain, and the other case is symmetrical.

Figure 7: Computing the endpoints of the $(p,q)$-bisector.

It will be convenient to think of the portion of $\ell$ that lies outside of
$\Omega$ as being oriented from right to left, first starting at $z_{0}$, then
progressing to left to $-\infty$ and wrapping around to $+\infty$, and finally
progressing to $z_{1}$. Observe that as a point $x$ travels from $z_{0}$ to
$z_{1}$ along the upper chain of $\partial\kern 1.0pt\Omega$, the associated
supporting line intersects $\ell$ at a point that moves monotonically along
$\ell$ from $z_{0}$ to $z_{1}$ (wrapping around when it hits the topmost
vertex of $\Omega$).

By Lemma 3.1(ii), we seek a point $x^{*}$ on the upper chain that satisfies
the conditions of the lemma (see Figure 7(a)). Given any point $x$ on the
upper chain, consider the endpoints $p^{\prime}$ and $q^{\prime}$ of the
respective chords $\overline{xp}$ and $\overline{xq}$ (see Figure 7(b)). Let
$z$ denote the point where a supporting line at $x$ intersects $\ell$. Let
$z^{\prime}$ denote the point where the line
$\overleftrightarrow{p^{\prime}q^{\prime}}$ intersects $\ell$. In the limit,
when $x\rightarrow z_{0}$, we have $z=z_{0}$ and $z^{\prime}=z_{1}$, and in
the other limit, when $x\rightarrow z_{1}$, we have $z=z_{1}$ and
$z^{\prime}=z_{0}$. The key to the search process is that as $x$ moves
monotonically along the upper chain from $z_{0}$ to $z_{1}$, the point $z$
moves monotonically along $\ell$ from $z_{0}$ to $z_{1}$. On the other hand,
the points $p^{\prime}$ and $q^{\prime}$ move monotonically along the lower
chain from $z_{1}$ and $z_{0}$, and hence $z^{\prime}$ moves in the reverse
direction along $\ell$. It follows that there is a unique point $x$ on the
upper chain where $z$ and $z^{\prime}$ coincide along $\ell$. This will be the
desired point $x^{*}$.

We will show that $x^{*}$ can be found through a binary search along the upper
chain of $\Omega$, where each probe takes $O(\log m)$ time. We maintain two
spokes through $p$, which delimit the portion along the upper chain where the
desired point $x^{*}$ resides (see Figure 7(c)). This partitions the upper and
lower chain into two polygonal arcs containing the endpoints of the candidate
chords. With each probe, we sample the median vertex $x$ from the arc that
contains the larger number of vertices. This will allow us to eliminate at
least one quarter of the remaining vertices from consideration.

Depending on the nature of the spoke, the point $x$ may be a vertex of the
upper chain or it may lie in the interior of an edge. If it is in the interior
of an edge, the supporting line at $x$ is uniquely determined, implying that
$z$ is uniquely determined as well. In $O(\log m)$ time, we can determine the
points $p^{\prime}$ and $q^{\prime}$, and $z^{\prime}$ can be computed in
$O(1)$ additional time. Depending on where $z^{\prime}$ lies relative to $z$
in the order defined along $\ell$, we eliminate either the region before or
after $x$ on the upper chain, thus eliminating a constant fraction of the
points from further consideration. On the other hand, if $x$ is a vertex,
there are two supporting lines, which define an interval along $\ell$. Again,
we can determine $p^{\prime}$, $q^{\prime}$, and $z^{\prime}$ in $O(\log m)$
time. If $z^{\prime}$ lies within this interval, then $x$ is the desired final
answer. Otherwise, we eliminate either the region before or after $x$,
depending on which side of the interval $z^{\prime}$ lies.

Each probe takes $O(\log m)$ time (to determine the points $p^{\prime}$ and
$q^{\prime}$). Since there are $O(m)$ spokes through $p$, and we eliminate a
constant fraction with each probe, the overall time is $O(\log^{2}m)$, as
desired.

##  4 Hilbert Circumcircles

In the Euclidean plane, any triple of points that are not collinear lie on a
unique circle, that is, the boundary of an Euclidean ball. This is not true in
the Hilbert geometry, however. Since Delaunay triangulations are based on an
empty circumcircle condition, it will be important to characterize when a
triangle admits a Hilbert circumcircle and when it does not. In this section
we explore the conditions under which such a ball exists.

###  4.1 Balls at infinity

We begin by introducing the concept of a Hilbert ball centered at a point on
the boundary of $\Omega$. Let $x$ be any point on $\partial\kern 1.0pt\Omega$.
Given any point $p\in\operatorname{int}(\Omega)$, we are interested in
defining the notion of a Hilbert ball centered at $x$ whose boundary contains
$p$. If we think of the points on the boundary of $\Omega$ as being infinitely
far from any point in $\operatorname{int}(\Omega)$, this gives rise to the
notion of a ball at infinity.

Given $x\in\partial\kern 1.0pt\Omega$ and $p\in\operatorname{int}(\Omega)$,
let $u$ be any unit vector directed from $x$ into the interior of $\Omega$
(see Figure 8(a)). Given any sufficiently small positive $\delta$, let
$x_{\delta}=x+\delta u$ denote the point at distance $\delta$ from $x$ along
vector $u$, and let $B(x_{\delta},p)$ denote the Hilbert ball centered at
$x_{\delta}$ of radius $d_{\Omega}(x_{\delta},p)$. The following lemma shows
that as $\delta$ approaches $0$, this ball approaches a shape, which we call
the _Hilbert ball at infinity_ determined by $x$ and $u$ and passing through
$p$, denoted $B_{u}(x,p)$ (see Figure 8(b)).

Figure 8: Constructing the Hilbert ball $B_{u}(x,p)$ at infinity.

######  Lemma 4.1.

As $\delta$ approaches $0$, $B(x_{\delta},p)$ approaches a convex polygon
lying within $\Omega$ having $x$ on its boundary.

######  Proof 4.2.

We will provide a construction of the limiting shape by deriving the point at
which any chord through $x$ passes through the boundary of the shape. First,
consider the chord $\overline{px_{\delta}}$, letting $p^{\prime}$ and
$p^{\prime\prime}$ denote the two points of intersection with $\partial\kern
1.0pt\Omega$. Next, take any point $q^{\prime}\in\partial\kern 1.0pt\Omega$
different from $p^{\prime}$. Consider the chord
$\overline{q^{\prime}x_{\delta}}$, letting $q^{\prime\prime}$ denote the other
endpoint of the chord. Label the chord endpoints so that $p^{\prime}$ and
$q^{\prime}$ are on the same side as $p$ with respect to $x_{\delta}$ and
$p^{\prime\prime}$ and $q^{\prime\prime}$ are on the other (see Figure 8(b)
and (c)). Let $v$ denote the point where the lines
$\overleftrightarrow{p^{\prime}q^{\prime}}$ and
$\overleftrightarrow{p^{\prime\prime}q^{\prime\prime}}$ intersect. Let $q$
denote the point of intersection between the chord
$\overline{q^{\prime}q^{\prime\prime}}$ and the line
$\overleftrightarrow{pv}$. (It is easy to see that $q$ exists and is unique by
convexity.)

Observe that the point sequences $\langle
p^{\prime\prime},x_{\delta},p,p^{\prime}\rangle$ and $\langle
q^{\prime\prime},x_{\delta},q,q^{\prime}\rangle$ form a projectivity through
$v$, and hence their cross ratios are equal. Therefore, $p$ and $q$ are
equidistant from $x_{\delta}$, and therefore they both lie on the boundary of
a Hilbert ball centered at $x_{\delta}$.

All of these quantities depend on the value of $\delta$. Observe that in the
limit, as $\delta\rightarrow 0$, $x_{\delta}$ approaches $x$, and thus the
chords $x_{\delta}p$ and $x_{\delta}q^{\prime}$ approach $xp$ and
$xq^{\prime}$, respectively. Also observe that for all sufficiently small
$\delta$, the edges on which $p^{\prime\prime}$ and $q^{\prime\prime}$ lie do
not vary. Therefore, for these and all smaller values of $\delta$, the
triangles $\triangle xp^{\prime\prime}q^{\prime\prime}$ are similar, and hence
the slope of the line $\overleftrightarrow{p^{\prime\prime}q^{\prime\prime}}$
is fixed. (Fixing this slope is the reason that the limiting ball depends not
only on $x$ but on $u$ as well.) It follows that this line approaches the
supporting line $\ell$ through $x$ with this same slope. (Figure 8(c) shows
the case where $x$ is a vertex. When $x$ lies on the interior of an edge, the
choice of $u$ does not matter, and the supporting line is the unique
supporting line through this edge.)

Therefore, the point where lines $\overleftrightarrow{p^{\prime}q^{\prime}}$
and $\ell$ intersect approaches a limiting point $v$. It follows that $q$
approaches a point along the chord $xq^{\prime}$ in the limit. The set of all
such points taken over all choices of $q^{\prime}$ on the boundary of $\Omega$
yields the boundary of the desired shape. Since all the intermediate balls
$B(x_{\delta},p)$ are convex polygons lying within $\Omega$ containing
$x_{\delta}$, it follows that the limiting shape is also a convex polygon
containing $x$, which completes the proof.

A useful utility, which we will apply in Section 7, involves computing the
largest empty ball centered at any boundary point $x$ with respect to a point
set $P$.

######  Lemma 4.3.

Given a set of $n$ points $P$ in the interior of $\Omega$, and any point
$x\in\partial\kern 1.0pt\Omega$, in $O(n\log m)$ time, it is possible to
compute a point $p\in P$, such that there is a ball at infinity centered at
$x$ whose boundary passes through $p$, and which contains no points of $P$ in
its interior.

######  Proof 4.4.

Let $P=\\{p_{1},\ldots,p_{n}\\}$, and let $\ell$ be any supporting line
passing through $x$, and fix an arbitrary point $y\in\partial\kern
1.0pt\Omega$ other than $x$ (see Figure 9(a)). For each point $p\in P$, shoot
a ray from $x$ through $p$, letting $p^{\prime}$ denote the point where it
hits the boundary. Let $v$ denote the intersection of $\ell$ and the line
$\overleftrightarrow{p^{\prime}y}$, and let $q$ denote the intersection of
line $\overleftrightarrow{pv}$ with the chord $\overline{xy}$ (see Figure
9(b)).

Figure 9: Proof of Lemma 4.3.

It follows from Lemma 3.1 that both $p$ and $q$ lie on the boundary of a ball
at infinity centered at $x$. Applying this construction to each $p_{i}\in P$
yields a set of points $\\{q_{1},\ldots,q_{n}\\}$ along $\overline{xy}$, which
by Lemma 2.2 are ordered by inclusion (see Figure 9(c)). The desired point
$p_{i}\in P$ is the one whose corresponding point $q_{i}$ is closest to $x$
along this chord. The running time for each point $p_{i}$ is $O(\log m)$ to
determine where the ray hits the boundary, which yields an overall running
time of $O(n\log m)$.

Balls at infinity are not proper aspects of Hilbert geometry, but they will be
convenient for our purposes. Given two points
$p,q\in\operatorname{int}(\Omega)$, let $x$ denote the endpoint of the
$(p,q)$-bisector on $\partial\kern 1.0pt\Omega$, oriented so that $x$ lies to
the left of the directed line $\overrightarrow{pq}$. Let $u$ denote the
tangent vector of the bisector at $x$. Define $B(p{\,:\,}q)=B_{u}(x,p)$. Note
that this (improper) ball is both centered at and passes through $x$. In this
sense it circumscribes the triangle $\triangle pqx$. Define $B(q{\,:\,}p)$
analogously for the opposite endpoint of this bisector (see Figure 10(a) and
(b)).

We can now characterize the set of points $r$ that admit a Hilbert
circumcircle with respect to two given points $p$ and $q$. This
characterization is based on two regions, called the _overlap_ and _outer
regions_ (see Figure 10(c)).

######  Definition 4.5 (Overlap/Outer Regions).

Given two points $p,q\in\operatorname{int}(\Omega)$:

Overlap Region:

denoted $Z(p,q)$, is $B(p{\,:\,}q)\cap B(q{\,:\,}p)$.

Outer Region:

denoted $W(p,q)$, is $\Omega\setminus(B(p{\,:\,}q)\cup B(q{\,:\,}p))$.

Figure 10: The overlap region $Z(p,q)$ and the outer region $W(p,q)$.

######  Lemma 4.6.

A triangle $\triangle pqr\subseteq\operatorname{int}(\Omega)$ admits a Hilbert
circumcircle if and only if $r\notin Z(p,q)\cup W(p,q)$.

######  Proof 4.7.

Let $x,x^{\prime}\in\partial\kern 1.0pt\Omega$ denote the endpoints of the
$(p,q)$-bisector. Since the bisector is a curve, we can parameterize points
along this curve by $s:[0,1]\rightarrow\Omega$ so that $s(0)=x$ and
$s(1)=x^{\prime}$. For each $t\in(0,1)$ there exists a unique Hilbert ball
centered at $s(t)$ and passing through $p$ and $q$, denoted
$B_{t}(p{\,:\,}q)$. Lemma 2.2 implies that any point $r$ other than $p$ and
$q$ that lies on the boundary of one of these circles, can lie on at most one
such circle. That is, the set $B_{t}(p{\,:\,}q)$ defines a “pencil” of Hilbert
circles (see Figure 10(d)) ranging from $B(p{\,:\,}q)=B_{0}(p{\,:\,}q)$ to
$B(q{\,:\,}p)=B_{1}(p{\,:\,}q)$. Conversely, any Hilbert circumcircle passing
through $p$ and $q$ must be generated by a point on the $(p,q)$-bisector, and
therefore every point $r$ that completes a Hilbert circumcircle lies on one of
the circles in this pencil. It is evident therefore that $r$ cannot lie within
either $Z(p,q)$ nor $W(p,q)$, for otherwise the boundary of the resulting ball
would cross the boundary of either $B(p{\,:\,}q)$ or $B(q{\,:\,}p)$, thus
violating Lemma 2.2 at this crossing point.

As shown in Figure 1, the Delaunay triangulation need not cover the convex
hull of the set of sites. We refer to the region that is covered as the
_Hilbert hull_ of the sites. Later, we will present an algorithm for computing
the Hilbert hull. The following lemma will be helpful. It establishes a
nesting property for the overlap regions.

######  Lemma 4.8.

If $r\in Z(p,q)$ then $Z(p,r)\subset Z(p,q)$.

######  Proof 4.9.

Given two points $p$ and $q$, let $r$ be a point in $Z(p,q)$. Consider that
the bisector between $p$ and $r$ and the bisector between $q$ and $r$ split
$\Omega$ into three pieces that can only connect at infinity. Let $u$ be a
point in $Z(p,r)$. Then this would necessarily require that the bisectors
between $p$ and $u$ and $u$ and $r$ would split $\Omega$ again such that they
only can connect at infinity. Note that every point on the $u,r$ bisector is
necessarily closer to $u$ than to $r$ and therefore $q$. Hence $u$ must split
space between $p$ and $q$ as well. Since $u$ was defined to be in one of the
balls at infinity it must therefore by in $Z(p,q)$ as it cannot lie in
$W(p,q)$.

##  5 Computing Circumcircles

A fundamental primitive in the Euclidean Delaunay triangulation algorithm is
the so-called _in-circle test_ [14]. Given a triangle $\triangle pqr$ and a
fourth site $s$, the test determines whether $s$ lies within the
circumscribing Hilbert ball for the triangle (if such a ball exits). In this
section, we present an algorithm which given any three sites either computes
the Hilbert circumcircle for these sites, denoted $B(p{\,:\,}q{\,:\,}r)$, or
reports that no circumcircle exists. The in-circle test reduces to checking
whether $d_{\Omega}(s,c)<\rho$, which can be done in $O(\log m)$ time as
observed in Section 2.

######  Lemma 5.1.

Given a convex $m$-sided polygon $\Omega$ and triangle $\triangle
pqr\subset\operatorname{int}(\Omega)$, in $O(\log^{3}m)$ time it is possible
to compute $B(p{\,:\,}q{\,:\,}r)$ or to report that no ball exists.

The remainder of the section is devoted to the proof. The circumscribing ball
exists if and only if there is a point equidistant to all three, implying that
$(p,q)$\- and $(p,r)$-bisectors intersect at some point
$c\in\operatorname{int}(\Omega)$. The following technical lemma shows that
bisectors intersect crosswise.

######  Lemma 5.2.

Given three non-collinear points $p,q,r\in\operatorname{int}(\Omega)$, the
$(p,q)$\- and $(p,r)$-bisectors have endpoints lying on $\partial\kern
1.0pt\Omega$. If they intersect within $\operatorname{int}(\Omega)$, they
intersect transversely in a single point (see Figure 11).

######  Proof 5.3.

As shown in [13], star-shapedness of Voronoi cells implies that bisectors are
not closed loops. Since each bisector is a simple curve, a bisector has two
endpoints on the boundary of $\Omega$. If the $(p,q)$\- and $(p,r)$-bisectors
intersect, the intersection point is equidistant from all three, and hence is
the center of a circumcircle. By the above lemma implies that there cannot be
more than one intersection. The intersection must be transverse. If not, the
regions on opposite sides of the intersection point would be closer to $p$
than to $q$ or $r$, and this would violate the property the Voronoi cells are
star-shaped.

Figure 11: Circumcircles and bisector intersection.

It follows that the $(p,q)$\- and $(p,r)$-bisectors subdivide the interior of
$\Omega$ into either three or four regions (three if they do not intersect,
and four if they do). We can determine which is this case by invoking Lemma
3.4 to compute the endpoints of these bisectors in $O(\log^{2}m)$ time. If
they alternate between $(p,q)$ and $(p,r)$ along the boundary of $\Omega$,
then the bisectors intersect, and otherwise they do not. In the latter case,
there is no circumcircle for $p$, $q$, and $r$, and hence, no possibility of
violating the circumcircle condition. (Note that this effectively provides an
$O(\log^{2}m)$ test for Lemma 4.6.) Henceforth, we concentrate on the former
case.

Let us assume, without loss of generality, that $\triangle pqr$ is oriented
counterclockwise. Let $v_{q}$ denote the endpoint of the $(p,q)$-bisector that
lies to the right of the oriented line $\overrightarrow{pq}$. Let $v_{r}$
denote the endpoint of the $(p,r)$-bisector that lies to the left of the
oriented line $\overrightarrow{pr}$ (see Figure 12(a)).

Figure 12: Computing the center of $\triangle{pqr}$.

Because Voronoi cells are star-shaped, it follows that the vector from $p$ to
the desired circumcenter $c$ lies in the counterclockwise angular interval
from $\overrightarrow{pv_{q}}$ to $\overrightarrow{pv_{r}}$ (see Figure
12(b)). The key to the search is the following observation. In the angular
region from $\overrightarrow{pv_{q}}$ to $\overrightarrow{pc}$, any ray shot
from $p$ intersects the $(p,q)$-bisector before hitting the $(p,r)$-bisector
(if it hits the $(p,r)$-bisector at all). On the other hand, in the angular
region from $\overrightarrow{pc}$ to $\overrightarrow{pv_{r}}$, any ray shot
from $p$ intersects the $(p,r)$-bisector before hitting the $(p,q)$-bisector
(if it hits the $(p,q)$-bisector at all). We say that the former type of ray
is _early_ and the latter type is _late_.

Let $v_{-}$ and $v_{+}$ denote the points on $\partial\kern 1.0pt\Omega$ that
bound the current search interval about $p$ (see Figure 12(c)). We will
maintain the invariant that the ray $\overrightarrow{pv_{-}}$ is early and
$\overrightarrow{pv_{+}}$ is late. Initially, $v_{-}=v_{q}$ and $v_{+}=v_{r}$.
Let $v^{\prime}_{-}$ and $v^{\prime}_{+}$ denote the opposite endpoints of the
chords $\overline{pv_{-}}$ and $\overline{pv_{+}}$. Consider the portion of
the boundary of $\Omega$ that lies in the counterclockwise from $v_{-}$ and
$v_{+}$. If the angle $\angle v_{-}pv_{+}$ is smaller than $\pi$, also
consider the portion of the boundary of $\Omega$ that lies in the
counterclockwise from $v^{\prime}_{-}$ and $v^{\prime}_{+}$. With each probe,
we sample the median vertex $v$ from whichever of these two boundary portions
that contains the larger number of vertices. If $v$ comes from the interval
$[v_{-},v_{+}]$, we probe along the ray $\overrightarrow{pv}$, and if it comes
from the complementary interval, $[v^{\prime}_{-},v^{\prime}_{+}]$, we shoot
the ray from $p$ in the opposite direction from $v$. We then apply Lemma 3.2
twice to determine in $O(\log^{2}m)$ time where this ray hits the $(p,q)$\-
and $(p,r)$-bisectors (if at all). Based on the results, we classify this ray
as being early or late and recurse on the appropriate angular subinterval.
When the search terminates, we have determined a pair of consecutive spokes
about $p$ that contain $c$.

Because each probe eliminates at least half of the vertices from the larger of
the two boundary portions, it follows that at $O(\log m)$ probes, we have
located $c$ to within a single pair of consecutive spokes around $p$. Since
each probe takes $O(\log^{2}m)$ time, the entire search takes $O(\log^{3}m)$
time.

We repeat this process again for $r$ and $q$. The result is three double
wedges defined by consecutive spokes, one about each site. It follows from our
earlier remarks from Section 2.2 that within the intersection of these
regions, the bisectors are simple conics. We can compute these conics in
$O(1)$ times [8] and determine their intersection point, thus yielding the
desired point $c$. The radius $\rho$ of the ball can also be computed in
$O(1)$ time.

##  6 Building the Triangulation

In this section we present our main result, a randomized incremental algorithm
for constructing the Delaunay triangulation $\operatorname{DT}(P)$ for a set
of $n$ sites $P$ in the interior of an $m$-sided convex polygon $\Omega$. Our
algorithm is loosely based on a well-known randomized incremental algorithm
for the Euclidean case [14, 10].

###  6.1 Orienting and Augmenting the Triangulation

In this section we introduce some representational conventions for the sake of
our algorithm. First, we will orient the elements of the triangulation. Given
$p,q\in\operatorname{int}(\Omega)$ define the _endpoint_ of the
$(p,q)$-bisector to be the endpoint that lies to the left of the directed line
$\overrightarrow{pq}$. (It follows from the star-shapedness of Voronoi cells
that the bisector endpoints intersect $\partial\kern 1.0pt\Omega$ on opposite
sides of this line.) The opposite endpoint will be referred to as the
$(q,p)$-bisector endpoint. When referring to a triangle $\triangle pqr$, we
will assume that the vertices are given in counterclockwise order. Also, edges
in the triangulation are assumed to be directed, so we can unambiguously
reference the left and right sides of a directed edge $pq$.

As mentioned earlier, the triangulation need not cover the convex hull of the
set of sites (see Figure 13(a)). For the sake of construction, it will be
convenient to augment the triangulation with additional elements, so that all
of $\Omega$ is covered. First, we add elements to include the endpoints of the
Voronoi bisectors on $\partial\kern 1.0pt\Omega$. For each edge $pq$ such that
the external face of the triangulation lies to its left, the endpoint of the
$(p,q)$-bisector intersects the boundary of $\Omega$. Letting $x$ denote this
boundary point, we add a new triangle $\triangle pqx$, called a _tooth_. (See
the green shaded triangles in Figure 13(b).)

Figure 13: The augmented triangulation.

Finally, the portion of $\Omega$ that lies outside of all the standard
triangles and teeth consists of a collection of regions, called _gaps_ , each
of which involves a single site, the sides of two teeth, and a convex
polygonal chain along $\partial\kern 1.0pt\Omega$. While these shapes are not
triangles, they are defined by three points, and we will abuse the notation
$\triangle pxy$ to refer to the gap defined by the site $p$ and boundary
points $x$ and $y$ (see Figure 13(b)).

Even when sites are in general position, multiple teeth can share the same
boundary vertex. For example, in Figure 13(b) three teeth meet at the same
boundary vertex $v$. In our augmented representation, it will be convenient to
treat these as three separate teeth, meeting at three distinct (co-located)
vertices, separated by two degenerate (zero-width) gaps. This allows us to
conceptualize the region outside of the standard triangulation as consisting
of an alternating sequence of teeth and gaps. The following lemma summarizes a
few useful facts about the augmented representation.

######  Lemma 6.1.

Given a set of $n$ point sites in the interior of a convex polygon $\Omega$:

  * •

The augmented Delaunay triangulation has complexity $O(n)$.

  * •

The region covered by standard triangles is connected.

  * •

Each standard triangle and each tooth satisfies the empty circumcircle
property.

  * •

The region outside the standard triangles consists of an alternating sequence
of teeth and gaps.

  * •

For any $p\in\Omega$, membership in any given triangle, tooth, or gap can be
determined in $O(1)$ time.

######  Proof 6.2.

  * •

The standard triangulation is clearly $O(n)$ by standard results on
triangulations (based on Euler’s formula). The number of teeth is clearly
$O(n)$, since each is generated by a standard edge.

  * •

The connectivity of the triangulation was proved in Lemma 2.4.

  * •

Through the addition of degenerate gaps, the teeth and gaps alternate.
Membership in a standard triangle or tooth can be computed in $O(1)$ time.

  * •

Since each gap is bounded by two edges and the boundary of $\Omega$,
membership for any $p\in\Omega$ can be determined in $O(1)$ time.

###  6.2 Local and Global Delaunay

Given a set $P$ of point sites, we say that an augmented triangulation
$\mathcal{T}(P)$ is _globally Delaunay_ (or simply _Delaunay_) if for each
standard triangle and each tooth $\triangle pqr$, the circumscribing ball,
denoted $B(p{\,:\,}q{\,:\,}r)$, exists and has no site within its interior.
The incremental Euclidean Delaunay triangulation algorithm employs a local
condition, which only checks this for neighboring triangles [14, 10]. Given a
directed edge $pq$ of the triangulation that is incident to two triangles or
to a triangle and tooth, let $a$ and $b$ be the vertices of the incident
triangles lying to the left and right of the $pq$, respectively. We say that
$\mathcal{T}(P)$ is _locally Delaunay_ if
$a\notin\operatorname{int}(B(q{\,:\,}p{\,:\,}b))$ and
$b\notin\operatorname{int}(B(p{\,:\,}q{\,:\,}a))$. The following lemma shows
that it suffices to certify the Delaunay properties locally. The proof follows
the same structure as the Euclidean case, but additional care is needed due to
the existence of teeth and gaps.

######  Lemma 6.3.

Given a set $P$ of point sites, an augmented triangulation $\mathcal{T}(P)$ is
globally Delaunay if and only if it is locally Delaunay.

######  Proof 6.4.

Our proof follows the same structure as that given in [10] for the Euclidean
case. Suppose to the contrary that $\mathcal{T}(P)$ is locally but not
globally Delaunay. That is, there exist sites $a,b,c,d\in P$ such that
$\mathcal{T}(P)$ has a triangle or tooth $\triangle abc$ such that the
circumscribing ball $B(a{\,:\,}b{\,:\,}c)$ contains $d$ in its interior (see
Figure 14(a)). Let $ab$ be the edge such that $c$ and $d$ lie on opposite
sides of $\overline{ab}$. Among all triangle-site pairs $(\triangle abc,d)$,
let this be the one that maximizes $\angle adb$.

Figure 14: Proof of Lemma 6.3.

To complete the proof, we consider two cases. First, if $\triangle abc$ is a
standard triangle or if it is a tooth and both $a$ and $b$ are sites, consider
the adjacent triangle $\triangle bae$ along edge $ab$. Since $\mathcal{T}(P)$
is locally Delaunay, $e$ does not lie within the interior of
$B(a{\,:\,}b{\,:\,}c)$. By Lemma 2.2, the circumscribing ball
$B(b{\,:\,}a{\,:\,}e)$ completely contains the portion of
$B(a{\,:\,}b{\,:\,}c)$ that lies on the same side of the chord $\overline{ab}$
as $e$. (If it failed to do so, there would be an additional point of
intersection of the boundaries of these balls, other than $a$ and $b$,
violating the lemma.) This implies that $d\in B(b{\,:\,}a{\,:\,}e)$. Now, let
$be$ denote the edge of $\triangle bae$ such that $a$ and $d$ lie on the
opposite sides of the chord $\overline{be}$. It follows that $\angle
edb>\angle adb$, contradicting the definition of the pair $(\triangle abc,d)$.

In the second case, $\triangle abc$ is a tooth, and either $a$ or $b$ is a
boundary vertex, implying that there is a gap on the other side of $ab$. In
this case, relabel the vertices so that $a$ is a site and $b$ is on the
boundary. Let $\triangle b^{\prime}ae$ denote the tooth adjacent to the
opposite side of this gap (see Figure 14(b)). Now, apply the above, using
$\triangle b^{\prime}ae$ in place of $\triangle bae$. The key property still
holds, namely that circumscribing ball $B(b^{\prime}{\,:\,}a{\,:\,}e)$
completely contains the portion of $B(a{\,:\,}b{\,:\,}c)$ that lies on the
same side of the chord $\overline{ab}$ as $e$.

###  6.3 Incremental Construction

The algorithm operates by first randomly permuting the sites. It begins by
inserting two arbitrary sites $a$ and $b$ and generating the resulting
augmented triangulation. In $O(\log^{2}m)$ time, we can compute the endpoints
of the $(a,b)$-bisector. We add the edge $ab$ and create two teeth by
connecting $a$ and $b$ to the bisector endpoints (see Figure 15(a)). We then
insert the remaining points one by one, updating the triangulation
incrementally after each insertion (see Algorithm 1). Later, we will discuss
how to determine which element each new site lies in, but for now, let us
focus on how the triangulation is updated.

Algorithm 1 Constructs the Hilbert Delaunay triangulation of a point set $P$

procedure Delaunay($P$) $\triangleright$ Build the Delaunay triangulation of
point set $P$

$\mathcal{T}\leftarrow\text{empty triangulation}$

Randomly permute $P$

$a,b\leftarrow\text{any two points of $P$}$ $\triangleright$ Initialize with
sites $a$ and $b$

$x,y\leftarrow\text{endpoints of the $(a,b)$-bisector}$ $\triangleright$ See
Figure 15(a)

Add edges $ab$, $ax$, $ay$, $bx$, $by$ to $\mathcal{T}$

for all $p\in P\setminus\\{a,b\\}$ do $\triangleright$ Add all remaining sites

Insert($p,\mathcal{T}$)

end for

return $\mathcal{T}$

end procedure

Figure 15: (a) Initialization and insertion into a (b) standard triangle, (c)
tooth, (d) gap.

When we insert a point, there are three possible cases, depending on the type
of element that contains the point, a standard triangle (three sites), a tooth
(two sites), or a gap (one site). The case for standard triangles is the same
as for Euclidean Delaunay triangulations [14], involving connecting the new
site to the triangle’s vertices (see Figure 15(b)). Insertion into a tooth
involves removing the boundary vertex, connecting the new site to the two
existing site vertices, and creating two new teeth based on the bisectors to
these sites (see Figure 15(c)). Finally, insertion into a gap involves
connecting the new site to the existing site vertex, and creating two new
teeth based on the bisectors to this site (see Figure 15(d)).

Algorithm 2 Site insertion

procedure Insert($p,\mathcal{T}$) $\triangleright$ Insert a new site $p$ into
triangulation $\mathcal{T}$

$\triangle abc\leftarrow\text{the triangle of $\mathcal{T}$ containing $p$}$

if $\triangle abc$ is a standard triangle then $\triangleright$ See Figure
15(b)

Add edges $ap$, $bp$, $cp$ to $\mathcal{T}$

FlipEdge($ab,p,\mathcal{T}$); FlipEdge($bc,p,\mathcal{T}$);
FlipEdge($ca,p,\mathcal{T}$)

else if $\triangle abc$ is a tooth then $\triangleright$ See Figure 15(c)

Let $a$ and $b$ be sites, and $c$ be on boundary

$x,y\leftarrow\text{endpoints of the $(a,p)$- and $(p,b)$-bisectors,
respectively}$

Remove $c$ and edges $ca$ and $cb$ from $\mathcal{T}$

Add edges $ap$, $bp$, $ax$, $px$, $by$, $py$ to $\mathcal{T}$

FlipEdge($ab,p,\mathcal{T}$);

FixTooth($\triangle apx,p,\mathcal{T}$); FixTooth($\triangle
pby,p,\mathcal{T}$)

else$\triangleright$ $\triangle abc$ is a gap; see Figure 15(d)

Let $a$ be the site, and let $b$ and $c$ be on the boundary

$x,y\leftarrow\text{endpoints of $(a,p)$- and $(p,a)$-bisectors,
respectively}$

Add edges $ap$, $ax$, $px$, $ay$, $py$ to $\mathcal{T}$

FixTooth($\triangle apx,p,\mathcal{T}$); FixTooth($\triangle
pay,p,\mathcal{T}$)

end if

end procedure

On return from Insert, $\mathcal{T}$ is a topologically valid augmented
triangulation, but it may fail to satisfy the Delaunay empty circumcircle
conditions, and may even fail to be geometrically valid. The procedures
FlipEdge and FixTooth repair any potential violations of the local Delaunay
conditions.

The procedure FlipEdge is applied to all newly generated standard triangles
$\triangle pab$. It is given the directed edge $ab$ of the triangle, such that
$p$ lies to the left of this edge. It accesses the triangle $\triangle abc$
lying to the edge’s right. This may be a standard triangle or a tooth. In
either case, it has an associated circumcircle, which we assume has already
been computed.111In the Euclidean case, the in-circle test may be run from
either $\triangle abc$ or $\triangle pab$, but this is not true for the
Hilbert geometry. In light of Lemma 4.6, we do not know whether the $\triangle
pab$ has a circumcircle, so we run the test from $\triangle abc$. We test
whether $p$ encroaches on this circumcircle, and if so, we remove the edge
$ab$. If $\triangle{abc}$ is a standard triangle, then we complete the edge-
flip by adding edge $pc$, and then continue by checking the edges $ac$ and
$cb$ (see Figure 16(a)). Otherwise, we remove vertex $c$, and compute the
endpoints $x$ and $y$ of the $(p,a)$\- and $(b,p)$-bisectors, respectively. We
create two new teeth, $\triangle pax$ and $\triangle bpy$ (see Figure 16(b)).
This creates a new (possibly degenerate) gap $\triangle pxy$. Later, we will
show (see Lemma 6.5) that no further updates are needed. The algorithm is
presented in Algorithm 3.

Figure 16: (a) Standard-triangle edge flip, (b) tooth edge flip, and (c)
fixing a tooth. Algorithm 3 In-circle test and edge flip.

procedure FlipEdge($ab,p,\mathcal{T}$) $\triangleright$ In-circle test. New
site $p$ lies to the left of edge $ab$.

$c\leftarrow\text{vertex of triangle to right of $ab$ in $\mathcal{T}$}$

if $p\in B(a{\,:\,}c{\,:\,}b)$ then $\triangleright$ $p$ fails the in-circle
test for $\triangle acb$

Remove edge $ab$ from $\mathcal{T}$

if $\triangle acb$ is a standard triangle then $\triangleright$ See Figure
16(a)

Add edge $pc$ to $\mathcal{T}$

FlipEdge($ac,p,\mathcal{T}$); FlipEdge($cb,p,\mathcal{T}$) $\triangleright$
Create $\triangle pac$ and $\triangle bpc$

else$\triangleright$ $\triangle acb$ is a tooth. See Figure 16(b)

$x,y\leftarrow\text{endpoints of $(p,a)$- and $(b,p)$-bisectors,
respectively}$

Remove $c$ and edges $cb$ and $ca$ from $\mathcal{T}$

Add edges $ax$, $px$, $by$, $py$ to $\mathcal{T}$ $\triangleright$ Create
teeth $\triangle pax$ and $\triangle bpy$

end if

end if

end procedure

The “else” clause creates two teeth $\triangle pax$ and $\triangle bpy$. The
following lemma shows that there is no need to apply FixTooth to them.

######  Lemma 6.5.

On return from FlipEdge the teeth $\triangle pax$ and $\triangle bpy$
generated in the “else” clause are both valid.

######  Proof 6.6.

Consider the two teeth lying immediately clockwise and counterclockwise from
$c$ along the boundary of $\Omega$. Let $x^{\prime}$ and $y^{\prime}$ denote
the respective boundary vertices of these two teeth. We will show that $x$
lies within the portion of the boundary of $\Omega$ running clockwise from $c$
to $x^{\prime}$. It will follow from a symmetrical argument that $y$ lies in
the counterclockwise portion of the boundary from $c$ to $y^{\prime}$.

By the preconditions of FlipEdge, we know that $p$ lies to the left of edge
$ab$, and since we have executed the “else” clause, we know that $p$
encroaches on the circumscribing ball $B(a{\,:\,}c{\,:\,}b)$. Let $\triangle
adx^{\prime}$ denote the tooth lying clockwise from $c$ prior to the insertion
of $p$. If $x$ lies between the $c$ and $x^{\prime}$, then the circumscribing
balls for these triangles can be arranged as shown in Figure 17(a).

Figure 17: Proof of Lemma 6.5.

Suppose, to the contrary, that $x$ lies on the other side of $x^{\prime}$ (see
Figure 17(b)), observe that the circumcircle $B(p{\,:\,}a{\,:\,}x)$ must cross
through the circumcircle $B(a{\,:\,}d{\,:\,}x^{\prime})$. It follows that, in
addition to $a$ and $d$, these two Hilbert circumcircles intersect in at least
one other point $z$. However, this means that there are two different Hilbert
circumcircles passing through the three points $a$, $d$, and $z$, which cannot
happen according to Lemma 2.2.

Finally, we present FixTooth. To understand the issue involved, consider
Figure 15(d). In the figure, the newly created teeth $\triangle pax$ and
$\triangle apy$ both lie entirely inside the existing gap. But, this need not
generally be the case, and the newly created boundary vertex may lie outside
the current gap, resulting in two or more teeth that overlap each other (see,
e.g., Figure 16(c)). The procedure FixTooth is given a tooth $\triangle abx$,
where $a$ and $b$ are sites, and $x$ is on $\Omega$’s boundary. Whenever it is
invoked the new site $p$ is either $a$ or $b$. Recalling our assumption that
teeth and gaps alternate around the boundary of $\Omega$, we check the teeth
that are expected to be adjacent to the current tooth on the clockwise and
counterclockwise sides. (Only the clockwise case is presented in the
algorithm, but the other case is symmetrical, with $a$ and $b$ swapped.)

Let $\triangle bcy$ denote the tooth immediately clockwise around the
boundary. We test whether these triangles overlap (see Figure 16(c)). If so,
we know that the $(a,b)$-bisector (which ends at $x$) and the $(b,c)$-bisector
(which ends at $y$) must intersect. This intersection point is the center of a
Hilbert ball circumscribing $\triangle abc$. We replace the two teeth
$\triangle abx$ and $\triangle bcy$ with the standard triangle $\triangle abc$
and the tooth $\triangle acz$, where $z$ is the endpoint of the
$(a,c)$-bisector. If $p=a$, then we invoke FlipEdge on the opposite edge $bc$
from $p$, and we invoke FixTooth on the newly created tooth. When the
algorithm terminates, all the newly generated elements have been locally
validated. The algorithm is presented in Algorithm 4.

Algorithm 4 Check for and fix overlapping teeth.

procedure FixTooth($\triangle abx,p,\mathcal{T}$) $\triangleright$ Fix tooth
$\triangle abx$ where $p=a$ or $p=b$

Let $\triangle bcy$ be the tooth clockwise adjacent to $\triangle abx$

if $\triangle abx$ overlaps $\triangle bcy$ then $\triangleright$ See Figure
16(c)

$z\leftarrow\text{endpoint of the $(a,c)$-bisector}$

Remove $x$ and $y$ and edges $ax$, $bx$, $by$, and $cy$ from $\mathcal{T}$

Add edges $ac$, $az$, and $cz$ to $\mathcal{T}$ $\triangleright$ Create
triangle $\triangle abc$ and tooth $\triangle acz$

if $p=a$ then

FlipEdge($bc,p,\mathcal{T}$) $\triangleright$ Check edge flip with triangle
opposite $bc$

FixTooth($\triangle acz,p,\mathcal{T}$) $\triangleright$ Check for further
overlaps

end if

else

Repeat the above for the triangle counterclockwise from $\triangle abx$
swapping $a\leftrightarrow b$

end if

end procedure

Based on our earlier remarks, it follows that our algorithm correctly inserts
a site into the augmented triangulation.

######  Lemma 6.7.

Given an augmented triangulation $\mathcal{T}(P)$ for a set of sites $P$, the
procedure Insert correctly inserts a new site $p$, resulting in the augmented
Delaunay triangulation for of $P\cup\\{p\\}$.

The final issue is the algorithm’s expected running time. Our analysis follows
directly from the analysis of the randomized incremental algorithm for the
Euclidean Delaunay triangulation. We can determine which triangle contains
each newly inserted point in amortized $O(\log n)$ expected time either by
building a history-based point location data structure (as with Guibas, Knuth,
and Sharir [14]) or by bucketing sites (as with de Berg et al. [10]). The
analysis in the Euclidean case is based on a small number of key facts, which
apply in our context as well. First, sites are inserted in random order.
Second, the conflict set for any triangle or tooth consists of the points
lying in the triangle’s circumcircle. Both of these clearly hold in the
Hilbert setting. Third, the number of structural updates induced by the
insertion of site $p$ is proportional to the degree of $p$ following the
insertion. This holds for our algorithm, because each modification to the
triangulation induced by $p$’s insertion results in a new edge being added to
$p$ (and these edges are not deleted until future insertions). Fourth, the
structure is invariant to the insertion order. Finally, the triangulation
graph is planar, and hence it has constant average degree. The principal
difference is that the in-circle test involves computing a Hilbert
circumcircle, which by Lemma 5.1 can be performed in $O(\log^{3}m)$ time.
These additional factors of $O(\log n)$ and $O(\log^{3}m)$ are performed a
constant number of times in expectation, for each of the $n$ insertions. This
implies our main result.

######  Theorem 6.8.

Given a set of $n$ points in the Hilbert geometry defined by a convex $m$-gon
$\Omega$, it is possible to construct the augmented Delaunay triangulation in
randomized expected time $O(n(\log^{3}m+\log n))$.

##  7 The Hilbert Hull

As observed earlier, the triangles of the Delaunay triangulation of $P$ do not
necessarily cover the convex hull of $P$. The region covered by these
triangles is called the _Hilbert hull_ (the blue region of Figure 13). In this
section, we present a simple algorithm for computing this hull for a set of
$n$ points in the Hilbert distance defined by a convex $m$-gon $\Omega$. Our
approach is roughly based on the Jarvis march algorithm for computing convex
hulls.

Henceforth, when we talk about a site being _on the Hilbert hull_ , we mean
that it is a vertex of the hull’s boundary. These are the sites whose Voronoi
cells extend to the boundary of $\Omega$. Our algorithm will output a
couterclockwise enumeration of the sites that appear on the hull. By Lemma
2.4, the Hilbert hull is simply connected, but its boundary may collapse
around a vertex or an edge. If so, a counterclockwise traversal of the
boundary may encounter a given site multiple times. (For example, Figure 18
shows an example of four sites whose Hilbert hull degenerates to a tree. The
final hull output would be $\langle a,b,c,b,d,b,a\rangle$.)

Figure 18: A set of four sites whose Hilbert hull degenerates to a tree.

Let $p$ be a point on the Hilbert hull of $P$. Since its Voronoi cell extends
to the boundary, there exists a point $x\in\partial\kern 1.0pt\Omega$ such
that the Hilbert ball at infinity, $B(x,p)$, centered at $x$ and passing
through $p$ contains no other points of $P$ (see Figure 19(a)). The point $x$
serves as a witness to the fact that $p$ is on the Hilbert hull.

Figure 19: Computing the Hilbert hull.

The next point of $P$ counterclockwise on the Hilbert hull is the site
$p^{\prime}$ associated with the neighboring Voronoi cell in counterclockwise
order around the boundary of $\Omega$ (see Figure 19(b)). Let $x^{\prime}$
denote the endpoint of the $(p^{\prime},p)$-bisector. In general, multiple
Voronoi cells may meet at the same point $x^{\prime}$, and if so, we can use
the angle formed between $x$, $x^{\prime}$ and $p^{\prime}$ to break ties (see
Figure 19(c)).

Based on this observation, our adaptation of Jarvis march involves ordering
the sites of $p^{\prime}\in P\setminus\\{p\\}$ based on the a lexicographical
ordering the bisector endpoints. The points are ordered first by the endpoint
of the $(p^{\prime},p)$-bisector, in counterclockwise order starting at $x$,
and ties are broken by considering the clockwise angle $\angle
xx^{\prime}p^{\prime}$ (ranging over the interval $[0,2\pi)$). Call this the
_ordering induced by $p$ and $x$_.

######  Lemma 7.1.

Given a set $P$ of $n\geq 2$ sites in the Hilbert geometry defined by a convex
body $\Omega$, a site $p\in P$ that lies on the boundary of $P$’s Hilbert
hull, and a point $x\in\partial\kern 1.0pt\Omega$ such that $B(x,p)$ is empty,
the next site counterclockwise on the boundary of the Hilbert hull is the site
in $P\setminus\\{p\\}$ that is smallest in the order induced by $p$ and $x$.

######  Proof 7.2.

As observed earlier, the fact that the ball $B(x,p)$ is empty implies that $x$
lies in the (closure) of the Voronoi cell of $p$. Let $x^{\prime}$ (possibly
equal to $x$) be the extreme counterclockwise point of
$\operatorname{Vor}(p)\cap\partial\kern 1.0pt\Omega$, and let $p^{\prime}$
denote the site such that $\operatorname{Vor}(p^{\prime})$ and
$\operatorname{Vor}(p)$ share a common bisector whose endpoint is $x$. The
site $p^{\prime}$ is next on the Hilbert hull, and thus, it suffices to show
that $p^{\prime}$ is the smallest site in the order induced by $p$ and $x$.

For each site $p^{\prime\prime}\in P\setminus\\{p\\}$, let
$V(p,p^{\prime\prime})$ denote the set of points of $\operatorname{int}\Omega$
that are closer to $p$ than to $p^{\prime\prime}$ in terms of Hilbert
distance. This is just the set of points that lie on $p$’s side of the
$(p^{\prime\prime},p)$-bisector. By definition,
$\operatorname{Vor}(p)=\bigcap_{p^{\prime\prime}\neq p}V(p,p^{\prime\prime})$.
Clearly then, no site $p^{\prime\prime}\neq p^{\prime}$ can precede
$p^{\prime}$ in the order induced by $p$ and $x$, for otherwise the
$(p^{\prime\prime},p)$-bisector would intersect the interior of
$\operatorname{Vor}(p)$, contradicting the definition of
$\operatorname{Vor}(p)$.

The above lemma can be viewed as a natural way to generalize the ordering used
by the standard Jarvis march to the Hilbert context. In the standard setting,
consider any point $p$ on the convex hull. The Jarvis ordering is based on the
angular order of the directed segments $\overleftrightarrow{pp^{\prime}}$, for
$p^{\prime}\in P\setminus\\{p\\}$. If instead, we were to consider the left-
directed perpendicular bisectors of each of the pairs $(p^{\prime},p)$, the
ordering would be exactly the same. These left-directed perpendicular
bisectors are each directed to a point at infinity, and hence both the Jarvis
condition and ours correspond to a sorting of bisector endpoints at infinity.

We can now present our algorithm for computing the Hilbert hull. To start the
process off, we identify a pair $(p_{0},x_{0})$, where $p_{0}\in P$,
$x_{0}\in\partial\kern 1.0pt\Omega$, and $p_{0}$ lies on the boundary of an
empty ball centered at $x_{0}$. By Lemma 4.3, such a pair can be computed in
$O(n\log m)$ time. Assuming that the algorithm has generated an initial
sequence points on the hull, arriving at a pair $(p_{i},x_{i})$ satisfying the
above invariant, we compute the next pair as follows. By applying Lemma 3.4,
we compute the bisector endpoints for all the sites of
$P\setminus\\{p_{i}\\}$. Using the ordering induced by $p_{i}$ and $x_{i}$, we
take the next point $p_{i+1}$ to be the smallest in this ordering, and let
$x_{i+1}$ be the associated bisector endpoint. We repeat this until we get
back to (or go beyond) the starting pair $(p_{0},x_{0})$. The pseudocode is
presented in Algorithm 5.

Algorithm 5 Hilbert Jarvis march

procedure HilbertJarvisMarch($P$) $\triangleright$ Compute the Hilbert hull of
$P$

Compute $p_{0}\in P$ and $x_{0}\in\partial\kern 1.0pt\Omega$ such that
$B(p_{0},x_{0})$ is empty (by Lemma 4.3)

$i\leftarrow 0$

repeat

$p_{\text{best}}=\text{null}$

for $p_{\text{test}}\in P\setminus\\{p_{i}\\}$ do $\triangleright$ Find the
next site on the hull

$x_{\text{test}}\leftarrow\text{endpoint of $(p_{\text{test}},p_{i})$ bisector
(by Lemma\leavevmode\nobreak\ \ref{lem:bisector-endpoint})}$

if ($p_{\text{best}}=\text{null}$) or
($(p_{\text{test}},x_{\text{test}})<(p_{\text{best}},x_{\text{best}})$ in the
$(p_{i},x_{i})$ ordering) then

$(p_{\text{best}},x_{\text{best}})\leftarrow(p_{\text{test}},x_{\text{test}})$

end if

end for

$(p_{i+1},x_{i+1})\leftarrow(p_{\text{best}},x_{\text{best}})$

$i\leftarrow i+1$

until $(p_{i},x_{i})$ goes beyond $(p_{0},x_{0})$ $\triangleright$ We’ve
returned to the starting point

return $\langle p_{0},p_{1},\ldots p_{i}\rangle$ $\triangleright$ Return the
final hull sequence

end procedure

The correctness of this algorithm follows directly from Lemma 7.1. With each
iteration, the algorithm invokes Lemma 3.4 on each of the $n-1$ sites other
than $p_{i}$, for a running time of $O(n\log^{2}m)$. Letting $h$ denote the
number of edges on the Hilbert hull, the algorithm terminates after $h$
iterations. This yields the following result.

######  Lemma 7.3.

Given a set of $n$ points in the Hilbert geometry defined by a convex $m$-gon
$\Omega$, it is possible to construct the Hilbert hull in time
$O(nh\log^{2}m)$, where $h$ is the number of edges on the hull’s boundary.

## References

  * [1] Ahmed Abdelkader, Sunil Arya, Guilherme Dias da Fonseca, and David M. Mount.  Approximate nearest neighbor searching with non-Euclidean and weighted distances.  In Proc. 30th Annu. ACM-SIAM Sympos. Discrete Algorithms, pages 355–372, 2019.  doi:10.1137/1.9781611975482.23.
  * [2] Ahmed Abdelkader and David M. Mount.  Economical Delone sets for approximating convex bodies.  In Proc. 16th Scand. Workshop Algorithm Theory, pages 4:1–4:12, 2018.  doi:10.4230/LIPIcs.SWAT.2018.4.
  * [3] Rahul Arya, Sunil Arya, Guilherme Dias da Fonseca, and David M. Mount.  Optimal bound on the combinatorial complexity of approximating polytopes.  ACM Trans. Algorithms, 18:1–29, 2022.  doi:10.1145/3559106.
  * [4] Sunil Arya, Guilherme Dias da Fonseca, and David M. Mount.  Near-optimal $\varepsilon$-kernel construction and related problems.  In Proc. 33rd Internat. Sympos. Comput. Geom., pages 10:1–15, 2017.  doi:10.4230/LIPIcs.SoCG.2017.10.
  * [5] Sunil Arya, Guilherme Dias da Fonseca, and David M. Mount.  On the combinatorial complexity of approximating polytopes.  Discrete Comput. Geom., 58(4):849–870, 2017.  doi:10.1007/s00454-016-9856-5.
  * [6] Sunil Arya, Guilherme Dias da Fonseca, and David M. Mount.  Approximate polytope membership queries.  SIAM J. Comput., 47(1):1–51, 2018.  doi:10.1137/16M1061096.
  * [7] Stephen Boyd and Lieven Vandenberghe.  Convex Optimization.  Cambridge University Press, 2004.  doi:10.1017/CBO9780511804441.
  * [8] Madeline Bumpus, Caesar Dai, Auguste H. Gezalyan, Sam Munoz, Renita Santhoshkumar, Songyu Ye, and David M. Mount.  Software and analysis for dynamic Voronoi diagrams in the Hilbert metric, 2023.  arXiv:2304.02745.
  * [9] Yongxin Chen, Tryphon Georgiou, and Michele Pavon.  Entropic and displacement interpolation: A computational approach using the Hilbert metric.  SIAM J. Appl. Math., 76:2375–2396, 2016.  doi:10.1137/16M1061382.
  * [10] Mark de Berg, Otfried Cheong, Marc van Kreveld, and Mark Overmars.  Computational Geometry: Algorithms and Applications.  Springer, 3rd edition, 2010.  doi:10.1007/978-3-540-77974-2.
  * [11] Friedrich Eisenbrand, Nicolai Hähnle, and Martin Niemeier.  Covering cubes and the closest vector problem.  In Proc. 27th Annu. Sympos. Comput. Geom., pages 417–423, 2011.  doi:10.1145/1998196.1998264.
  * [12] Friedrich Eisenbrand and Moritz Venzin.  Approximate CVPs in time $2^{0.802n}$.  J. Comput. Sys. Sci., 124:129–139, 2021.  doi:10.1016/j.jcss.2021.09.006.
  * [13] Auguste H. Gezalyan and David M. Mount.  Voronoi diagrams in the Hilbert metric, 2021.  arXiv:2112.03056.
  * [14] Leonidas J. Guibas, Donald E. Knuth, and Micha Sharir.  Randomized incremental construction of Delaunay and Voronoi diagrams.  Algorithmica, 7:381–413, 1992.  doi:10.1007/BF01758770.
  * [15] D. Hilbert.  Ueber die gerade Linie als kürzeste Verbindung zweier Punkte.  Math. Annalen, 46:91–96, 1895.  doi:10.1007/BF02096204.
  * [16] S. Kullback and R. A. Leibler.  On information and sufficiency.  Annals. Math. Stat., 22:79–86, 1951.  doi:10.1214/aoms/1177729694.
  * [17] Bas Lemmens and Roger Nussbaum.  Birkhoff’s version of Hilbert’s metric and its applications in analysis, 2013.  arXiv:1304.7921.
  * [18] Márton Naszódi and Moritz Venzin.  Covering convex bodies and the closest vector problem.  Discrete Comput. Geom., 67:1191–1210, 2022.  doi:10.1007/s00454-022-00392-x.
  * [19] Frank Nielsen and Laetitia Shao.  On balls in a Hilbert polygonal geometry.  In Proc. 33rd Internat. Sympos. Comput. Geom., pages 67:1–67:4, 2017.  (Multimedia contribution).  doi:10.4230/LIPIcs.SoCG.2017.67.
  * [20] Frank Nielsen and Ke Sun.  Clustering in Hilbert’s projective geometry: The case studies of the probability simplex and the elliptope of correlation matrices.  In Frank Nielsen, editor, Geometric Structures of Information, pages 297–331. Springer Internat. Pub., 2019.  doi:10.1007/978-3-030-02520-5_11.
  * [21] Frank Nielsen and Ke Sun.  Non-linear embeddings in Hilbert simplex geometry, 2022.  arXiv:2203.11434.
  * [22] Athanase Papadopoulos and Marc Troyanov.  Handbook of Hilbert geometry, volume 22 of IRMA Lectures in Mathematics and Theoretical Physics.  European Mathematical Society Publishing House, 2014.  doi:10.4171/147.
  * [23] David Reeb, Michael J. Kastoryano, and Michael M. Wolf.  Hilbert’s projective metric in quantum information theory.  J. Math. Physics, 52(8), 2011.  doi:10.1063/1.3615729.
  * [24] Thomas Rothvoss and Moritz Venzin.  Approximate CVP in time $2^{0.802n}$ – Now in any norm!  In Proc. 23rd Internat. Conf. on Integ. Prog. and Comb. Opt. (IPCO 2022), pages 440–453, 2022.  doi:10.1007/978-3-031-06901-7_33.
  * [25] Godfried T. Toussaint.  The relative neighbourhood graph of a finite planar set.  Pattern Recogn., 12:261–268, 1980.  doi:10.1016/0031-3203(80)90066-7.
  * [26] Constantin Vernicos.  On the Hilbert geometry of convex polytopes.  In Handbook of Hilbert geometry, volume 22 of IRMA Lectures in Mathematics and Theoretical Physics, pages 111–126. European Mathematical Society Publishing House, 2014.  doi:10.48550/arXiv.1406.0733.
  * [27] Constantin Vernicos and Cormac Walsh.  Flag-approximability of convex bodies and volume growth of Hilbert geometries, 2018.  arXiv:1809.09471.