\section{\bf Introduction}\label{sec:Intro}
\subsection{Background}
Classical objects in computational geometry, such as polytopes, 
arrangements, and Delaunay triangulations are defined by explicit relations 
\cite{ComputeGeoBook_BCKO,Edelsbrunner-book-1987,Mulmuley1994}. 
A few years ago an interesting family of geometric objects defined by implicit 
relations was introduced in the pioneering works of T. Asano, J. Matou{\v{s}}ek 
and T. Tokuyama \cite{AMT2,AMTn}. An important member in this family is a zone diagram. 
In order to understand this object better, consider first the 
more familiar concept of a Voronoi diagram, another classical object which is defined explicitly. 
Given a set $X$, a distance function $d$, and a collection of subsets $(P_k)_{k\in K}$ in $X$  (called the sites or the  generators),  we associate with each site $P_k$ a corresponding Voronoi cell, that is, the set $R_k$ of all $x\in X$ whose distance to $P_k$ is not
greater than its distance to the union $\bigcup_{j\neq k}P_j$. The collection $(R_k)_{k\in K}$ 
of Voronoi regions is the Voronoi diagram. 
On the other hand, in the case of a zone diagram $R=(R_k)_{k\in K}$ we associate with each site $P_k$ the set $R_k$ of all $x\in X$ whose distance to $P_k$ is not greater than its distance to the union of the other sets $R_j$, $j\neq k$. Figures \ref{fig:IntroVoronoi} and  \ref{fig:IntroZone} show the Voronoi and zone diagrams, respectively, corresponding to the same ten singleton sites in the Euclidean plane.
\begin{figure}[t]
\begin{minipage}[t]{0.45\textwidth}
\begin{center}
{\includegraphics[scale=0.57]{2D-10Sites-l2-0001-New.png}
\end{center}
 \caption{The Voronoi diagram of 10 point sites in a square in the Euclidean plane.}
\label{fig:IntroVoronoi}
\end{minipage}
\hfill
\begin{minipage}[t]{0.48\textwidth}
\begin{center}
{\includegraphics[scale=0.57]{ZD-5Iterations-2D-10Sites-l2-0005-New.png}}
\end{center}
 \caption{The zone diagram of the same 10 points as in
Figure \ref{fig:IntroVoronoi}.}
\label{fig:IntroZone}
\end{minipage}
\end{figure}

At  first glance it seems that the definition of a zone diagram is
circular, because the definition of each $R_k$ depends on $R_k$ itself via 
the definition of the other cells  $R_j$, $j\neq k$. On a second thought, we
see that, in fact, a zone diagram is defined to be a fixed point of a
certain mapping (called the $\textnormal{Dom}$ mapping), that is, a solution of the equation $R=\textnormal{Dom}(R)$. While the Voronoi diagram is explicitly
defined and, hence, its existence and uniqueness are obvious, neither the existence nor the uniqueness of a zone diagram are obvious in advance. 
In addition, even if some existence (or uniqueness) results are proved, one still faces the problem of finding algorithms for computing
zone diagrams.

In the original works \cite{AMT2,AMTn}, zone diagrams were introduced and 
studied in the case where $(X,d)$ was the Euclidean plane, each site $P_k$ was a
single point, and all these (finitely many) points were different. The existence 
and uniqueness of a zone diagram was proved in this case, and an iterative 
algorithm for approximating it was suggested. The existence of a zone diagram is not known in general but it is known in 
a few particular cases, some of them involve quite mild assumptions. 
 For instance, in   \cite{ReemReichZone} it was shown that a zone diagram of two general sites 
 in any metric space always exists. In fact, the proof holds in a more  general setting 
 called $m$-spaces, in which $X$ is an arbitrary nonempty set and the ``distance'' function $d$ 
 should only satisfy the condition $d(x,x)\leq d(x,y)\,\,\,\,
\forall x,y\in X$ and can take any value in the interval
$[-\infty,\infty]$. Simple examples given there show that in general uniqueness of the zone 
diagram does not necessarily hold.  In \cite{KopeckaReemReich} it is shown, in particular, 
that a zone diagram of any finite  number  of sites exists, assuming these sites are compact 
and positively  separated (that is, there is a positive lower bound on the distance between 
any pair of them) and that they are located in a large compact subset of $\mathbb{R}^n$ with a strictly 
convex norm (the exact  result is more general). Another existence (and also uniqueness) 
result is discussed in \cite{KMT2012}. Here the setting is the $n$-dimensional Euclidean space 
$\mathbb{R}^n$, or, more generally, finite dimensional normed spaces which are both strictly convex 
and smooth  \cite{KMT2012}. The sites are again positively separated. 



In addition to zone diagrams, other implicit geometrical objects were introduced and studied, 
partly in order to understand better zone diagrams. One such an object is called a double zone diagram  
\cite{ReemReichZone}. Formally, it is  defined as the fixed point of the second iteration of 
$\textnormal{Dom}$, the mapping which defines zone diagrams. As proved in \cite{ReemReichZone}, double zone diagrams 
exist in a relatively general setting (in any $m$-space, for infinitely many arbitrary sites). 
The proof is based on a nonconstructive argument (the Knaster-Tarski fixed point theorem \cite{Knaster,Tarski}). 
As a matter of fact, it was shown that there exist a least and a greatest double zone diagrams.   
Because of the nonconstructive nature of the proof, no general way was suggested to compute them. 
The importance of double zone diagrams to the computation of zone diagrams will become clear later. 
Additional implicit objects related to zone diagrams are  trisectors \cite{AsanKirk,AMT2,AsanoTokuyama, ChunOkadeTokuyama2010}, $k$-sectors and $k$-gradations \cite{ImaiKawamuraMatousekReemTokuyamaCGTA}, 
and territory diagrams  \cite{DeBiasiKalantaris2011} (also called mollified zone diagrams). 


 

One of the main challenges regarding zone diagrams is their computation. This is evident already 
from the title of the original work \cite{AMTn}. 
In a continuous setting this task has been addressed  so far only by  Asano, Matou{\v{s}}ek, and Tokuyama 
\cite{AMTn} in the case of the Euclidean plane with finitely many point sites. While the formal 
claims and proofs mentioned there are insightful, 
the discussion about an actual method for the approximate computation of zone diagrams 
was very brief with almost no theoretical or  practical details. More precisely, it was written there 
that  one can use convex polygons with many sides for approximating the components of these sequences, 
and that each iteration is computationally  demanding, but with the exception of a few (interesting) 
pictures, no additional details were given. This (approximate) computation was restricted to the above setting. 
As will be explained in Section  \ref{sec:QualitativeDescription}, one of  the main reasons for this restriction 
was the lack of a method for computing Voronoi diagrams in a general setting. Another 
reason is the lack of any known representation of the boundaries of the involved regions (conjectured to 
be non-algebraic in many cases; this is supported by the recent preprint \cite{MonterdeOngayPre2013}). 

In a discrete setting ($X$ is a finite set of points) there has been a limited discussion in the case of two sites: a point and a ``line segment'' or a ``curve'' in the digital Euclidean plane \cite{COT} (it was claimed there on the last page that the algorithm can be generalized to finitely many line segments but no details were given)  and two finite sites in any finite $m$-space \cite{ReemReichZone}. The actual computation in the first case is restricted to the above setting and no convergence    analysis was given; time complexity analysis was quite brief. In the second case the actual computation is by brute force. 
In addition, it is  not clear in which sense the resulting discrete zone diagrams approximate the continuous ones, although  intuitively they may approximate them with respect to the Hausdorff distance. The computation of double zone diagrams was discussed briefly in \cite{ReemReichZone} only in a discrete setting (finitely many finite sites in a  finite $m$-space).


As mentioned above, a few additional geometrical objects have been introduced and studied. 
In particular, some of them have been studied from the computational point of view, 
as in the case of trisectors and $k$-sectors \cite{AsanKirk,AMT2,AsanoTokuyama, ChunOkadeTokuyama2010,ImaiKawamuraMatousekReemTokuyamaCGTA} and territory diagrams  \cite{DeBiasiKalantaris2011}. 
The method of computation is approximate. One can find interesting pictures and ideas in these works, 
but in many cases a corresponding precise analysis of convergence or of the level of approximation or 
of time complexity is very brief or absent. It is interesting to note however that 
 the main computational ingredient suggested in this work, namely the algorithm for computing 
 Voronoi diagrams in a general setting (see below for more details), enables the computation 
 of some of the above objects due to the approximating schemes suggested in the corresponding works. 
 This issue will not be dealt with here, but see, e.g.,  \cite[p. 719]{ImaiKawamuraMatousekReemTokuyamaCGTA} 
 for a short related discussion. 

 
\subsection{Contribution of this paper} 
 This paper considers the question of computing zone diagrams in a wide class of metric spaces, called 
  proper geodesic metric spaces which have the geodesic inclusion property. 
 This class, which is introduced here, includes, in particular, Euclidean spheres (of any dimension) 
 and finite  dimensional strictly convex spaces. 
 The considered sites can be any positively separated closed sets (possibly infinitely many). 
 It is shown that a generalization of the iterative algorithm suggested by Asano, Matou{\v{s}}ek, and Tokuyama  converges to the least and the greatest double zone diagrams. In various cases  a zone diagram can be obtained from the resulting double zone diagrams, and as a matter of fact, in many of these cases  (such as the Euclidean case) the limit is the unique zone diagram.  In the normed case the suggested way to compute approximately the corresponding iterations and the resulting (double)  zone diagram is by using a new  algorithm for computing Voronoi diagrams \cite{ReemISVD09, ReemPhD} which enables their approximate computation in any normed space (possibly infinite dimensional) and with sites of a general form. Many pictures of (double) zone diagrams or their  approximating sets produced by the method described here are given. Since the algorithm for computing  Voronoi  diagrams enables the computation of each cell independently of the other ones, the above iterative algorithm can be easily  parallelized. Finally, certain interesting topological properties of Voronoi cells in the above mentioned setting are derived along the way, e.g., that their boundaries must coincide with their bisectors. 

\subsection{The structure of this paper} In Section \ref{sec:definitions}
 the basic definitions and notation are presented. In Section \ref{sec:Geodesic} 
geodesic metric spaces and particular classes of them are discussed. In Section \ref{sec:QualitativeDescription} 
 a qualitative discussion regarding the iterative algorithm is given. In Section \ref{sec:ConvergenceResult} 
 the main convergence result is presented. In Section \ref{sec:ActualCompute} a few 
 theoretical and practical issues are discussed, among them the issue of time complexity and some 
 details regarding the actual computation of the corresponding  iterations. 
 The topological properties of Voronoi cells which follow as a by-product 
 of the approach given here are described in Section \ref{sec:TopologicalProperties}. 
 The paper is concluded in Section \ref{sec:end}, which contains several interesting  questions and unexplained phenomena. 
  Proofs of the convergence theorems and related claims are given in Section \ref{sec:proofs}. 
\section{Notation and definitions}\label{sec:definitions}

This section presents the notation and basic definitions. The important 
concepts of geodesic metric spaces and a few particular classes of such spaces 
are discussed in Section \ref{sec:Geodesic}. 


Throughout the text we will make use of tuples, the components of which are sets (which are subsets of the given set $X$). Every operation or relation between such tuples, or on a single tuple, is done component-wise. Hence, for example, if $K\neq \emptyset$ is a set of indices, and if $R=(R_k)_{k\in K}$ and $S=(S_k)_{k\in K}$ are two tuples  of sets, then $R\bigcap S=(R_k\bigcap S_k)_{k\in K}$,  $\overline{R}=(\overline{R_k})_{k\in K}$, and $R\subseteq S$ means $R_k\subseteq S_k$ for each $k\in K$. When $R$ is a tuple, the  notation $(R)_k$ is the $k$-th component of $R$, i.e, $(R)_k=R_k$. 
When dealing with a subset $A$ of a given metric space $(X,d)$, we denote by 
$\overline{A}$, $\textnormal{int}(A)$ and $\partial A$ its closure, interior, and the  boundary, respectively.
We note that the definitions given below can be generalized almost word for 
word to the case where $(X,d)$ is a general  metric space and sometimes even beyond (e.g., $m$-spaces \cite{ReemReichZone}).
However, in practice $X$ will often be nothing more than a large rectangle or box in $\mathbb{R}^n$ 
and the metric will be induced by a given norm. 


\begin{defin}\label{def:dom}
Let $(X,d)$ be a metric space. Given two nonempty sets $P,A\subseteq X$, the dominance region
$\textnormal{dom}(P,A)$ of $P$ with respect to $A$ is the set of all $x\in X$
the distance of which to $P$ is not greater than its distance to $A$,
that is, \begin{equation*}
\textnormal{dom}(P,A)=\{x\in X: d(x,P)\leq d(x,A)\}.
\end{equation*}
Here $d(x,A)=\inf\{d(x,a): a\in A\}$.
\end{defin}
\begin{defin}\label{def:Voronoi}
Let $(X,d)$ be a metric space. Let $K$ be a set of at least 2 elements (indices), possibly
infinite. Given a  tuple $(P_k)_{k\in K}$ of nonempty subsets
$P_k\subseteq X$, called the generators or the sites, the Voronoi diagram  induced by this tuple is the tuple $(R_k)_{k\in K}$ of nonempty subsets
$R_k\subseteq X$, such that for all $k\in K$,
\begin{equation*}
R_k=\textnormal{dom}(P_k,{\underset{j\neq k}\bigcup P_j})=\{x\in X: d(x,P_k)\leq d(x,P_j)\,\,\forall j\neq k,\,j\in K \}.
\end{equation*}
 In other words, each $R_k$, called a Voronoi cell or a Voronoi region, is the set of
all $x\in X$ whose distance to the site $P_k$ is not greater than 
its distance to any other site $P_j$, $j\neq k$. 
\end{defin}
\begin{defin}\label{def:zone}
Let $(X,d)$ be a metric space. Let $K$ be a set of at least 2 elements (indices), possibly
infinite. Given a  tuple $(P_k)_{k\in K}$ of nonempty subsets
$P_k\subseteq X$, a zone diagram with respect to
that tuple is a tuple $R=(R_k)_{k\in K}$ of nonempty subsets
$R_k\subseteq X$ such that
\begin{equation*}
R_k=\textnormal{dom}(P_k,\textstyle{\underset{j\neq k}\bigcup R_j})\quad \forall k\in K.
\end{equation*}
In other words, if one defines $X_k=\{C: P_k\subseteq C\subseteq X\}$, then a zone diagram
is a fixed point of the mapping $\textnormal{Dom}:\underset{{k\in K}}\prod
X_k\to \underset{{k\in K}}\prod
X_k$, defined by
\begin{equation}\label{eq:TZoneDef}
\textnormal{Dom}(R)=(\textnormal{dom}(P_k,\textstyle{\underset{j\neq k}\bigcup R_j}))_{k\in K}.
\end{equation}
A tuple $R=(R_k)_{k\in K}$ is called a double zone diagram if it is a fixed point of the second iteration $\textnormal{Dom}\circ \textnormal{Dom}$, i.e., $R=\textnormal{Dom}^2(R)$.
\end{defin}
Note that if $R$ is a zone diagram, then it is a double zone diagram as implied by iterating the equation
 $R=\textnormal{Dom}(R)$ with $\textnormal{Dom}$. 

\begin{defin}\label{def:proper}
A metric space $(X,d)$ is called proper, or finitely compact, if closed balls in $X$ are compact, or, equivalently, 
 any bounded sequence has a convergent subsequence. 
\end{defin}
Typical examples of proper metric spaces are finite dimensional normed space and compact metric spaces. 
\begin{defin}\label{def:Hausdorff}
Given two nonempty sets $A_1,A_2$ in a metric space $(X,d)$, the Hausdorff distance between them is defined by
\begin{equation*}
D(A_1,A_2)=\max\{\sup_{a_1\in A_1}d(a_1,A_2),\sup_{a_2\in A_2}d(a_2,A_1)\}.
\end{equation*}
\end{defin}
Recall that the Hausdorff distance is different from the usual distance between two sets which is defined by 
\begin{equation*}
d(A_1,A_2)=\inf\{d(a_1,a_2): a_1\in A_1,\,a_2\in A_2\}. 
\end{equation*}

\section{Geodesic metric spaces of various types}\label{sec:Geodesic}
This section discusses geodesic metric spaces of various types, and in particular 
it introduces  the type which is important for the rest of the paper, namely 
geodesic metric spaces which have the geodesic inclusion property. 
\begin{defin}\label{def:GeodesicMetric}
Let $(X,d)$ be a metric space. 
\begin{enumerate}[(a)]
\item \label{def:GeodesicMetric:GeodesicSegment} 
Let $x,y\in S\subseteq X$.  The subset $S$ is
called a geodesic segment (or a metric segment) between $x$ and $y$  if there exists
 an isometric function $\gamma$ (that is, a distance preserving mapping) which maps 
 a  real line segment $[r_1,r_2]$ onto $S$ and satisfies  $\gamma(r_1)=x$ and  $\gamma(r_2)=y$. 
 We denote $S=[x,y]_{\gamma}$, or simply $S=[x,y]$. 
 If between all points $x,y\in X$ there exists a geodesic segment, 
 then $(X,d)$ is called  a  geodesic metric space. The sets $[x,y)=[x,y]\backslash\{y\}$, 
 $(x,y]=[x,y]\backslash\{x\}$, and $(x,y)=[x,y]\backslash\{x,y\}$ represent the half open segment and open  
 segments respectively. 
\item\label{def:GeodesicMetric:GeodesicInclusion} A geodesic metric space $(X,d)$ is said to have 
the geodesic inclusion property if 
a nontrivial intersection between two different geodesics can have no strict bifurcation   point. 
More precisely, the following holds: given $u,v,b,z\in X$, if  $b\in [u,z)_{\gamma_1}$, and $b\in [v,z)_{\gamma_2}$, then either $u\in [v,z]_{\gamma_2}$ or $v\in [u,z]_{\gamma_1}$. 
\end{enumerate}
\end{defin}

An illustration of Definition \ref{def:GeodesicMetric}\eqref{def:GeodesicMetric:GeodesicInclusion} is 
given in Figure \ref{fig:GeodesicInclusion}. Note that a simple consequence of the definition is  
that actually either $u\in [v,b]_{\gamma_2}$ or $v\in [u,b]_{\gamma_1}$ because if for 
instance $u\in (b,z]_{\gamma_2}$, then this and $b\in [u,z]_{\gamma_1}$ imply $d(u,z)<d(b,z)\leq d(u,z)$. 

\begin{figure}
\begin{minipage}[t]{1\textwidth}
\begin{center}
{\includegraphics[scale=0.51]{GeodesicInclusionPropertyNew.png}
\end{center}
 \caption{An illustration of the geodesic inclusion property (Definition \ref{def:GeodesicMetric}\eqref{def:GeodesicMetric:GeodesicInclusion}).}
\label{fig:GeodesicInclusion}
\end{minipage}
\end{figure}

\begin{expl}
 Simple and familiar examples of geodesic metric spaces are: the Euclidean plane, 
 any convex subset of a normed space [the geodesics are line segments, but in the non-strictly 
 convex case other types of geodesics exists, e.g., ``zigzag'' ones: 
 for instance, $([0,1]\times \{0\})\cup (\{1\}\times [0,1])$ connecting 
 $x=(0,0)$ and  $y=(1,1)$ in $(\mathbb{R}^2,\ell_1)$],  
 Euclidean spheres (a geodesic segment between two points is the short arc on the great circle on which the 
points are located), complete  Riemannian manifolds \cite[pp. 25-28]{Jost} (in particular, compact ones), 
and hyperbolic spaces \cite[pp. 538-9]{ReichShafrir}. 
\end{expl}
\begin{expl}
The geodesic inclusion property holds for Euclidean 
spheres of any dimension (because a nontrivial intersection between two geodesics, namely 
parts of large circles, can happen only when they coincide at the poles or when they 
are identical; but then the points $u,v$, and $b$ from Definition \ref{def:GeodesicMetric} coincide). 
It also holds (see Proposition \ref{prop:GeodesicInclusionStrictlyConvex}) for the class of normed spaces called strictly convex. These spaces have the property that  $|x+y|<|x|+|y|$ holds 
whenever $x$ and $y$ are arbitrary elements in the space which are  not on the 
 same ray (i.e., $x\neq 0, y\neq 0$, and $x/|x|\neq y/|y|$).  
 Equivalently,  the unit sphere of the space does not contain any line segment. 
 In particular, $\mathbb{R}^d$ with the $\ell_p$ norm, $p\in (1,\infty)$ is a typical example of a (proper) 
 strictly convex geodesic spaces and $\mathbb{R}^d$ with the $\ell_1$ or $\ell_{\infty}$ norms are 
  typical examples of  spaces which are not strictly convex. See  \cite{LindenTzafriri, Prus2001} for more details. 
 The geodesic inclusion property is always absent in any normed space 
which is not strictly convex: in each such a space (whose dimension must be at least 2)  
one can find 3 points $b,c,z$ such that $c\notin [b,z]$ and $|z-b|=|z-c|+|c-z|$; now one 
takes a point $u\notin [b,z]$ but on the ray emanating from $z$ in the direction of $b$ 
and another point $v\notin [b,c]$ but on the ray emanating from $c$ in the direction of $b$;  
 now one observes that $[u,b]\cup [b,z]$ and $[v,b]\cup [b,c]\cup [c,z]$ are two 
geodesic segments whose intersection includes a nontrivial bifurcation point ($b$).  
\end{expl}

\section{A qualitative description of the algorithm}\label{sec:QualitativeDescription}
As mentioned before, a tuple $R=(R_k)_{k\in K}$ of cells (subsets) is a zone diagram if it satisfies the fixed point equation $R=\textnormal{Dom}(R)$. A common and natural approach in fixed point theory for the computation of a fixed point of a given mapping $f$ is to use iterations  \cite{GranasDugondji,HB}. One starts with some point $y_0$ in the space $Y$ on which $f$  is defined, and starts iterating $f$. A sequence $y_1=f(y_0), \ldots, y_{n+1}=f(y_n),\ldots$ is generated, and one hopes that it converges in some sense to a fixed point of $f$. In general convergence is not guaranteed (just take a point 
on the unit circle and apply iteratively a rotation operator on it), but under some assumptions on 
the mapping and the space corresponding convergence results can be proved.   

Returning to our setting, the given mapping is the $\textnormal{Dom}$ mapping. The given  space on which it is defined is $Y=\prod_{k\in K}
X_k$, where $X_k=\{C: P_k\subseteq C\subseteq X\}$ for all $k\in K$. In other words, $Y$ is the collection of all tuples (vectors) whose $k$-th component is a subset $C$ of the given world $X$ such that $C$ contains the site $P_k$. As a result of the above,  a natural choice for the starting point $y_0\in Y$ is the collection $(P_k)_{k\in K}$ of the given sites. 

It turns out that it is convenient to denote $I^{(0)}:=(P_k)_{k\in K}=y_0$,  $O^{(0)}:=\textnormal{Dom} (I^{(0)})$, and to define inductively  the inner sequence $I^{(n)}:=\textnormal{Dom}(O^{(n-1)})$ and the outer one  $O^{(n)}:=\textnormal{Dom} (I^{(n)})$ for all integers  $n\geq 1$. The reason behind this choice is explained in the following simple lemma, a large part of it was observed in  \cite{AMTn} 
and a small part in \cite{ReemPhD,ReemReichZone}.

\begin{lem}\label{lem:DomPk} 
\begin{enumerate}[(a)
Let $(P_k)_{k \in K}$ be a tuple of  nonempty and closed subsets in $X$. 
\item $\textnormal{Dom}$ is antimonotone, i.e., $\textnormal{Dom}(R)\subseteq \textnormal{Dom}(S)$ whenever $S\subseteq R$; $\textnormal{Dom}^2$ is monotone, that is, $R\subseteq S\Rightarrow \textnormal{Dom}^2(R)\subseteq \textnormal{Dom}^2(S)$. 
\item  The sequence $(I^{(n)})_{n=0}^{\infty}$ is increasing, $(O^{(n)})_{n=0}^{\infty}$ is decreasing 
 and 
\begin{equation}\label{eq:InOn}
I^{(n)}\subseteq \bigcup_{\gamma=0}^{\infty}\textnormal{Dom}^{2\gamma}(P_k)_{k\in K}\subseteq \bigcap_{\gamma=0}^{\infty}\textnormal{Dom}^{2\gamma+1}(P_k)_{k\in K}\subseteq O^{(n)}
\end{equation}
for each nonnegative integer $n$. In addition, if $R$ is a zone or double zone diagram in $X$, then
\begin{equation}\label{eq:DomIterations}
I^{(n)}\subseteq\bigcup_{\gamma=0}^{\infty}\textnormal{Dom}^{2\gamma}(P_k)_{k\in K}\subseteq R \subseteq \bigcap_{\gamma=0}^{\infty}\textnormal{Dom}^{2\gamma+1}(P_k)_{k\in K}\subseteq O^{(n)}\quad \forall n\geq 0.
\end{equation}
\end{enumerate}
\end{lem}


Figures  \ref{fig:Zone2} and \ref{fig:Zone3} show $I^{(2)},O^{(2)}$ and $I^{(3)},O^{(3)}$ respectively, in a square in $(\mathbb{R}^2,\ell_2)$. The sites are as in Figure  \ref{fig:IntroZone} and thus these figures approximate 
the zone diagram of Figure \ref{fig:IntroZone}. For another example of these iterations, see Figure  \ref{fig:ZoneLinfty}.

\begin{figure}
\begin{minipage}[t]{0.45\textwidth}
\begin{center}
{\includegraphics[scale=0.27]{ZD-2Iterations-2D-10Sites-l2-0005-New.png}
\end{center}
 \caption{Approximation of the zone diagram of Figure  \ref{fig:IntroZone} using the boundaries of  $I^{(2)}$ 
 and $O^{(2)}$ (1600 rays in  each iteration).}
\label{fig:Zone2}
\end{minipage}
\hfill
\begin{minipage}[t]{0.45\textwidth}
\begin{center}
{\includegraphics[scale=0.27]{ZD-3Iterations-2D-10Sites-l2-0005-New.png}}
\end{center}
 \caption{The setting of Figure \ref{fig:Zone2}, but with $I^{(3)}$ and $O^{(3)}$ (1600 rays in  each iteration). }
\label{fig:Zone3}
\end{minipage}

\begin{minipage}[t]{0.45\textwidth}
\begin{center}
{\includegraphics[scale=0.55]{ZD-5Iterations-2D-10Sites-l2-01-Rays-New.png}}
\end{center}
 \caption{Approximation of the zone diagram of Figure \ref{fig:IntroZone} (80 rays in 
 each iteration). }
\label{fig:Zone80}
\end{minipage}
\hfill
\begin{minipage}[t]{0.45\textwidth}
\begin{center}
{\includegraphics[scale=0.55]{ZD-5Iterations-2D-10Sites-l2-005-Rays-New.png}}
\end{center}
 \caption{Approximation of the zone diagram of Figure \ref{fig:IntroZone} (160 rays in 
 each iteration).}
\label{fig:Zone160}
\end{minipage}
\end{figure}


From the lemma it can be seen that the goal of $I^{(n)}$ is to approximate the zone diagram from below, and the goal of $O^{(n)}$ is to approximate it from above.   However, at this stage several difficulties arise. First, it is not clear that these sequences converge, and if they do converge, it is not clear that both of them converge to the same limit and that the limit of at least one of them is indeed a zone diagram.  

Second, even if a convergence result is obtained, one faces the problem of the actual computation of $I^{(n)}$ and $O^{(n)}$ for $n\geq 1$. Even for $n=0$ this is not a simple task since $O^{(0)}=\textnormal{Dom}(I^{(0)})=(\textnormal{dom}(P_k,\bigcup_{j\neq k}P_j))_{k\in K}$ is the Voronoi diagram of the sites $(P_k)_{k\in K}$, and hence, for sites of a general form, or sites which possibly form a    degenerate configuration, or for space with a general norm, one has to find algorithms which  enable the computation of Voronoi diagrams in such a setting, and most of the familiar algorithms for computing Voronoi diagrams (see e.g., \cite{Aurenhammer,AurenhammerKlein,ComputeGeoBook_BCKO,Edelsbrunner-book-1987,Mulmuley1994,OBSC,SharirAgarwal} for some  reviews) are not able to achieve this task (or, in some cases, 
they are too complicated or too  slow).  Even if one restricts oneself to the familiar case of the Euclidean plane with point sites, one still faces problems starting from the iteration $n=1$. This is because one has to  know how to compute the components of $\textnormal{Dom}(S)$ for a given tuple $S=(S_k)_{k\in K}$, and hence one has to know a method for computing  $\textnormal{dom}(P_k,\bigcup_{j\neq k}S_j)$, i.e., the Voronoi cell of $P_k$ with respect to the set $\bigcup_{j\neq k}S_j$. Unfortunately, when $S=I^{(n)}$, $n\geq 1$, or $S=O^{(n)}$, $n\geq 0$, the components of $S$ are general sets, and again, it seems that most of the familiar algorithms for computing Voronoi diagrams  are not helpful here. 

The sequences $(I^{(n)})_{n=0}^{\infty}$ and $(O^{(n)})_{n=0}^{\infty}$ were introduced in \cite{AMTn} in the case of the Euclidean plane with finitely many point sites. It was shown that $\overline{\bigcup_{n=0}^{\infty} I^{(n)}}=\bigcap_{n=0}^{\infty}O^{(n)}$ and that this tuple is the unique zone diagram. 
The proof of this result is not at all obvious and it contains some useful ideas (see Section 
~\ref{sec:ConvergenceResult} below). Although this is not a pure  convergence result (no limits) and although 
no error estimates for the level of approximation were given, one 
still has the pleasant phenomenon that $(I^{(n)})_{n=0}^{\infty}$ increases to a zone diagram and $(O^{(n)})_{n=0}^{\infty}$ decreases to it.  As for the actual computation of $I^{(n)}$ and $O^{(n)}$, it was already mentioned in Section \ref{sec:Intro} that very few details were given  there. 

\section{The convergence theorems}\label{sec:ConvergenceResult}
As mentioned in Section \ref{sec:QualitativeDescription}, it was proved in \cite{AMTn} that in the case of the Euclidean plane with point sites one has the equality $m:=\overline{\bigcup_{n=0}^{\infty}I^{(n)}}=\bigcap_{n=0}^{\infty}O^{(n)}:=M$, and $m=M$ is the unique zone diagram. As Theorem ~\ref{thm:Metric} below shows, it is always true that $M=\textnormal{Dom}(m)$ 
in any metric space and for arbitrary sites. It is not known that $m=\textnormal{Dom}(M)$ in general, but, 
as Theorem ~\ref{thm:AlgUniConvex} shows,  under the assumption that the space is a  geodesic metric space with (possibly infinitely many general)  positively separated sites, the equality $m=\textnormal{Dom}(M)$ holds. However, in this latter case both $m$ and $M$ are double zone diagrams rather than  zone diagrams and hence they are not necessarily equal 
(see Example~ \ref{ex:mM}). 
In addition, Theorem ~\ref{thm:AlgUniConvex} also discusses another way to obtain a zone diagram in the special case of two sites.

A corollary to the theorem ensures that  when the world $X$ is compact, then $(I^{(n)})_{n=0}^{\infty}$ and $(O^{(n)})_{n=0}^{\infty}$ converge to $m$ and $M$ respectively  with respect to the Hausdorff distance. As a result, since in practice $X$ is taken to be compact, if $I^{(n)}$ and $O^{(n)}$ are shown experimentally to be almost the same for some $n$, then one has a good approximation to both $m$ and $M$. When it is known that a zone diagram exists, 
e.g., under the assumptions of \cite{KopeckaReemReich}, then one has a good approximation to this zone diagram, and in fact this also shows that probably this zone diagram is unique and coincides with both $m$ and $M$.  Another corollary to the above ensures that whenever the least and the greatest double zone diagrams coincide, then $m=M$ and  they both coincide with the unique zone diagram. This is true in particular in the Euclidean case.  

\begin{thm}\label{thm:Metric}
Let $(X,d)$ be a metric space and let $P=(P_k)_{k\in K}$ be a tuple of nonempty sets in $X$. 
Let
\begin{equation}\label{eq:mM}
M=\!\bigcap_{\gamma=0}^{\infty}\textnormal{Dom}^{2\gamma+1}(P)=\bigcap_{n=0}^{\infty}O^{(n)}, \quad m\!=\!\overline{\bigcup_{\gamma=0}^{\infty}\textnormal{Dom}^{2\gamma}(P)}=\overline{\bigcup_{n=0}^{\infty}I^{(n)}}.
\end{equation}
Then $M=\textnormal{Dom}(m)$.
\end{thm}

\begin{thm}\label{thm:AlgUniConvex}
Let $(X,d)$ be a proper metric space which has the geodesic 
inclusion property. Suppose that $(P_k)_{k\in K}$ is a given tuple 
of closed subsets of $X$ satisfying 
\begin{equation}\label{eq:d(P_k,P_j)>0}
 \inf\{d(P_k,P_j): j,k\in K,\, j\neq k\}>0. 
\end{equation}
Let $m$ and $M$ be defined by \eqref{eq:mM}. 
Then $M=\textnormal{Dom}(m)$, $m=\textnormal{Dom}(M)$, and $m$ and $M$ are, respectively, the least  and  greatest double zone diagrams. 
In addition, if $|K|=2$, then by letting  $m=(m_1,m_2)$ and $M=(M_1,M_2)$, each of the 
pairs $(m_1,M_2)$ and $(M_1,m_2)$ is a zone diagram. 
\end{thm}
\begin{cor}\label{cor:CompactLimit}
Under the setting of Theorem \ref{thm:AlgUniConvex}, let 
$M=(M_k)_{k\in K}$ and $m=(m_k)_{k\in K}$ be defined by \eqref{eq:mM}. Assume also 
that $(X,d)$ is compact. Then  
\begin{equation}\label{eq:Mm_k}
M_k=\lim_{\gamma\to\infty}(\textnormal{Dom}^{2\gamma+1}(P))_k
\quad \textnormal{and} \quad
m_k\!=\lim_{\gamma\to\infty}(\textnormal{Dom}^{2\gamma}(P))_k\quad \forall k\in K
\end{equation}
with respect to the Hausdorff distance. 
\end{cor}
\begin{cor}\label{cor:ZoneDoubleZone}
Under the setting of Theorem \ref{thm:AlgUniConvex}, if it is known that 
 the least and the greatest double zone diagrams coincide, then $m=M$ and they both coincide with the unique 
  zone  diagram. In particular this is true when $(X,d)$ is a finite dimensional Euclidean space. 
\end{cor}


The proof of the above mentioned theorems and corollaries is quite long and technical and it can be found in Section ~\ref{sec:proofs}.  
The proof of  Corollary ~\ref{cor:CompactLimit} actually follows from a quite general argument not related to zone diagrams. As already mentioned, the proof given in \cite{AMTn} for the 2D Euclidean case of point sites is 
far from being immediate and it contains certain useful ideas which can be modified to the 
setting considered here. In particular, this is true for \cite[Lemma 3.1, Lemma ~5.1]{AMTn}. 
However, the generalization of some of the arguments given there to the setting of this paper is definitely 
not immediate and one has to pay attention to certain subtle points, among them the verification of the equality  
\begin{equation}\label{eq:UnionIntersectBody}
\bigcup_{j\neq  k}\bigcap_{\gamma=1}^{\infty}(\textnormal{Dom}^{2\gamma+1}(P_k)_{k\in K})_j=\bigcap_{\gamma=1}^{\infty}\bigcup_{j\neq k} (\textnormal{Dom}^{2\gamma+1}(P_k)_{k\in K})_j 
\end{equation}
One may wonder regarding additional differences between the proof in the setting of \cite{AMTn} 
and the one considered here. Essentially, the main tools established here 
are the derivation of certain properties of $\textnormal{dom}$ and $\textnormal{Dom}$ (e.g.., Lemma \ref{lem:Dom3}),  
all of them seem to be new but some of them (see below) generalize known results.  
We feel that when known results are generalized, the approach given here either illuminates 
certain implicit or explicit key points discussed elsewhere or establishes new tools not 
discussed elsewhere. Again, even in the case where a generalization is made, 
because of the general setting considered here there are several difficulties and 
subtle points that should be handled correctly and one example was mentioned 
above regarding \eqref{eq:UnionIntersectBody}. One may also wonder whether it is 
possible to simplify the proof in the specific case of the Euclidean plane (or at least 
for Euclidean spaces). Unfortunately, it seems that the answer is no, partly because (as mentioned 
above) already the Euclidean proof is not simple. But this can be regarded as an advantage 
because it shows that the arguments given here use only certain mild but important properties 
which certain spaces have, among them Euclidean spaces. 

As a final remark regarding this issue, it may be interesting to note that there are 
some connections between some auxiliary tools used here and in  \cite{ImaiKawamuraMatousekReemTokuyamaCGTA},  
e.g., Proposition 5. In addition, the setting of the existence assertion in 
\cite{ImaiKawamuraMatousekReemTokuyamaCGTA} (but not of the computational part, namely 
\cite[Section 4]{ImaiKawamuraMatousekReemTokuyamaCGTA}, in which the setting is a finite dimensional 
strictly convex normed space) is 
 proper geodesic metric spaces. Nevertheless 
 the proofs of \cite[Proposition 5]{ImaiKawamuraMatousekReemTokuyamaCGTA} and the assertions  
 mentioned above (e.g., Theorem \ref{thm:AlgUniConvex}) are different since the involved mappings are 
  not the same  and hence one needs to find a strategy which is appropriate for each case separately. An important feature used in the proof of Theorem \ref{thm:AlgUniConvex} is the fact that the sites are positively separated (this also allows us to consider infinitely many sites) while in  \cite[Proposition ~5]{ImaiKawamuraMatousekReemTokuyamaCGTA} the sites satisfy the more general condition of being merely disjoint, but now one must consider only two sites (by definition) and the mapping must have finitely many components. In addition, it seems that some of the auxiliary 
assertions described in \cite[Section 4]{ImaiKawamuraMatousekReemTokuyamaCGTA} can be 
actually generalized to proper metric spaces having the geodesic inclusion property using  some 
of the tools mentioned in Section  \ref{sec:proofs}. See Section \ref{sec:TopologicalProperties} for more details 
in a specific case. 


\begin{expl}\label{ex:mM}
An illustration of Theorem \ref{thm:AlgUniConvex} for the case where $m\neq M$ is given in Figures~ \ref{fig:GreatestDoubleZoneStrangeNorm}-\ref{fig:LeastDoubleZoneStrangeNorm}. 
Here the sites are $P_1=\{(0,1)\}$, $P_2=\{(0,-1)\}$ in a rectangle in $\mathbb{R}^2$, but the 
norm is $\|(x_1,x_2)\|=\delta\sqrt{\alpha^2|x_1|^2+|x_2|^2}+(1-\delta\alpha)|x_1|+(1-\delta)|x_2|$, where 
$\alpha=\delta=0.1$. This norm is strictly convex but not smooth (the unit sphere contains points which have 
several supporting lines). 
A closely related example was discussed in \cite[Section 5]{KMT2012} in the context of 
non-uniqueness of zone diagrams (double zone diagrams were not mentioned) and it actually follows 
from the discussion presented there that indeed $m\neq M$ in this case. Two different zone diagrams can be obtained 
now: either $(m_1,M_2)$ or $(M_1,m_2)$. The first of them is shown in Figure \ref{fig:ZD-StrangeNorm-2Sites-0002}. 


Figure \ref{fig:ZD-Iteration2-StrangeNorm-2Sites-0002} shows an approximation of $m$ and $M$ using $I^{(2)}$ and $O^{(2)}$.
For producing this figure (see also Section \ref{sec:ActualCompute}) the endpoints of 4000 rays emanating from 
each site were computed in the Voronoi algorithm stage. However, it can be seen that still some parts of the boundaries 
of the components of $I^{(2)}$ and $O^{(2)}$ are not full. 
The reason is that many rays should be produced in a very small angle (i.e., a to belong to the intersection 
of a very narrow cone with the unit sphere) because of the location of the site with respect to these boundary parts (these 
parts and the site are almost located on the same line, thus the rays should be produced in directions 
very close to the two possible directions of the line). A simple way to overcome 
the problem (after detected) is to distribute the rays in the unit sphere in a non-uniform way such 
that many rays will emanate in the corresponding sector but relatively few in other sectors, in contrast to  
the current way of producing the rays in a roughly uniform way. Alternatively (or complementary), 
one can also connect neighbor endpoints by line segments. 
\end{expl}

\begin{figure}
\begin{minipage}[htb]{0.45\textwidth}
\begin{center}
{\includegraphics[scale=0.51]{ZD-Max-StrangeNorm-2Sites-0002.png}}
\end{center}
 \caption{The components of the greatest double zone diagram $M$ from Example ~\ref{ex:mM}.}
\label{fig:GreatestDoubleZoneStrangeNorm}
\end{minipage}
\begin{minipage}[htb]{0.48\textwidth}
\begin{center}
{\includegraphics[scale=0.51]{ZD-Min-StrangeNorm-2Sites-0002.png}
\end{center}
 \caption{The components of the least double zone diagram $m$ from Example ~\ref{ex:mM}.}
\label{fig:LeastDoubleZoneStrangeNorm}
\end{minipage}


\begin{minipage}[t]{0.45\textwidth}
\begin{center}
{\includegraphics[scale=0.51]{ZD-Iteration2-StrangeNorm-2Sites-0002.png}}
\end{center}
 \caption{$I^{(2)}$ and $O^{(2)}$ for the double zone diagrams from Example~\ref{ex:mM}.}
\label{fig:ZD-Iteration2-StrangeNorm-2Sites-0002}
\end{minipage}
\begin{minipage}[t]{0.48\textwidth}
\begin{center}
{\includegraphics[scale=0.49]{ZD-StrangeNorm-2Sites-0002.png}
\end{center}
 \caption{The zone diagram $(m_1,M_2)$ from Example ~\ref{ex:mM}.}
\label{fig:ZD-StrangeNorm-2Sites-0002}
\end{minipage}
\end{figure}



\section{The actual computation of $I^{(n)}$ and $O^{(n)}$, time complexity, and clarifications}\label{sec:ActualCompute}
In this section we provide details about practical and theoretical aspects of the method which allows the 
approximate computation of the corresponding regions. The section has three parts: in the first part practical details 
regarding the method are given. In the second part the issue of time complexity is considered. In the 
last part several related issues are clarified. Only the case where the geodesic 
metric space is actually a compact subset of a normed space is going to be considered, 
because in this case we have a working implementation. In principle the details can be generalized, but, 
however, no attempt has been done so far to implement the algorithm in this setting.  

\subsection{\bf The actual computation of $I^{(n)}$ and $O^{(n)}$:}
As mentioned in Section \ref{sec:QualitativeDescription}, in order to compute the iterations $I^{(n)}$ and $O^{(n)}$ one has to know how to compute or approximate dominance regions $\textnormal{dom}(P,A)$ of general sets $P,A$. The way we choose to overcome this difficulty is to use the approximation algorithm for computing Voronoi diagrams of general sites in general normed spaces which was introduced in \cite{ReemISVD09,ReemPhD}. In a nutshell, for approximating $\textnormal{dom}(P,A)$ one uses the fact that $\textnormal{dom}(P,A)$ can be represented as a union of  rays emanating from the points of $P$ (see Figures~  \ref{fig:Zone80}-\ref{fig:Zone160} and  \ref{fig:VoronoiL2718281828}-\ref{fig:ZoneL2718281828} for an illustration). 

First, the world $X$ is assumed to be a large compact subset, e.g., a square or a  hypercube. Now one approximates $P$ using a finite collection of points; this is always possible if  $P$ is compact. After that a  finite collection of directions (unit vectors) is chosen,  then one approximates the endpoints of the line segments emanating from the points of $P$ in these directions. These endpoints can be approximated to any required precision, and then they are stored. At this point $\textnormal{dom}(P,A)$ is represented by this collection of endpoints and it is regarded as computed. 

For computing $I^{(n)}$ and $O^{(n)}$ one  computes their corresponding components iteratively: the components are $\textnormal{dom}(P_k,A_k)$ where $k$ runs over all the (finitely many in practice) indices in $K$ and $A_k$ depends on the iteration and on $k$. In general, $A_k=\bigcup_{j\neq k}(I^{(n)})_k$ or $A_k=\bigcup_{j\neq k}(O^{(n)})_k$ where $(I^{(n)})_k$ is the $k$-th component of $I^{(n)}$. For instance, the components of $I^{(1)}$ are $\textnormal{dom}(P_k,A_k)$ where $A_k=\bigcup_{j\neq k} (O^{(0)})_k$ and $(O^{(0)})_k=\textnormal{dom}(P_k,\bigcup_{j\neq k}P_j)$ is known, as a  collection of endpoints,  from a previous computation.   An illustration of the process is given in Figures  \ref{fig:Zone2},\ref{fig:Zone3}, \ref{fig:Zone80},\ref{fig:Zone160},\ref{fig:ZoneLinfty}, and \ref{fig:ZoneL1}.  

\begin{figure}
\begin{minipage}[t]{0.45\textwidth}
\begin{center}
{\includegraphics[scale=0.6]{Voronoi2D-4Site-3inSite-01-L2718281828.png}}
\end{center}
 \caption{The Voronoi diagram of 4 sites in $(\mathbb{R}^2,\ell_p)$, where $p\approx 2.71$; 
 each consists of 3 points; 80 rays  emanate from each point of each site.}
\label{fig:VoronoiL2718281828}
\end{minipage}
\hfill
\begin{minipage}[t]{0.45\textwidth}
\begin{center}
{\includegraphics[scale=0.6]{Zone2D-4Site-3inSite-5Iter-01-L2718281828.png}}
\end{center}
 \caption{The zone diagram (and hence the double zone diagram) of the sites of Figure \ref{fig:VoronoiL2718281828}.}
\label{fig:ZoneL2718281828}
\end{minipage}
\end{figure}


An important parameter which determines the level of approximation is the number of rays (or, actually, the density of the corresponding unit vectors in the unit sphere of the space) used in the construction of the dominance regions. In Figures \ref{fig:Zone2} and \ref{fig:Zone3} about 1600 rays  emanate from each site, and only the endpoints are displayed. In Figure \ref{fig:ZD-Iteration2-StrangeNorm-2Sites-0002} about  
4000 rays emanate from each site and only endpoints are displayed. 
In Figures \ref{fig:Zone80},\ref{fig:VoronoiL2718281828},\ref{fig:ZoneL2718281828} eighty rays have been used in each iteration, and the whole ray is displayed. 160 rays (emanating from each site in each iteration) 
have been used to produce Figure \ref{fig:Zone160}. Despite the difference in the level of 
approximation, all of the  approximations of the zone diagram of Figure \ref{fig:IntroZone} 
(namely, Figures \ref{fig:IntroZone}, \ref{fig:Zone80}, and \ref{fig:Zone160}) give roughly the same shapes. 



\subsection{\bf Time complexity:}\label{sec:complexity}
As for the time complexity of the algorithm, it is determined by the time complexity of $\textnormal{dom}(P_k,A_k)$ for each $k$. 
Starting from iteration $n=1$, the number of endpoints in $A_k$ is $C_1 C_2^{d-1}G(|K|-1)$ where  $C_1$ is a universal positive constant depending on the norm, $C_2$ is a positive constant depending on the level of approximation of the unit sphere of the space 
(roughly speaking, it is a bound on the distance between a point on the unit sphere and some unit vector from the list of unit vectors created by the user), $d$ is the dimension, $G$ is a bound on the number of points in $P_k$ (for each $k$), $C_1C_2^{d-1}G$ is the number of rays emanating from the points of $P_k$, and $|K|$ is the number of sites. For determining an endpoint up to some error parameter one makes $C_3 |A_k|$ comparisons where $C_3$ is a positive constant depending on the level of approximation. Hence the time complexity for the entire number of components (for a given iteration) is 
$O(C_3(C_1C_2^{d-1} G |K|)^2)$. Note that here $G$ and $|K|$ are parameters depending on the input, $C_1$ and $d$ are global parameters, and $C_2$ and $C_3$ are parameters depending on the user (the level of approximation of the computation). Under some assumptions, e.g., that 
the sites are points which are uniformly distributed, it is possible to reduce significantly the number of calculations.  


\subsection{\bf Clarifications:}
We finish this section with a discussion on certain issues which perhaps have not been so clear so far. 
It should be emphasized that the goal of this section  is to  describe schematically a 
practical way for  approximating the regions which appear in each iteration and to evaluate roughly the number 
of calculations in each iteration. The fact that the algorithm can approximate a given dominance region up to any desired 
precision is a consequence of the stability of the algorithm presented in \cite{ReemISVD09}. 
However, the proof of this property is not obvious and it will be discussed in a paper 
devoted to this algorithm. The proof is in the spirit of the proof of the geometric stability 
of Voronoi diagrams with respect to small changes of the sites \cite{ReemGeometricStabilityArxiv}. 

Roughly speaking, if the users want to approximate (in terms 
of the Hausdorff distance) a given region up to some error parameter $\epsilon$, then they need to approximate 
well enough the sites, to choose in advance enough approximating rays (this is determined by the error parameter related 
to the unit sphere), and to fix a small enough error parameter for the endpoints of the rays. If the users decide 
in advance how many iterations they want to perform for approximating the (double) zone diagram and what is 
the level of approximation of the regions in the final iteration, then they can estimate in advance the number 
of calculations by iterative ``reverse engineering'': using the target error parameter $\epsilon$, one estimates 
the error parameters needed as input for the final iteration, and from them the error parameters needed for 
the previous iteration, and so on, until the initial iteration. This gives an estimate on the initial error parameters. 

In the above description the number of iterations $n$ was chosen by the user but there was no guarantee that the real 
(double) zone diagram will be approximated well by $I^{(n)}$ or $O^{(n)}$. However, since the algorithm 
converges to the double zone diagram by Theorem \ref{thm:AlgUniConvex}, then given $\epsilon>0$, there is a number $n_0$, depending only on $\epsilon$, such that for any $n\geq n_0$, the regions of any iteration $I^{(n)}$ and $O^{(n)}$ will be at Hausdorff  distance of at most $\epsilon$ from the limit regions. Hence, if one can compute $n_0(\epsilon)$, 
then one can know in advance the corresponding needed initial error parameters. Unfortunately, it is not clear how 
to estimate $n_0(\epsilon)$ and this is a major open problem.


\section{A byproduct: Topological properties of Voronoi cells}\label{sec:TopologicalProperties}
It may be of interest to give more details regarding a certain interesting 
auxiliary tool which is established along the way. Roughly speaking, it says 
that the corresponding Voronoi cells of the given sites can be written as 
the closure of their strict inequality part (the set of all $x$ in the cells strictly 
closer to their site than to the other sites). This is closely related to the fact 
that the boundaries of the cells must coincide with their bisectors (no ``fat'' bisectors can occur). 
Both properties are formulated below, namely Proposition ~\ref{prop:StrictInt} and Proposition ~\ref{prop:BoundaryInterior}. 
The property mentioned in Proposition ~\ref{prop:StrictInt} was not used in \cite{AMTn} but it was used in the 
related paper \cite[Lemma 6]{ImaiKawamuraMatousekReemTokuyamaCGTA} for proving a related convergence 
theorem \cite[Proposition 5]{ImaiKawamuraMatousekReemTokuyamaCGTA}. 
Proposition ~\ref{prop:StrictInt} 
actually generalizes \cite[Lemma 6]{ImaiKawamuraMatousekReemTokuyamaCGTA} 
(whose setting is finite dimensional strictly convex normed spaces) 
and its proof is inspired from a somewhat different  proof of 
\cite[Lemma 6]{ImaiKawamuraMatousekReemTokuyamaCGTA} rather 
than directly from \cite[Lemma 6]{ImaiKawamuraMatousekReemTokuyamaCGTA}. A careful analysis 
of this second proof enabled us to generalize it to geodesic metric spaces having the geodesic inclusion property 
and actually to discover this class of spaces. 

\begin{prop}\label{prop:StrictInt}
Let $(X,d)$ be a geodesic metric space which has the geodesic inclusion 
property. Let $P,A\subseteq X$ be nonempty. Suppose that  $P\bigcap A=\emptyset$ 
and that the distances $d(x,P)$ and $d(x,A)$ are attained for all $x\in X$. Then 
\begin{equation}
\begin{array}{lll}\label{eq:closure}
\textnormal{dom}(P,A)&=&\overline{\{x\in X: d(x,P)<d(x,A)\}}.
\end{array}
\end{equation}
 In particular, the  above conclusion holds when all the above assumptions hold 
 with the exception of the distance attainment one 
but instead $(X,d)$ is proper and both $P$ and $A$ are closed.  
\end{prop}
\begin{prop}\label{prop:BoundaryInterior}
Under the assumptions of Proposition ~\ref{prop:StrictInt} the following hold: 
\begin{equation}
\partial(\textnormal{dom}(P,A))=\{x\in X: d(x,P)=d(x,A)\},\label{eq:boundary}
\end{equation}
\begin{equation}
\textnormal{int}(\textnormal{dom}(P,A))=\{x\in X: d(x,P)<d(x,A)\},\label{eq:interior} 
\end{equation}
\end{prop}
The properties mentioned in 
Proposition ~\ref{prop:BoundaryInterior} have been essentially known for a long time before their present 
formulation (or before \cite[Lemma 6]{ImaiKawamuraMatousekReemTokuyamaCGTA}) 
at least in the case of point sites in certain finite dimensional spaces. 
But finding them in the formulation given in Proposition ~\ref{prop:BoundaryInterior} 
does not seem to be an easy task (no such a formulation has been found). 

In general, understanding bisectors, especially in high dimensions, 
is not a well explored area. But see e.g., 
 \cite{CorbalanMazonRecio1996,Horvath2000,HorvathMartini2013,
IKLM1995,Le1997,LevenSharir1987,Ma2000,MartiniSwanepoel2004} for several works 
discussing somewhat related properties of bisectors (either topological or combinatorial) 
in certain settings. Recently results similar to 
  Proposition ~\ref{prop:BoundaryInterior} have been established in two cases. 
The first is in \cite{ReemVorStabilityNonUC2012}, where the setting is a general 
(possibly infinite dimensional) normed space assuming the distance to the sites is 
attained and assuming that the sites are aligned in a certain way with respect to the 
structure of the unit sphere of the space. The second 
case is for arbitrary sites satisfying $d(P,A)>0$ and located 
in infinite dimensional uniformly convex normed spaces \cite{ReemTopologicalPropertiesPreprint2013}.  

As a final remark, it should be emphasized that there are well known counterexamples 
to \eqref{eq:closure}-\eqref{eq:interior} 
even in the case of $(\mathbb{R}^2,\ell_{1})$ with two point sites. 
See, e.g., \cite[p. 390, Figure 37]{Aurenhammer}, \cite[p. 191, Figure 3.7.2]{OBSC}. 
For additional counterexamples (of a somewhat different nature), 
see  \cite{ReemGeometricStabilityArxiv,ReemVorStabilityNonUC2012,ReemTopologicalPropertiesPreprint2013}. 

\section{Open problems}\label{sec:end}

\begin{figure}
\begin{minipage}[t]{0.42\textwidth}
\begin{center}
{\includegraphics[scale=0.55]{ZoneLinfty0005.png}}
\end{center}
 \caption{$I^{(1)},\ldots, I^{(4)}$ and $O^{(1)},\ldots, O^{(4)}$ for the sites  $\{(0,3)\},\,\{(0,-3)\}$ in a square in $(\mathbb{R}^2,\ell_{\infty})$; the bold line is $I^{(3)}\approx I^{(4)}\approx O^{(3)}\approx O^{(4)}$. See also \cite[Example 2.4]{ReemReichZone}. }
\label{fig:ZoneLinfty}
\end{minipage}
\hfill
\begin{minipage}[t]{0.4\textwidth}
\begin{center}
{\includegraphics[scale=0.56]{Zone2D-I2-O2-Site30-0025-L1.png}}
\end{center}
 \caption{$I^{(2)}$ and $O^{(2)}$ for 30 point sites in a square in $(\mathbb{R}^2,\ell_1)$.}
\label{fig:ZoneL1}
\end{minipage}
\end{figure}

The non-proof part of the paper is concluded with several interesting open problems and unexplained phenomena.
First, perhaps the most interesting open problem is to establish error estimates for  the convergence speed of $(I^{(n)})_{n=0}^{\infty}$ and  $(O^{(n)})_{n=0}^{\infty}$. One of the reasons that this problem is interesting, is the observed fast convergence of these sequences (usually 4  iterations suffice). We feel that the approach and error estimates developed in  \cite{ReemGeometricStabilityArxiv} may help here. Second, although in general a zone diagram is not necessarily  unique, even in the case of $\mathbb{R}^2$ with two point sites (see e.g., \cite[Sections 1,5]{KMT2012} or Example \ref{ex:mM}), and although our computational results are limited to strictly convex normed spaces, experiments show  that the algorithm actually converges to a unique zone diagram most of the times also in the case of  non-strictly convex normed spaces; 
see Figures \ref{fig:ZoneLinfty}-\ref{fig:ZoneL1} for an illustration. 


Third, it is interesting to find out whether Theorem \ref{thm:AlgUniConvex} can be generalized to other settings, for instance to all proper geodesic metric spaces, or at least to all finite dimensional 
normed spaces (even ones which are not strictly convex) and to some ``nice'' manifolds. 
Removing the properness assumption  is of interest too, even in the 
particular case of infinite dimensional normed spaces. 
Proving Theorem \ref{thm:AlgUniConvex} in the full generality of any metric space is, of course, 
a much stronger statement than what described above, but it seems out of reach now. In fact, 
it is not even clear that this is true (we do not have any counterexample, but perhaps 
one can construct counterexamples in a discrete setting; perhaps eventually normed spaces which are 
not strictly convex will provide counterexamples; this is true at least partially in the sense that
 there are simple counterexamples to \eqref{eq:AntiCont} even in $(\mathbb{R}^2,\ell_{\infty})$: just take $P=\{(0,-3)\}$ and $C_{\gamma}=[-1/\gamma,0]\times\{(0,3)\}$). On the other hand, since 
$\textnormal{Dom}(m)=M$ holds in general, this may imply that also $m=\textnormal{Dom}(M)$ 
holds in general, at least for positively separated sites. 


\section{Proof of the convergence theorems and related claims}\label{sec:proofs}
This section presents the proofs of the convergence theorems and related claims. 

The first lemmas describe simple properties of dom and Dom. 
 
\begin{lem}\label{lem:GeneralDistance} Let $(X,d)$ be a metric space. Then
\begin{enumerate}[(a)]



\item\label{item:d_x_Agamma} $d(x,\bigcup_{\gamma\in \Gamma}A_{\gamma})=\inf\{d(x,A_{\gamma}): \gamma\in \Gamma\}$ for any $x\in X$ and any collection $\{A_{\gamma}\}_{\gamma\in \Gamma}$ of nonempty subsets in $X$.
\item\label{item:UnionIntersection} $\textnormal{dom}(P,\bigcup_{\gamma\in \Gamma} A_{\gamma})=\bigcap_{\gamma\in \Gamma} \textnormal{dom}(P,A_{\gamma})$
for any collection $\{A_{\gamma}\}_{\gamma\in \Gamma}$ of nonempty subsets in $X$ and any $P\subseteq X$ nonempty.
\item\label{item:FirstComponent} $\textnormal{dom}(\bigcup_{\gamma\in \Gamma} P_{\gamma},A)\subseteq\bigcup_{\gamma\in \Gamma} \textnormal{dom}(P_{\gamma},A)$
for any collection $\{P_{\gamma}\}_{\gamma\in \Gamma}$ of nonempty subsets in $X$ and any $A\subseteq X$ nonempty. If, in addition, for each $x\in X$ the distance $d(x,\bigcup_{\gamma\in \Gamma}P_{\gamma})$ is attained at some $P_{\gamma_0}$, i.e., $d(x,\bigcup_{\gamma\in \Gamma}P_{\gamma})=d(x,P_{\gamma_0})$ for some $\gamma_0\in \Gamma$, then the equality $\textnormal{dom}(\bigcup_{\gamma\in \Gamma} P_{\gamma},A)=\bigcup_{\gamma\in \Gamma} \textnormal{dom}(P_{\gamma},A)$ holds.
\item \label{item:DomUnionIntersect}$\textnormal{Dom}(\bigcup_{\gamma\in \Gamma} R^{\gamma})=\bigcap_{\gamma\in\Gamma} \textnormal{Dom}(R^{\gamma})$ for any collection $\{R^{\gamma}\}_{\gamma\in \Gamma}$
of tuples, where each $R^{\gamma}$ is indexed by the same set of
indices $K$, i.e., $R^{\gamma}=(R_k^{\gamma})_{k\in K}$. Here $\textnormal{Dom}$ is defined with respect  to some tuple $(P_k)_{k \in K}$ of nonempty subsets in $X$.


\end{enumerate}
\end{lem}

\begin{proof}
\begin{enumerate}[(a)]

\item Let $x\in X$, $\alpha:=d(x,\bigcup_{\gamma\in \Gamma}A_{\gamma})$ and $\beta:=\inf\{d(x,A_{\gamma}): \gamma\in \Gamma\}$. Then $\alpha\leq d(x,A_{\gamma})$ for all $\gamma\in \Gamma$ by the definition of $\alpha$, 
so $\alpha\leq \beta$.  If $\alpha<\beta$, then there is $y\in \bigcup_{\gamma\in \Gamma}A_{\gamma}$ such that $d(x,y)<\beta$, and since $y\in A_{\gamma}$ for some $\gamma\in \Gamma$, we have $d(x,A_{\gamma})\leq d(x,y)<\beta$, a contradiction with the definition of $\beta$.
\item By the antimonotonicity of $\textnormal{dom}(P,\cdot)$ we have $\textnormal{dom}(P,\bigcup_{i\in \Gamma} A_{i})\subseteq
\textnormal{dom}(P,A_{\gamma})$ for any ${\gamma\in \Gamma}$. Consequently, we have 
$\textnormal{dom}(P,\bigcup_{\gamma\in \Gamma} A_{\gamma})\subseteq
\bigcap_{\gamma\in \Gamma}\textnormal{dom}(P,A_{\gamma})$. Conversely, 
suppose that $x\in
\bigcap_{\gamma\in \Gamma}\textnormal{dom}(P,A_{\gamma})$. If $x\notin
\textnormal{dom}(P,\bigcup_{\gamma\in \Gamma} A_{\gamma})$, then there is $y\in
\bigcup_{\gamma\in \Gamma} A_{\gamma}$ such that $d(x,P)>d(x,y)$. But $y\in
A_{\gamma}$ for some ${\gamma}$, and $x\in
\textnormal{dom}(P,A_{\gamma})$, so $d(x,P)\leq d(x,A_{\gamma})\leq
d(x,y)$, a contradiction.
\item By monotonicity of $\textnormal{dom}(\cdot, A)$, i.e., $\textnormal{dom}(P,A)\leq \textnormal{dom}(Q,A)$ whenever $P\subseteq Q$ (which is easily verified), 
we have $\textnormal{dom}(P_{\gamma},A)\subseteq \textnormal{dom}(\bigcup_{i\in \Gamma} P_{i},A)$ for any ${\gamma}$. As a result, we have 
$\bigcup_{\gamma\in \Gamma}\textnormal{dom}(P_{\gamma},A)\subseteq \textnormal{dom}(\bigcup_{\gamma\in \Gamma}
P_{\gamma},A)$. Now suppose that $d(x,\bigcup_{\gamma\in \Gamma}P_{\gamma})=d(x,P_{\gamma_0})$ for some $\gamma_0\in \Gamma$. Therefore,  if $x\in\textnormal{dom}(\bigcup_{\gamma\in \Gamma} P_{\gamma},A)$, then $d(x,P_{\gamma_0})=d(x,\bigcup_{\gamma\in \Gamma}P_{\gamma})\leq d(x,A)$, i.e., $x\in  \textnormal{dom}(P_{\gamma_0},A)\subseteq\bigcup_{\gamma\in \Gamma}\textnormal{dom}(P_{\gamma},A)$.
\item Since the union and intersection are taken component-wise, this follows from the   definition of the $\textnormal{Dom}$ mapping (equation \eqref{eq:TZoneDef})
and part (\ref{item:UnionIntersection}):
\begin{multline*}
\quad \textnormal{Dom}(\bigcup_{\gamma\in
\Gamma}R^{\gamma})=(\textnormal{dom}(P_k,\bigcup_{j\neq
k}(\bigcup_{\gamma\in \Gamma}{R_j^{\gamma}})))_{k\in
K}=(\textnormal{dom}(P_k,\bigcup_{\gamma\in \Gamma}(\bigcup_{j\neq
k}{R_j^{\gamma}})))_{k\in
K}\\=(\bigcap_{\gamma\in \Gamma}\textnormal{dom}(P_k,\bigcup_{j\neq
k}{R_j^{\gamma}}))_{k\in K}=\bigcap_{\gamma\in \Gamma}\textnormal{Dom}(R^{\gamma}).
\end{multline*}
\end{enumerate}
\end{proof}
 
\begin{lem}\label{lem:DomPk} 
\begin{enumerate}[(a)
Let $(P_k)_{k \in K}$ be a tuple of  nonempty and closed subsets in a metric space $(X,d)$. 
\item $\textnormal{Dom}$ is antimonotone, i.e., $\textnormal{Dom}(R)\subseteq \textnormal{Dom}(S)$ whenever $S\subseteq R$; $\textnormal{Dom}^2$ is monotone, that is, $R\subseteq S\Rightarrow \textnormal{Dom}^2(R)\subseteq \textnormal{Dom}^2(S)$. 
 
\item  The sequence $(I^{(n)})_{n=0}^{\infty}$ is increasing, $(O^{(n)})_{n=0}^{\infty}$ is decreasing 
and we have
\begin{equation}\label{eq:InOn}
I^{(n)}\subseteq \bigcup_{\gamma=0}^{\infty}\textnormal{Dom}^{2\gamma}(P_k)_{k\in K}\subseteq \bigcap_{\gamma=0}^{\infty}\textnormal{Dom}^{2\gamma+1}(P_k)_{k\in K}\subseteq O^{(n)}
\end{equation}
for each nonnegative $n$. In addition, if $R$ is a zone or double zone diagram in $X$, then
\begin{equation}\label{eq:DomIterations}
I^{(n)}\subseteq\bigcup_{\gamma=0}^{\infty}\textnormal{Dom}^{2\gamma}(P_k)_{k\in K}\subseteq R \subseteq \bigcap_{\gamma=0}^{\infty}\textnormal{Dom}^{2\gamma+1}(P_k)_{k\in K}\subseteq O^{(n)}\quad \forall n\geq 0.
\end{equation}
\end{enumerate}
\end{lem}
\begin{proof}
\begin{enumerate}[(a)
\item The antimonotonicity of $\textnormal{Dom}$ is a simple consequence of the easily verified  fact that $\textnormal{dom}(P,B)\subseteq \textnormal{dom}(P,A)$ for any nonempty subset $P$ and any nonempty subsets $A,B$ satisfying  $A\subseteq B$. The composition of two antimonotone mapping is monotone, and hence $\textnormal{Dom}^2$ is monotone. 
\item Let $(X)_{k\in K}$ be the tuple whose components are $X$. An immediate check shows that $P_k=\textnormal{dom}(P_k,X)$ and $\textnormal{Dom}(X)_{k\in K}=(P_k)_{k \in K}$. From the  inclusion $(P_k)_{k\in K}\subseteq (X)_{k\in K}$, the definition of the space on which $\textnormal{Dom}$ acts, and the monotonicity of $\textnormal{Dom}^2$ we have 
\begin{equation*}
(P_k)_{k\in K}\subseteq \textnormal{Dom}^2 (P_k)_{k\in K}\subseteq \textnormal{Dom}^2 (X)_{k\in K}=\textnormal{Dom}(P_k)_{k\in K}\subseteq (X)_{k\in K}.
\end{equation*}
As a result, by iterating with $\textnormal{Dom}^2$, using its monotonicity, and using the definition of $I^{(n)}$ and $O^{(n)}$ we see that $(I^{(n)})_{n=0}^{\infty}$ is increasing and $(O^{(n)})_{n=0}^{\infty}$ is decreasing. Hence $I^{(n)}$ is contained in $\bigcup_{\ell=0}I^{(\ell)}$ and $\bigcap_{\ell=0}O^{(\ell)}$ is contained in $O^{(n)}$ for each $n$. 
Finally, after iterating the inclusion $(P_k)_{k\in K}\subseteq R\subseteq (X)_{k\in K}$ using $\textnormal{Dom}^2$ we obtain  \eqref{eq:DomIterations} when $R$ is a double zone diagram. But a zone diagram is a double  zone diagram 
[$\textnormal{Dom}(\textnormal{Dom}(R))=\textnormal{Dom}(R)=R$ whenever $R=\textnormal{Dom}(R)$], so  \eqref{eq:DomIterations} is true in this case too.
\end{enumerate}
\end{proof}

\begin{lem}\label{lem:dom}
Let $(X,d)$ be a metric space and let $P,A\subseteq X$ be nonempty.  Then 
\begin{enumerate}[(a)]
\item\label{item:closure} $\textnormal{dom}(P,A)$ is a closed set.
\item\label{item:Aclosure} $\textnormal{dom}(P,A)=\textnormal{dom}(P,\overline{A})=\textnormal{dom}(\overline{P},A)=\textnormal{dom}(\overline{P},\overline{A})$.
\item\label{item:DomClosure} For a tuple $R=(R_k)_{k\in K}$ of nonempty subsets, let $\overline{R}=(\overline{R_k})_{k\in K}$.
Then $\textnormal{Dom}(\overline{R})=\textnormal{Dom}(R)=\overline{\textnormal{Dom}(R)}$ (with respect to a given  tuple $(P_k)_{k\in K}$).
\end{enumerate}
\end{lem}
\begin{proof}
\begin{enumerate}[(a)]
\item If $x=\lim_{n\to \infty}x_n$ and $x_n\in \textnormal{dom}(P,A)$, then $d(x_n,P)\leq
d(x_n,A)$, and this inequality is preserved in the limit because the
function $x\mapsto d(x,P)-d(x,A)$ is continuous (with respect to the topology 
induced on $X$ by $d$).
\item This follows from the general fact that $d(x,B)=d(x,\overline{B})$ for each  $x\in X$ and each $B\subseteq X$.
\item By part \eqref{item:Aclosure} and by Lemma \ref{lem:GeneralDistance}\eqref{item:UnionIntersection}, we have
\begin{equation*}
\quad\quad \textnormal{dom}(P_k,\bigcup_{j\neq k}\overline{R_j})=\bigcap_{j\neq k}\textnormal{dom}(P_k,\overline{R_j})
=\bigcap_{j\neq k}\textnormal{dom}(P_k,R_j)=\textnormal{dom}(P_k,\bigcup_{j\neq k}R_j).
\end{equation*}
Thus $\textnormal{Dom}(\overline{R})=(\textnormal{dom}(P_k,\bigcup_{j\neq k}\overline{R_j}))_{k\in K}=(\textnormal{dom}(P_k,\bigcup_{j\neq k}R_j))_{k\in K}=\textnormal{Dom}(R)$. Finally, by part \eqref{item:closure} and by definition, we have
\begin{equation*}
\overline{\textnormal{Dom}(R)}=(\overline{\textnormal{dom}(P_k,\bigcup_{j\neq k}P_j)})_{k\in K}=
(\textnormal{dom}(P_k,\bigcup_{j\neq k}P_j))_{k\in K}=\textnormal{Dom}(R).
\end{equation*}
\end{enumerate}
\end{proof}


The following propositions and lemma (followed by some proofs 
of related claims) discuss issues related to geodesic metric spaces. 


\begin{prop}\label{prop:GeodesicInclusionStrictlyConvex}
Any strictly convex normed space $(X,d)$ (or any nontrivial convex subset of such a space) 
has the geodesic inclusion property. 
\end{prop}
\begin{proof}
To see this, first note that a geodesic segment 
$[x,y]_{\gamma}$ connecting $x=\gamma(r_1)$ and $y=\gamma(r_2)$ (where $r_2> r_1$ 
and $\gamma:[r_1,r_2]\to X$ preserves distances) must be the line segment $[x,y]$.  
Indeed let $b=\gamma(r)\in [x,y]_{\gamma}$  for some $r\in (r_1,r_2)$. Then 
\begin{multline*}
d(x,y)\leq d(x,b)+d(b,y)=d(\gamma(r_1),\gamma(r))+d(\gamma(r),\gamma(r_2))\\
=r-r_1+r_2-r=r_2-r_1=d(x,y), 
\end{multline*}
showing that $d(x,b)+d(b,y)=d(x,y)$. This is in contrast to the assumed strict convexity of the space 
unless  $b$ is in the line segment $[x,y]$. Thus $[x,y]_{\gamma}\subseteq [x,y]$. 
On the other hand, since $\gamma$ is continuous, its image $\gamma([r_1,r_2])=[x,y]_{\gamma}$ 
is a (nontrivial) connected subset of $[x,y]$. But the only nontrivial connected subsets of a line segment 
are line segments. We conclude that $[x,y]_{\gamma}$ is a line segment containing 
$x$ and $y$ and the equality $[x,y]_{\gamma}=[x,y]$ follows.

Returning to the geodesic inclusion property,  given $u,v,b,z\in X$ satisfying 
$b\neq z$, $b\in [u,z]_{\gamma_1}$, and $b\in [v,z]_{\gamma_2}$, we already know from 
the above that the geodesic segments are linear ones. The property obviously holds if $u=b$ or $v=b$. 
Otherwise  $(b-z)/d(b,z)=(u-b)/d(u,b)$ and $(b-z)/d(b,z)=(v-b)/d(v,b)$, 
showing that $[b,u]$ and $[b,v]$ are on the same ray emanating from $b$ in the joint 
direction $(b-z)/d(b,z)$. Hence the shorter interval among $[b,u]$, $[b,v]$ 
is contained in the longer one and this proves the assertion.  
\end{proof}

\begin{lem}\label{lem:GeodesicEquality}
Let $(X,d)$ be a geodesic metric space and let $a,b,c\in X$. 
Then equality in the triangle inequality $d(a,c)\leq d(a,b)+d(b,c)$ holds 
if and only if  $b\in [a,c]_{\gamma}$ for some isometric mapping $\gamma$. 
\end{lem}
\begin{proof}
Suppose first that $b\in [a,c]_{\gamma}$ for some $b\in X$ and some 
isometric mapping $\gamma: [r_1,r_3]\to X$. We can write $a=\gamma(r_1),b=\gamma(r_2)$, and $c=\gamma(r_3)$ 
for some real number $r_2\in [r_1,r_3]$. Since $\gamma:[r_1,r_3]\to X$ preserves 
distances we have 
\begin{multline*}
d(a,b)+d(b,c)=d(\gamma(r_1),\gamma(r_2))+d(\gamma(r_2),\gamma(r_3))=
(r_2-r_1)+(r_3-r_2)=r_3-r_1=d(a,c),
\end{multline*}
namely equality in the triangle inequality. 

Assume now the reverse direction, i.e., that $b\in X$ satisfies $d(a,b)+d(b,c)=d(a,c)$. 
Since $(X,d)$ is a geodesic metric space we can find isometric mappings  $\gamma_1:[r_1,r_2]\to X$ 
and $\gamma_2:[r_2,r_3]\to X$ such that $\gamma_1(r_1)=a$, $\gamma_1(r_2)=b$, $d(a,b)=r_2-r_1$, 
 $\gamma_2(r_2)=b$, $\gamma_2(r_3)=c$, and $d(b,c)=r_3-r_2$ (for some real numbers $r_1\leq r_2\leq r_3$). 
Define $\gamma:[r_1,r_3]\to X$ by 
\begin{equation*}
\gamma(r)=\left\{ \begin{array}{ll}
\gamma_1(r) & \textnormal{if}\,\, r\in [r_1,r_2],\\ 
\gamma_2(r) & \textnormal{if}\,\, r\in [r_2,r_3].
\end{array}
\right. 
\end{equation*}
It will be proved that $\gamma$ is an isometric mapping and that $b\in [a,c]_{\gamma}$. 
By the definition of $\gamma$ we have $\gamma(r_2)=b$ and hence $b\in \gamma([r_1,r_3])$. 
For proving that $\gamma$ preserves 
distances, one observes that when restricted to $[r_1,r_2]$ or to $[r_2,r_3]$ this holds 
because $\gamma_i,i=1,2$ preserve distances. It remains to show that if 
 $t_1\in [r_1,r_2]$ and $t_2\in [r_2,r_3]$, then the equality $d(\gamma(t_1),\gamma(t_2))=t_2-t_1$ holds. 
 Indeed, by the triangle inequality, by the fact that the restriction of $\gamma$ 
 to the intervals $[r_1,r_2]$ and $[r_2,r_3]$ is an isometric mapping,  
and using the properties of $a,b,c$ that were assumed above, it follows that 
\begin{multline*}
d(a,c)\leq d(a,\gamma(t_1))+d(\gamma(t_1),\gamma(t_2))+d(\gamma(t_2),c)\\ 
\leq d(\gamma(r_1),\gamma(t_1))+d(\gamma(t_1), \gamma(r_2))+d(\gamma(r_2), \gamma(t_2))+d(\gamma(t_2),\gamma(r_3))\\
=(t_1-r_1)+(r_2-t_1)+(t_2-r_2)+(r_3-t_2)=(r_2-r_1)+(r_3-r_2)\\
=d(a,b)+d(b,c)=d(a,c).
\end{multline*}
Thus there is equality all over the way. This and 
the fact that $\gamma$ preserves distances when restricted to the intervals $[r_1,r_2]$ and $[r_2,r_3]$ imply  
\begin{equation*} 
d(\gamma(t_1),\gamma(t_2))=d(\gamma(t_1), \gamma(r_2))+d(\gamma(r_2), \gamma(t_2))=(r_2-t_1)+(t_2-r_2)=t_2-t_1, 
\end{equation*}
as required. 
\end{proof}

 
\begin{prop}\label{prop:StrictSegment}
Let $(X,d)$ be a geodesic metric space which has the geodesic 
inclusion property.  Let $A\subseteq X$ be nonempty. 
Let $p,z\in X$ and suppose that $d(z,p)\leq d(z,A)$ and $p\notin A$. 
Suppose also that $d(x,A)$ is attained for all $x\in [p,z)$. 
Then $d(x,p)<d(x,A)$ for all $x\in [p,z)$. In particular, the 
 above conclusion holds when all the above assumptions hold 
 with the exception of the distance attainment one 
but instead $(X,d)$ is proper and $A$ is closed. 
\end{prop}
\begin{proof}
It can be assumed that  $z\neq p$, because otherwise the assertion is obvious (void).
Fix $x\in [p,z)$ and let $a_x\in A$ be chosen such
that $d(x,A)=d(x,a_x)$. 
Since $p\notin A$, $z\neq p$, and $d(z,p)\leq d(z,A)$, it follows that $z,p$ and $a_x$ are all different. 
Using the facts that $[p,z]$ is a geodesic segment, that $x\in [p,z]$, and that $d(z,p)\leq d(z,A)$, it follows 
from the triangle inequality and Lemma \ref{lem:GeodesicEquality} that 
\begin{equation}\label{eq:zxpA}
\quad \quad d(z,x)+d(x,p)=d(z,p)\leq d(z,A)
\leq d(z,a_x)\leq d(z,x)+d(x,a_x),
\end{equation}
so $d(x,p)\leq d(x,a_x)=d(x,A)$. Assume for a contradiction that 
\begin{equation}\label{eq:d(x,p)d(x,a_x)}
d(x,p)=d(x,a_x). 
\end{equation}
Then there is equality in \eqref{eq:zxpA} and hence, because of Lemma \ref{lem:GeodesicEquality} 
it follows that $x\in [a_x,z]_{\gamma_1}$ for some isometric mapping $\gamma_1$. 
Because $x\neq z$ the inclusion $x\in [a_x,z)_{\gamma_1}$ holds.  
However, by assumption $x\in [p,z)_{\gamma_2}$ for another isometric mapping $\gamma_2$. 
Thus, using the geodesic inclusion property of the space 
(where $x$ is the bifurcation point), either $p\in [a_x,z]_{\gamma_1}$ 
or $a_x\in [p,z]_{\gamma_2}$. As explained after the definition of 
the geodesic inclusion property, in the first case actually $p\in[a_x,x]_{\gamma_1}$ (otherwise 
$p\in (x,z]_{\gamma_1}$ and hence, since $\gamma_1$ preserves distances, $d(p,z)<d(x,z)$;  
but $d(x,z)\leq d(p,z)$ because $x\in [p,z]_{\gamma_2}$, a contradiction). Consequently, $p\in[a_x,x]_{\gamma_1}$ and $p\neq a_x$ imply that 
$d(x,a_x)=d(x,p)+d(p,a_x)>d(x,p)$, a contradiction to \eqref{eq:d(x,p)d(x,a_x)}.
Thus the second case implied by the geodesic inclusion property holds, i.e., $a_x\in [p,z]_{\gamma_2}$. 
But then actually $a_x\in [p,x]_{\gamma_2}$ by a similar argument as used above in the first case. 
Therefore $d(p,x)=d(p,a_x)+d(a_x,x)$, and so \eqref{eq:d(x,p)d(x,a_x)} implies that $d(p,a_x)=0$, 
a contradiction. 

Finally, if all the assumptions in the formulation of the proposition 
hold with the exception of the assumption that $d(x,A)$ is attained for all $x\in [p,z)$,  
but instead it is known that $(X,d)$ is proper and $A$ is closed, then a simple well known 
argument based on compactness (without referring to any geodesic assumption on $X$) 
shows that $d(x,A)$ is actually attained for all $x\in X$ and hence the conclusion in the formulation 
of the proposition holds. 
\end{proof}

\begin{proof}[{\bf Proof of Proposition \ref{prop:StrictInt} } ]
From the inclusion $\{x\in X: d(x,P)<d(x,A)\}\subseteq \textnormal{dom}(P,A)$ and the fact that $\textnormal{dom}(P,A)$ 
is closed (Lemma ~\ref{lem:dom}\eqref{item:closure}) it follows that the inclusion $\overline{\{x\in X: d(x,P)<d(x,A)\}}\subseteq \textnormal{dom}(P,A)$ 
holds. For the second inclusion, let $z\in \textnormal{dom}(P,A)$. If $d(z,P)<d(z,A)$, then obviously the point $z$ is in the set $\overline{\{x\in X: d(x,P)<d(x,A)\}}$. Now suppose that $d(z,P)=d(z,A)$.
By assumption there is $p\in P$ such that $d(z,P)=d(z,p)$, and since $p\in P$, it  follows that  $p\notin A$. This 
and the fact that $d(p,A)$ is attained imply that $d(p,A)>0$. Hence if $z=p$, then $z\in \overline{\{x\in X: d(x,P)<d(x,A)\}}$.  If $z\neq p$, then $[p,z)\neq \emptyset$, and Proposition  \ref{prop:StrictSegment} implies that every $x\in [p,z)$ (arbitrary close to $z$) satisfies the inequality $d(x,P)\leq d(x,p)<d(x,A)$. Thus again $z\in \overline{\{x\in X: d(x,P)<d(x,A)\}}$.
\end{proof}

\begin{proof}[{\bf Proof of Proposition \ref{prop:BoundaryInterior}}]
Let $f:X\to \mathbb{R}$ be defined by  $f(x)=d(x,P)-d(x,A)$ for all $x\in X$.  
Because of the continuity of $f$ and because $\textnormal{dom}(P,A)=\{x\in X: f(x)\leq 0\}$  
one obtains that  that $\textnormal{dom}(P,A)$ is closed and that all the points 
in $\{x\in X: d(x,P)<d(x,A)\}$ and $\{x\in X: d(x,A)<d(x,P)\}$ 
are interior points. It follows that $\partial(\textnormal{dom}(P,A))\subseteq \{x\in X: d(x,P)=d(x,A)\}$ 
without any assumption on the sites or on $X$. For the reverse inclusion, let 
 $z$ satisfy $d(z,P)=d(z,A)$. Then $z\in \textnormal{dom}(P,A)$. Let $a\in A$ satisfy $d(z,a)=d(z,A)$. 
 Then $a\neq z$ since otherwise the equality $d(z,P)=d(z,A)$ and the fact that $d(z,P)$ is 
 attained would imply that $z\in P\cap A$, a contradiction. Hence $[a,z)\neq\emptyset$. The 
 inclusion  $[a,z)\subseteq \{x\in X: d(x,A)<d(x,P)\}$ holds  
because of Proposition  ~\ref{prop:StrictSegment} (with $P$ instead of $A$ 
and $a$ instead of $p$) and it proves that arbitrary close to $z$ 
there are points outside $\textnormal{dom}(P,A)$. Thus 	\eqref{eq:boundary} holds. 
Finally, \eqref{eq:boundary}, the equality 
\begin{multline*}
\partial(\textnormal{dom}(P,A))\bigcup \textnormal{int}(\textnormal{dom}(P,A))=\textnormal{dom}(P,A)\\= \{x\in X: d(x,P)=d(x,A)\}\bigcup\{x\in X: d(x,P)<d(x,A)\},    
\end{multline*}
 and the fact that the terms in both unions are disjoint, all imply \eqref{eq:interior}. 
\end{proof}

The next two lemmas have a somewhat independent nature but they are needed for later use. 
They are  probably known, and their proofs are given for the sake of completeness. 
\begin{lem}\label{lem:DistIntersect}
Let $(X,d)$ be a proper metric space and let $(A_{\gamma})_{\gamma=1}^{\infty}$ be a decreasing sequence 
of closed subsets of $X$ such that $A:=\bigcap_{\gamma=1}^{\infty}A_{\gamma}$ is nonempty.  
Then $d(x,A)=\lim_{\gamma\to\infty}d(x,A_{\gamma})$  for each $x\in X$. 
\end{lem}
\begin{proof}
Let $x\in X$. The definition of properness implies that any bounded sequence in $X$ has a convergent 
subsequence.  Thus the distance between $x$ and any nonempty closed subset is attained. Therefore for each $\gamma$ there exists $x_{\gamma}\in A_{\gamma}$ such that
$d(x,A_{\gamma})=d(x,x_{\gamma})$. Since $(A_{\gamma})_{\gamma=1}^{\infty}$ is decreasing and $A_{\gamma}$ is closed
for each $\gamma$, it follows that any limit point of the sequence $(x_{\gamma})_{\gamma\in \mathbb{N}}$ is in $A_{\gamma}$ for each $\gamma$ and hence in  $A$. Since 
$(A_{\gamma})_{\gamma=1}^{\infty}$ is decreasing, the sequence $(d(x,x_{\gamma}))_{\gamma\in \mathbb{N}}$ is increasing and its limit is bounded by $d(x,A)$. Hence $(x_{\gamma})_{\gamma\in \mathbb{N}}$ is bounded and  there exists $x_{\infty}\in A$ and a subsequence such that $x_{\infty}=\lim_{\beta\to\infty}x_{\gamma_{\beta}}$. From the continuity of the distance function we deduce that $d(x,x_{\infty})=\lim_{\beta\to\infty}d(x,x_{\gamma_{\beta}})\leq d(x,A)$. But
$d(x,A)\leq d(x,x_{\infty})$ by definition. Therefore  $\lim_{\gamma\to\infty}d(x,A_{\gamma})=\lim_{\beta\to\infty}d(x,A_{\gamma_{\beta}})=d(x,x_{\infty})=d(x,A)$.
\end{proof}

\begin{lem}\label{lem:HausdorffUnionIntersection}
Let $(X,d)$ be a compact metric space. Let $(A_n)_{n=1}^{\infty}$ be a decreasing sequence  of nonempty closed sets of $X$, and let $(B_n)_{n=1}^{\infty}$ be an increasing  sequence of nonempty sets of $X$. Let $A=\bigcap_{n=1}^{\infty}A_n$ and $B=\overline{\bigcup_{n=1}^{\infty}B_n}$. Then
$A=\lim_{n\to \infty}A_n$ and $B=\lim_{n\to\infty}B_n$ with respect to the Hausdorff distance.
\end{lem}
\begin{proof}
Suppose by way of contradiction that it is not true that $\lim_{n\to\infty}D(A_n,A)=0$, where $D$ is the Hausdorff distance. Then there exists some $\epsilon>0$ and a subsequence $n_m\in \mathbb{N}$ such that $\epsilon\leq D(A_{n_m},A)$ for each $m\in\mathbb{N}$. Since $A\subseteq A_n$ for each $n$ and because $X$ is compact, the definition of the Hausdorff  distance implies that for each $m$ there exists $x_m\in A_{n_m}$ such that $\epsilon \leq d(x_m,A)$. Let $(x_{m_l})_{l\in\mathbb{N}}$  be a convergent subsequence of the sequence $(x_m)_{m=1}^{\infty}$, and let $x$ be its limit. Since the sequence $(A_n)_n$ is decreasing, it follows that $x_{m_{l'}}\in A_{n_{m_l}}$ whenever $l\leq l'$. Therefore $x\in A_{n_{m_l}}$ because $A_{n_{m_l}}$ is closed. Thus $x\in \bigcap_{l=1}^{\infty}A_{n_{m_l}}$.  But $\bigcap_{l=1}^{\infty}A_{n_{m_l}}=A$ because the intersection is decreasing. Hence  $0=d(x,A)=\lim_{l\to\infty}d(x_{m_l},A)\geq \epsilon$, a contradiction.

Now suppose by way of contradiction that it is not true that $\lim_{n\to\infty}D(B_n,B)=0$. Then there exists some $\epsilon>0$ and a subsequence $n_m\in \mathbb{N}$ such that $\epsilon\leq D(B_{n_m},B)$ for each $m\in\mathbb{N}$. Since $B_n\subseteq B$ for each $n$ and $X$ is compact, the definition of the Hausdorff  distance implies that for each $m$ there exists $x_m\in B$ such that $\epsilon \leq d(x_m,B_{n_m})$.  Let $(x_{m_l})_{l\in\mathbb{N}}$  be a convergent subsequence of the sequence $(x_m)_{m=1}^{\infty}$, and let $x$ be its limit. Then $d(x,x_{m_l})<\epsilon/4$ for all $l$ large enough. Since $B$ is closed it follows that $x\in B$, so by the definition of $B$ there exists $n_0\in \mathbb{N}$ such that $d(x,B_{n_0})<\epsilon/2$. But $(B_n)_{n=1}^{\infty}$ is increasing, so  $d(x,B_n)<\epsilon/2$ for all $n\in \mathbb{N}$ such that $n\geq n_0$, and, in particular,  $d(x,B_{n_{m_l}})<\epsilon/2$ for $l$ large enough. Thus
\begin{equation*}
\epsilon\leq d(x_{m_l},B_{n_{m_l}})\leq d(x_{m_l},x)+d(x,B_{n_{m_l}})<3\epsilon/4,
\end{equation*}
a contradiction.
\end{proof}

In the next claims several subtler properties of $\textnormal{dom}$ and $\textnormal{Dom}$ will be derived.  Lemma ~\ref{lem:AntiIntersec_dom} and
Proposition ~\ref{prop:mDomM} generalize \cite[Lemma ~3.1]{AMTn} and \cite[Lemma ~5.1]{AMTn} respectively, and are partly inspired by them. 
\begin{lem}\label{lem:Dom3}
Let $(X,d)$ be a metric space and let $P=(P_k)_{k\in K}$ be a tuple of nonempty subsets of $X$ such that for all $k\in K$,
\begin{equation*}
r_k:=\inf\{d(P_k,P_j):j\neq k\}>0.
\end{equation*}
For $Q\subseteq X$ nonempty and $r>0$, let $B(Q,r)=\{x\in X:d(x,Q)<r\}$. We denote by $(\textnormal{Dom} P)_k$ the $k$-th  component of $\textnormal{Dom} P$. Then:
\begin{enumerate}[(a)]
\item\label{item:B(P,r)Dom2} $B(P_k,r_k/4)\subseteq (\textnormal{Dom}^2 P)_k\subseteq (\textnormal{Dom}^{\gamma} P)_k$ for each integer $\gamma\geq 1$ and any $k\in K$.
\item \label{item:r} If $X$ is a geodesic metric space, then $d(x,B(Q,r))+r\leq d(x,Q)$ for each $Q\subseteq X$ nonempty, $r>0$ and $x\notin B(Q,r)$.
\item\label{item:Dom3Disjoint} If $d(x,B(P_k,r))<d(x,P_k)$ for all $k\in K$, $r>0$, and  $x\notin \overline{P_k}$, then
for each integer $\gamma\geq 2$ the components of $\textnormal{Dom}^{\gamma}P $ are disjoint. Moreover, if $X$ is a geodesic metric 
 space, then for all $j,k\in K$,$j\neq k$ and all $\gamma\geq 2$ we have $r_k/8+r_j/8\leq d((\textnormal{Dom}^{\gamma}P)_k,(\textnormal{Dom}^{\gamma} P)_j)$.  
\item\label{item:r_k} Suppose  that $X$ is a geodesic metric space. Given $j,k\in K$, $j\neq k$, and an integer $\gamma\geq 2$, if $x\in (\textnormal{Dom}^{\gamma}P)_j$, then $r_k/4\leq d(x,P_k)$.
\end{enumerate}
\end{lem}
\begin{proof}
\begin{enumerate}[(a)]
\item Let $k\in K$ and suppose that $x\in B(P_k,r_k/4)$, i.e., $d(x,P_k)< r_k/4$.
By definition we have $(\textnormal{Dom}^2 P)_k=\textnormal{dom}(P_k,\bigcup_{j\neq k}\textnormal{dom}(P_j,\bigcup_{i\neq j}P_i))$, so in order to
prove that $x\in (\textnormal{Dom}^2 P)_k$ it suffices to prove that $r_k/4\leq d(x,y)$ for all $y\in \bigcup_{j\neq k}\textnormal{dom}(P_j,\bigcup_{i\neq j}P_i)$. Given $y$ in this union, it is in $\textnormal{dom}(P_j,\bigcup_{i\neq j}P_i)$ for some $j\neq k$. 
Hence 
\begin{equation*}
d(y,P_j)\leq d(y,\bigcup_{i\neq j}P_i)\leq d(y,P_k)\leq d(y,x)+d(x,P_k).
\end{equation*}
Therefore, by the definition of $r_k$,
\begin{equation*}
r_k\leq d(P_j,P_k)\leq d(P_j,y)+d(y,x)+d(x,P_k)\leq 2d(x,y)+2r_k/4,
\end{equation*}
i.e., $d(x,P_k)<r_k/4\leq d(x,y)$ and the assertion follows. Finally, Lemma \ref{lem:DomPk}  implies that $\textnormal{Dom}^2 P\subseteq \textnormal{Dom}^{\gamma} P$ for any integer $\gamma\geq 1$, and hence $(\textnormal{Dom}^2 P)_k\subseteq (\textnormal{Dom}^{\gamma} P)_k$ for each $k\in K$.
\item
Let $p\in Q$. Since $x\notin B(Q,r)$, 
a simple argument  shows that the intersection of the compact segment $[x,p]$  with the closed set $\partial B(Q,r)$  is attained at some point $y\in [x,p]$, e.g., at $y=\gamma(t)$ for $t=\inf\{u\in [0,d(x,p)]: \gamma(u)\in B(Q,r)\}$ 
(where $\gamma:[0,d(x,p)]\to X$, is the isometric function which maps $[0,d(x,p)]$ onto $[x,p]$). It must be that $s:=r-d(y,p)\leq 0$. If this is not true, then for any $z$ in the open ball $B(y,s)$ we have $d(z,Q)\leq d(z,p)\leq d(z,y)+d(y,p)<s+d(y,p)=r$,  and hence $y$ is in the interior of $B(Q,r)$, a contradiction.  
Thus, since $y\in [x,p]$, we obtain $d(x,p)=d(x,y)+d(y,p)\geq d(x,y)+r$, and so
\begin{equation*}
d(x,B(Q,r))=d(x,\overline{B(Q,r)})\leq d(x,y)\leq d(x,p)-r,
\end{equation*}
for any $p\in Q$. As a result, $d(x,B(Q,r))+r\leq d(x,Q)$.
\item Suppose that we know that the components of $\textnormal{Dom}^{3}P $ are disjoint.
Then for each integer $\gamma\geq 2$ the components of $\textnormal{Dom}^{\gamma}P $ are disjoint,  because $\textnormal{Dom}^{\gamma}P\subseteq \textnormal{Dom}^{3}P$ by Lemma  \ref{lem:DomPk}.
We now give a proof for the case $\gamma=3$. 
Let $k_1\neq k_2$ be two indices in $K$, and assume to the contrary that $x\in (\textnormal{Dom}^3 P)_{k_1}\bigcap (\textnormal{Dom}^3 P)_{k_2}$.  
By definition and by part \eqref{item:B(P,r)Dom2} we have
\begin{multline*}
\quad\quad d(x,P_{k_1})\leq d(x,\bigcup_{j\neq k_1}(\textnormal{Dom}^2 P)_j) \leq d(x,(\textnormal{Dom}^2 P)_{k_2})\\
\leq d(x,B(P_{k_2},r_{k_2}/4))\leq d(x,P_{k_2}).
\end{multline*}
Hence $x\notin \overline{P_{k_2}}$, because otherwise also $x\in \overline{P_{k_1}}$, a contradiction with the assumption that $0<r_{k_1}\leq d(P_{k_1},P_{k_2})$. Thus $d(x,B(P_{k_2},r_{k_2}/4))< d(x,P_{k_2})$ by assumption, so actually $d(x,P_{k_1})<d(x,P_{k_2})$. In the same way $d(x,P_{k_2})< d(x,P_{k_1})$, a contradiction.

Finally, suppose that $X$ is a geodesic metric space. 
Let $j,k\in K$ be different. Since $(\textnormal{Dom}^{\gamma}P)_i\subseteq (\textnormal{Dom}^3 P)_i$ for each integer $\gamma\geq 2$ and each $i\in K$, it suffices to show that $r_k/8+r_j/8\leq d(x,y)$ for any $x\in (\textnormal{Dom}^3 P)_{k}$ and $y\in(\textnormal{Dom}^3 P)_j$.
 By definition, the triangle inequality and parts \eqref{item:B(P,r)Dom2},\eqref{item:r},
\begin{multline*}
d(x,P_k)\leq d(x,\bigcup_{i\neq k}(\textnormal{Dom}^2 P)_i) \leq d(x,(\textnormal{Dom}^2 P)_j)
\leq d(x,B(P_j,r_{j}/4))\\
\leq d(x,P_j)-r_j/4\leq d(x,y)+d(y,P_j)-r_j/4.
\end{multline*}
In the same way $d(y,P_j)\leq d(y,x)+d(x,P_k)-r_k/4$. By adding these two inequalities we obtain the desired conclusion.
\item If $x\in (\textnormal{Dom}^{\gamma}(P))_j$, then $x\notin B(P_k,r_k/4)\subseteq (\textnormal{Dom}^{\gamma}(P))_k$ by parts  \eqref{item:B(P,r)Dom2} and \eqref{item:Dom3Disjoint}. Hence  $r_k/4\leq d(x,P_k)$ by part \eqref{item:r}.
\end{enumerate}
\end{proof}


\begin{lem}\label{lem:AntiIntersec_dom}
Let $(X,d)$ be a metric space, let $P\subseteq X$ be nonempty and suppose that $\{C_{\gamma}\}_{\gamma=1}^{\infty}$ 
is a family of subsets of $X$ such that $\bigcap_{\gamma=1}^{\infty} C_{\gamma}\neq
\emptyset$. If
\begin{equation}\label{eq:domInLess}
\textnormal{dom}(P,\bigcap_{\gamma=1}^{\infty} C_{\gamma})= \overline{\{x\in X: 
d(x,P)<d(x,\bigcap_{\gamma=1}^{\infty} C_{\gamma})\}},
\end{equation}
and
\begin{equation}\label{eq:ALimAk}
d(y,\bigcap_{\gamma=1}^{\infty} C_{\gamma})=\limsup_{\gamma\to \infty}d(y,C_{\gamma}),\quad \forall y\in X,
\end{equation} then
\begin{equation}\label{eq:AntiCont}
\textnormal{dom}(P,\bigcap_{\gamma=1}^{\infty} C_{\gamma})=\overline{\bigcup_{\gamma=1}^{\infty}\textnormal{dom}(P,C_{\gamma})}.
\end{equation}
\end{lem}
\begin{proof}
By antimonotonicity, $\textnormal{dom}(P,\bigcap_{j=1}^{\infty}
C_j)\supseteq\textnormal{dom}(P,C_{\gamma})$ for all $\gamma$, and since the left hand side 
is closed, the inclusion $\textnormal{dom}(P,\bigcap_{\gamma=1}^{\infty}
C_{\gamma})\supseteq\overline{\bigcup_{\gamma=1}^{\infty}\textnormal{dom}(P,C_{\gamma})}$ follows. For 
the reverse inclusion, let $\epsilon>0$ be given and suppose that
$x\in \textnormal{dom}(P,\bigcap_{\gamma=1}^{\infty} C_{\gamma})$. We should prove that there 
are $\gamma\in \mathbb{N}$ and $y\in \textnormal{dom}(P,C_{\gamma})$ such that $d(y,x)<\epsilon$.

By \eqref{eq:domInLess} there is $y\in X$ such that
$d(x,y)<\epsilon$ and $r:=d(y,\bigcap_{\gamma=1}^{\infty} C_{\gamma})-d(y,P)>0$,
and by \eqref{eq:ALimAk} there is $\gamma$ large enough such that
$|d(y,\bigcap_{j=1}^{\infty} C_j)-d(y,C_{\gamma})|<r/2$. Hence
\begin{equation*}
d(y,P)+r/2<d(y,\bigcap_{\gamma=1}^{\infty} C_{\gamma})<d(y,C_{\gamma})+r/2.
\end{equation*}
Thus $d(y,P)<d(y,C_{\gamma})$, and so $y\in \textnormal{dom}(P,C_{\gamma})$. 
\end{proof}


\begin{prop}\label{prop:mDomM} 
Let $(X,d)$ be a geodesic metric space and let $(P_k)_{k\in K}$ be a tuple of nonempty closed  sets in $X$ satisfying  $r:=\inf\{d(P_k,P_j): j,k\in K,\, j\neq k\}>0$. For each nonnegative integer $\gamma$ 
and each $k\in K$ let 
\begin{equation}\label{eq:R_gammaA_k}
R^{\gamma}=\textnormal{Dom}^{2\gamma+1}(P_k)_{k\in K}, \quad 
A_{\gamma,k}=\bigcup_{j\neq k} R_j^{\gamma},\quad A_k=\bigcap_{\gamma=1}^{\infty}A_{\gamma,k}.
\end{equation}
\begin{enumerate}[(a)]
\item\label{item:UnionIntersect} The following equality holds:
\begin{equation}\label{eq:UnionIntersect}
\bigcup_{j\neq  k}\bigcap_{\gamma=1}^{\infty}R_j^{\gamma}=\bigcap_{\gamma=1}^{\infty}\bigcup_{j\neq k} R_j^{\gamma}
\end{equation}
\item\label{item:AkNonEmpty} $A_{\gamma,k}$ is a closed set  for each $k\in K$ and each $\gamma\geq 1$. 
$A_k$ is nonempty, closed, and satisfies $d(P_k,A_k)\geq r/4$  for each $k\in K$. 
\item\label{item:DomCapCup} If for each $k\in K$ both \eqref{eq:domInLess} and \eqref{eq:ALimAk}  hold with $C_{\gamma}=A_{\gamma,k}$, then 
\begin{equation}\label{eq:DomCapCup}
\textnormal{Dom}(\bigcap_{\gamma=0}^{\infty} \textnormal{Dom}^{2\gamma+1}(P_k)_{k\in K})=
\overline{\bigcup_{\gamma=0}^{\infty}\textnormal{Dom}^{2\gamma}(P_k)_{k\in K}}.
\end{equation}
\item\label{item:StrictlyConvexProper} If  $(X,d)$ is  proper and has the geodesic 
inclusion property, then \eqref{eq:domInLess}-\eqref{eq:ALimAk} (with $C_{\gamma}=A_{\gamma,k}$, $k\in K$ arbitrary) 
and hence \eqref{eq:DomCapCup} hold. 
\end{enumerate}
\end{prop}


\begin{proof}
\begin{enumerate}[(a)]
\item 
If $x\in \bigcup_{j\neq  k}\bigcap_{\gamma=1}^{\infty}R_j^{\gamma}$, then
$x\in \bigcap_{\gamma=1}^{\infty}R_j^{\gamma}$ for some $j\neq k$. Since $R_j^{\gamma}\subseteq \bigcup_{i\neq k} R_i^{\gamma}$, it follows that
$\bigcap_{\gamma=1}^{\infty}R_j^{\gamma}\subseteq \bigcap_{\gamma=1}^{\infty}\bigcup_{i\neq k} R_i^{\gamma}$, so $x\in \bigcap_{\gamma=1}^{\infty}\bigcup_{i\neq k} R_i^{\gamma}$. On the other hand, let $x\in \bigcap_{\gamma=1}^{\infty}\bigcup_{j\neq k} R_j^{\gamma}$. Given  $\gamma\in \mathbb{N}$, Lemma \ref{lem:Dom3}\eqref{item:Dom3Disjoint} (note that $2\gamma+1\geq 3$ and 
$R^{\gamma}=\textnormal{Dom}^{2\gamma+1}(P_k)_{k\in K}$) implies that  $R_i^{\gamma}\bigcap R_j^{\gamma}=\emptyset$ for $j\neq i$, so there exists exactly one  index $j\neq k$ such that $x\in R_j^{\gamma}$. It must be that all these indices coincide.    If this is not true, then $x\in R_j^{\gamma}\bigcap R_{j'}^{\gamma'}$ for some $\gamma'>\gamma$ and corresponding indices $j'\neq j$. But  $R_{j'}^{\gamma'}\subseteq R_{j'}^{\gamma}$ by Lemma \ref{lem:DomPk}, so $x\in R_j^{\gamma}\bigcap R_{j'}^{\gamma}$, a contradiction to Lemma \ref{lem:Dom3}\eqref{item:Dom3Disjoint}. As a result
$x\in \bigcap_{\gamma=1}^{\infty}R_j^{\gamma}$ for some $j\neq k$, and this establishes  \eqref{eq:UnionIntersect}.

\item By Lemma \ref{lem:DomPk} we have  $P_j\subseteq \bigcap_{\gamma=1}^{\infty}R_j^{\gamma}$, and so  \eqref{eq:UnionIntersect} implies that $\emptyset\neq\bigcup_{j\neq k}P_j\subseteq A_k$. 
By Lemma \ref{lem:Dom3}\eqref{item:r_k} it follows 
that $0<r/4\leq d(P_k,A_{\gamma,k})\leq d(P_k,A_k)$ for each $\gamma\in \mathbb{N}$ and $k\in K$. By Lemma \ref{lem:DomPk} 
the intersection which defines $A_k$ is decreasing. By Lemma \ref{lem:dom}\eqref{item:closure} and Lemma \ref{lem:Dom3}\eqref{item:Dom3Disjoint} each $A_{\gamma,k}$ is a closed set, because it is a union of closed and disjoint sets with a positive distance (at least $r/4$)  between any two  different members in the union. 
 Therefore $A_k$ is the intersection of closed sets and thus closed.

\item By part \eqref{item:AkNonEmpty} we know that $A_k\neq \emptyset$ for any $k\in K$. Given $k\in K$, our assumption 
 implies that \eqref{eq:domInLess}-\eqref{eq:ALimAk} hold with $C_{\gamma}=A_{\gamma,k}$ for 
each $\gamma\in \mathbb{N}$. By part \eqref{item:UnionIntersect} we know that \eqref{eq:UnionIntersect} 
holds. This, together with  Lemma  \ref{lem:DomPk}, Lemma \ref{lem:AntiIntersec_dom},  the definition of $\textnormal{Dom}$, 
and the fact that intersection,union, and closure on tuples are taken component-wise, imply   
\begin{multline*}
\textnormal{Dom}\left(\bigcap_{\gamma=0}^{\infty}\textnormal{Dom}^{2\gamma+1}(P_k)_{k\in K}\right)=\textnormal{Dom}\left(\bigcap_{\gamma=0}^{\infty}R^{\gamma}\right)
=\textnormal{Dom}\left(\bigcap_{\gamma=1}^{\infty}R^{\gamma}\right)\\
=\left(\textnormal{dom}\left(P_k,\bigcup_{j\neq k}\bigcap_{\gamma=1}^{\infty}R_j^{\gamma}\right)\right)_{k\in K}
=\left(\textnormal{dom}\left(P_k,\bigcap_{\gamma=1}^{\infty}\bigcup_{j\neq k}R_j^{\gamma}\right)\right)_{k\in K}\\
=\left(\overline{\bigcup_{\gamma=1}^{\infty}\textnormal{dom}(P_k,\bigcup_{j\neq k}R_j^{\gamma})}\right)_{k\in K}
=\overline{\left(\bigcup_{\gamma=1}^{\infty}\left(\textnormal{Dom}(R^{\gamma})\right)_k\right)_{k\in K}}\\
=\overline{\bigcup_{\gamma=1}^{\infty}\textnormal{Dom}(R^{\gamma})}
=\overline{\bigcup_{\gamma=0}^{\infty}\textnormal{Dom}^{2\gamma}(P_k)_{k\in K}}.
\end{multline*}

\item Since $A_k$ is a decreasing intersection (by Lemma \ref{lem:DomPk}) of closed subsets (part \eqref{item:AkNonEmpty}), 
Lemma  \ref{lem:DistIntersect} implies that $d(x,A_k)=\lim_{\gamma\to\infty}d(x,A_{\gamma,k})$ for each $x\in X$ 
and hence \eqref{eq:ALimAk} holds. 
Because $d(P_k,A_k)\geq r/4>0$ by part \eqref{item:AkNonEmpty}, we conclude from Proposition ~\ref{prop:StrictInt} that \eqref{eq:domInLess} holds for $P=P_k$ and $A=A_k$. 
From part \eqref{item:DomCapCup} we conclude that \eqref{eq:DomCapCup} holds. 
\end{enumerate}
\end{proof}



After the above long preparation it is possible to prove the assertions formulated in 
Section ~\ref{sec:ConvergenceResult}. 

\begin{proof}[{\bf Proof of Theorem \ref{thm:Metric}}]
By Lemma~\ref{lem:GeneralDistance}\eqref{item:DomUnionIntersect} 
and Lemma  \ref{lem:dom}\eqref{item:DomClosure} it follows that 
\begin{equation}\label{eq:Dom(m)M}
\textnormal{Dom}(m)=\textnormal{Dom}\left(\overline{\bigcup_{\gamma=0}^{\infty}\textnormal{Dom}^{2\gamma}(P)}\right) 
=\bigcap_{\gamma=0}^{\infty}\textnormal{Dom}^{2\gamma+1}(P)=M.
\end{equation}
\end{proof}
\begin{proof}[{\bf Proof of Theorem \ref{thm:AlgUniConvex}}]
By assumption $(X,d)$ is a proper geodesic metric space which has 
the geodesic inclusion property and  the sites satisfy \eqref{eq:d(P_k,P_j)>0}. 
Proposition  \ref{prop:mDomM}\eqref{item:StrictlyConvexProper} implies that 
\begin{equation}\label{eq:Dom(M)m}
\textnormal{Dom}(M)=\textnormal{Dom}\left(\bigcap_{\gamma=0}^{\infty}\textnormal{Dom}^{2\gamma+1}(P)\right)
=\overline{\bigcup_{\gamma=0}^{\infty}\textnormal{Dom}^{2\gamma}(P)}=m.
\end{equation}
Thus $m=\textnormal{Dom}(M)$. This and Theorem ~\ref{thm:Metric} imply 
that $M=\textnormal{Dom}(m)=\textnormal{Dom}^2(M)$ and $m=\textnormal{Dom}(M)=\textnormal{Dom}^2(m)$. Hence both $m$ and $M$ are double zone diagrams. From Lemma ~\ref{lem:DomPk} 
and Lemma ~\ref{lem:dom}\eqref{item:DomClosure} it can be seen that $m$ and $M$ are 
respectively the least and greatest double zone diagrams. Finally, if $|K|=2$, then 
by \eqref{eq:Dom(m)M}, \eqref{eq:Dom(M)m}, and the definition of $\textnormal{Dom}$ it follows that  $(M_1,M_2)=\textnormal{Dom}(m)=((\textnormal{dom}(P_1,m_2),\textnormal{dom}(P_2,m_1))$ and $(m_1,m_2)=(\textnormal{dom}(P_1,M_2),\textnormal{dom}(P_2,M_1))$. 
Hence $m_1=\textnormal{dom}(P_1,M_2)$ and $M_2=\textnormal{dom}(P_2,m_1)$. But $\textnormal{Dom}(m_1,M_2)=(\textnormal{dom}(P_1,M_2),\textnormal{dom}(P_2,m_1))$ 
by the definition of $\textnormal{Dom}$. Therefore $\textnormal{Dom}(m_1,M_2)=(m_1,M_2)$, i.e., $(m_1,M_2)$ is a zone 
 diagram.  In the same way $(M_1,m_2)$ is a zone diagram. 
\end{proof}


\begin{proof}[\bf Proof of Corollary \ref{cor:CompactLimit}]
By Lemma \ref{lem:DomPk} the sequence  $A_{\gamma}=(\textnormal{Dom}^{2\gamma+1}(P))_k$ is decreasing and the sequence  $B_{\gamma}=(\textnormal{Dom}^{2\gamma}(P))_k$ is increasing for each $k\in K$. As a result, when $X$ is compact \eqref{eq:Mm_k} follows from Lemma  \ref{lem:HausdorffUnionIntersection}. 
\end{proof}


\begin{proof}[{\bf Proof of Corollary \ref{cor:ZoneDoubleZone}}]
If the least and the greatest double zone diagram coincide,  
then, without any restriction on the sites or the space, there exists a unique zone diagram  which 
is equal to both of them \cite[Corollary 6.2]{ReemReichZone}. 
 Since in the proof of \cite[Theorem 1.1]{KMT2012} 
it was proved that the  least double zone diagram coincides with the greatest one 
when $(X,d)$ is a finite dimensional Euclidean space, the assertion follows from Theorem \ref{thm:AlgUniConvex}. 
\end{proof}



\section*{Acknowledgments}
{\bf To be added.}



\bibliographystyle{amsplain}
