\input CS445Hdr
\def\file{CS445N8B.tex of November 7, 1982 7:54 PM}
\def\header{8. Closest-point problems}
\def\title{Lectures 17--18}

\sect{8.4. The Delaunay diagram.}

We can represent the adjacency structure of $\Vscr(P)$
by drawing a line segment $pq$ between the points $p$ and
$p$ whenever the corresponding regions $Vp̌$ and $Vq̌$
have an edge in common.
The set of these segments ({\it dual edges})
is called the {\it dual Voronoi diagram}, or {\it Delaunay
diagram}, of the set $P$, and will be denoted by
$\Dscr(P)$.\fig5
Note that there is a one-to-one correspondence between the
regions of $\Vscr(P)$ and the vertices of $\Dscr(P)$, and also
between the edges of the two.\par

The Delaunay diagram for $P$ seems to contain
less information than $\Vscr(P)$; for example,
the Voronoi vertices are not explicitly represented in the
former. For computational puproses, however, the two are
pretty much equivalent, since the missing information
can generally be computed from $\Dscr(P)$, without
increasing the running time or the space requirements by
more than a constant factor. Actually, algorithms
become simpler and more efficient when they
are rewritten so as to use the Delaunay diagram in place
of the full Voronoi, since the data structure overhead is
smaller and some information doesn’t have to be computed
at all. One major drawback of doing so is that
$\Dscr(P)$ is harder to visualize than $\Vscr(P)$,
and therefore the algorithms become less obvious.\par

The following lemma gives an important property of the
Delaunay diagram:

\lemma{8.5}{Let $p,q \in P$ be connected in $\Dscr(P)$, and
let $r,s\in P$ be such that $Q=\langle p,r,q,s\rangle$ is a
convex quadrilateral. Then the internal angles of $Q$
satisfy the relation $\hat p+\hat q > \pi > \hat r+\hat s$.}
\proof{If $p$ and $q$ are connected in $\Dscr(P)$, $Vp̌$
and $Vq̌$ share an edge $e{̌pq}$ in $\Vscr(P)$. Let $x$
be a point of $e{̌pq}$ that is not a Voronoi vertex,
and let $C$ be the circle with center $x$ and radius
$\bar{xp}$.\fig4
Since $x$ is on $B(p,q)$ but not a vertex, the circle must pass
through $q$, and all other data points, including $r$ and $s$,
must lie strictly outside $C$.\par

\hp Since $Q$ is convex, the points $r$ and $s$ must lie on
opposite sides of the line $pq$. Let $\alpha$ be the arc
of $C$ with endpoints $pq$ that is on the same side as $s$;
since $r$ is outside $C$, the internal angle $\hat r$ of $Q$ at
$r$ is strictly less than $\alpha\over 2$; similarly, the angle
$\hat s$ is stricly less than $\pi-{\alpha\over2}$, and
we have $\hat r + \hat s < \pi$. Since the internal angles
of $Q$ add up to $2\pi$, the angles at $p$ and $q$ must add
up to more than $\pi$.}

This lemma has a surprising consequence:

\theorem{8.6}{Two distinct edges $pq$ and $rs$ of
$\Dscr(P)$ are either disjoint or have a single endpoint in
common.}
\proof{Suppose $p$ lies on $rs$ but is distinct from $r$ and
$s$. Then $B(r,p)$ and $B(s,p)$ are parallel, and $C(r,p)\inte
C(s,p)=\ety$. \fig4
Since $Vř\subset C(r,p)$ and $Vš\subset C(s,p)$, this
contradicts the hypothesis that $r$ and $s$ are connected.\par

\hp Now suppose $pq$ and $rs$ meet at some point $o$
that is none of the four endpoints.\fig4
This implies the quadrilateral $Q=\langle p,r,q,s\rangle$ is
convex; by lemma 8.5, we should have both
$\hat p+\hat q>\hat r + \hat s$,
and vice-versa --- again, a contradiction. Therefore,
if $pq$ and $rs$ have a point in common, it must
be an endpoint of both segments.}

Therefore, the edges of the Delaunay diagram do not cross!
This fact is all the more surprising if we note
that $pq$ may not cross $e{̌pq}$, but may cross
several other edges $e{̌p\pr q\pr}$ ``unrelated’’ to
it.\fig4\par

\digress {By considering only the incidence relationships
between vertices and edges, and ignoring their coordinates,
the Voronoi diagram $\Vscr(P)$ can be viewed as an
undirected planar graph; indeed, if we take into account
also the incidence between faces, vertices and edges,
$\Vscr(P)$ is a specific embedding of that graph in the plane,
i.e. a {\it planar map}.\par

\dp Given an embedding $M$ of a planar graph $G$,
we define a {\it dual graph} $G*$ of $G$ by
taking the faces of $M$ as vertices of $G*$, and drawing an
edge between $f$ and $g$ in $G*$ whenever $f$ and $g$
share an edge in $M$. (The dual of a graph is not uniquely
defined, for different embeddings
may give rise to non-isomorphic duals.)
We know from Graph Theory that
$G*$ is also a planar graph, and therefore
can be drawn in the plane without crossing edges.\par

\dp In view of this, theorem 8.6 may not seem so exciting,
since $\Dscr(P)$, viewed as a graph,
is clearly the dual of $\Vscr(P)$. However, the planarity
theorems from Graph Theory assume that
the edges of the dual can be drawn with curved lines,
and/or that the vertices of the dual can be positioned
arbitrarily in the plane. If we require that the edges be
straight lines, and that the dual vertices are inside the
corresponding faces of the original map, it may be impossible
to avoid crossing edges, as is the case for the map
below.\fig4
If nothing else, this shows that not every convex subdivision
of the plane is the Voronoi diagram of some set $P$.\par

\dp Let $M$ be an embedding of a planar graph $G$,
let $G*$ be the dual graph constructed from $M$, and
let $M*$ be an embedding of $G*$. We say that
$M*$ is the {\it dual map} of $M$ if and only if
the neighbors of a given vertex $f$ of $M*$,
in counterclockwise order, are the faces
of $M$ adjacent to the face $f$, in the same order.
In general, not all embeddings of $G*$ have this property,
and those which do are topologically equivalent.
The next few theorems show that $\Dscr(P)$ is
the dual of $\Vscr(P)$ as a map, not just as a graph.}

The structure of the Delaynay diagram is further
clarified by the following result:

\theorem{8.7}{The neighbors of $p$ in $\Dscr(P)$ are
$q1̌, q2̌, \ldotss, qm̌$, in counterclockwise
order, if and only if $V{̌q1̌}, V{̌q2̌}, \ldotss,
V{̌qm̌}$ are the regions adjacent to $Vp̌$, in
the same order.\fig4}
\proof{It follows immediately from the
definition of $\Dscr(P)$ that the neighbors of $p$ in
$\Dscr(P)$ are precisely those $q$ for which $Vp̌$ and
$Vq̌$ are adjacent, so we have only to show that the
clockwise ordering of both sets is the same.\par

\hp Orient each Voronoi edge $e{̌pqǐ}$
of $Vp̌$ counterclockwise around $p$, and orient the
corresponding dual edge $pqǐ$ away from $p$.
Then the direction of $e{̌pqǐ}$ is that of $pqǐ$,
rotated by $\pi\over2$ counterclockwise; therefore,
the angle between $e{̌pq1̌}$ and $e{̌pqǐ}$
is exactly the same as the angle $q1̌pqǐ$,
for all $i$.\fig3
Since $Vp̌$ is a convex polygon, the angles between
$e{̌pq1̌}$ and $e{̌pqǐ}$ (measured in the interval
$\left[0,2\pi\right)$ are strictly increasing if and only if the
edges of $Vp̌$ are considered in clockwise order.
Similarly, the angles $q1̌pqǐ$ are strictly increasing
if and only if the points $qǐ$ appear in counterclockwise
order around $p$. The theorem follows immediately.}

The result below may be seen as the ``dual’’
of theorem 8.7:

\theorem{8.8}{The regions incident to a finite vertex $v$
of $\Vscr(P)$ are $V{̌p1̌}, V{̌p2̌}, \ldotss,
V{̌pm̌}$, in counterclockwise order, if and only if
$p1̌, p2̌, \ldotss, pm̌$ are the vertices of an internal
face of $\Dscr(P)$, in the same order.\fig4}
\proof{Let’s first assume that the regions
around $v$ are $\Vscr(P)$ are $V{̌p1̌}, V{̌p2̌}, \ldotss,
V{̌pm̌}$, in that order. From lemma 8.3, we see that
$p1̌, p2̌, \ldotss, pm̌$ lie on a $P$-free
circle with center $v$. We see also that $pǐ$ and
$pǰ$ will be adjacent in $\Dscr(P)$ if and only if
$\leftv i-j\rightv=1$\foot1{
To simplify the indexing of the vertices, we assume
$p0̌\eqv pň$, $p1̌\eqv p{̌n+1}$, etc.}. It follows
$F=\langle p1̌, p2̌, \ldotss, pm̌\rangle$ is a cycle
of $\Dscr(P)$ containing no other data point and no
diagonals, i.e. $F$ is an internal face of $\Dscr(P)$.\par

\hp Conversely, let’s assume $F=\langle p1̌, p2̌,
\ldotss, pm̌\rangle$ is an internal face of $\Dscr(P)$.
For any $i$, the points $p{̌i+1}$ and $p{̌i-1}$ are
(counterclockwise)
consecutive neighbors of $pǐ$ in $\Dscr(P)$; therefore,
by theorem 8.7, $e{̌pǐ p{̌i+1}}$ and
$e{̌pǐ p{̌i-1}}$ are consecutive edges of $V{̌pǐ}$.
It follows that $V{̌p{̌i-1}}$, $V{̌pǐ}$ and $V{̌p{̌i+1}}$
are consecutive regions incident, in that order,
to some Voronoi vertex $vǐ$, and the same holds for all $i$.
\fig4
The regions $V{̌p{̌i}}$ and $V{̌p{̌i+1}}$ are
both incident to $vǐ$ and $v{̌i+1}$, and in
both vertices they appear in
the same counterclockwise order. This is possible only
if $vǐ=v{̌i+1}$. It follows that all regions are incident to
the same Voronoi vertex, in the given order.}

Combining this result with lemma 8.3, we see that

\theorem{8.9}{An internal face of $\Dscr(P)$ with $m$ sides
is a convex polygon, whose points lie on a circle centered at
a Voronoi vertex of degree $m$.}

A result similar to theorem 8.8
(and whose proof we leave to the reader) relates
the external face of $\Dscr(P)$ to the Voronoi vertex at
infinity $\omega$:

\theorem{8.10}{The infinite regions of $\Vscr(P)$
are $V{̌p1̌}, V{̌p2̌}, \ldotss, V{̌pm̌}$, in
{\it clockwise} order, if and only if $p1̌, p2̌,
\ldotss, pm̌$ are the vertices of the external face of
$\Dscr(P)$, in {it counterclockwis} order.}

\sect {8.6. Using the Voronoi diagram.}

Let’s postpone for the moment the question of
how to construct the Voronoi diagram, and
let’s discuss how it can be used to solve the closest-point
problem and its relatives. To find the data point $p$
that is closest to a given $x\in \Sscr$, all
we have to do is to locate the region $Vp̌$ that
contains $x$.

The Lee-Preparata algorithm (??? 5.???) can be used
for this purpose, since $\Vscr(P)$ is a convex
(and therefore monotone) subdivision of the plane.
Its query time is not optimal ($O(\log n)↑2$), but is
significantly beter than that of exhaustive search.\par

Kirkpatrick’s algorithm (5.???) has a better asymptotic
performance, just $O(\log n)$ time per query. As it is
formulated, Kirkpatrick’s algorithm requires the Voronoi
regions to be subdivided even further, so as to
obtain a triangulation. Two possible ways to
implement this subdivision are shown below\foot1{Note that
these subdivisions are {\it not} related to the Delaunay
triangulation $\Dscr*(P)$ in any simple way: the edges here
connect Voronoi vertices, not data points.}.\fig4
Given $\Vscr(P)$, we can construct either
of these triangulations (and the data structure required
by algorithm 5.???) in linear time, so the
preprocessing cost is still dominated by the computation
of the Voronoi diagram.\par

{\bf ??? Is there a location algorithm that does not
require the full Voronoi, only the Delaunay triang.?}

A common variant of the closest-point problem is the
{\it nearest neighbor problem}: for each point $p\in P$,
find the point $q\in P$ that is closest to it. The
straightforward solution, computing all $n-1$ distances from
each $p$ and selecting the minimum, gives a $\Theta(n↑2)$ algorithm. If we were to construct the Voronoi diagram for
$P-\set{p}$ and locate $p$ in it, the total work would be
even worse, $\Theta(n↑2\log n)$ in the worst case.\par

However, it is the case that $\Vscr(P)$ already contains
enough information to allow us to solve the nearest-neighbor
problem in linear time, as shown by the following result:

\lemma{8.11}{The point $q\in P$ that is closest to a given
point $p\in P$ is adjacent to $p$ in the Delaunay diagram for
$P$, and the midpoint of $pq$ is on the Voronoi edge
$e{̌pq}$.}
\proof{Let $C$ be the circle with diameter $pq$, and $C\pr$
be the one with center $p$ and radius $pq$.\fig3
Since $q$ is the closest neighbor of $p$, there are no
points of $P$ interior to $C\pr$; therefore, $p$ and $q$ are
the only points of $P$ on $C$. It follows that the center
$o$ of $C$, which is the midpoint of $pq$, is a point of
the edge $e{̌pq}$ and is not a Voronoi vertex,
so $p$ and $q$ are adjacent in $\Dscr(P)$.}

Therefore, to find the closest neighbor of $p$ we have to
look only at the vertices $q$ adjacent to it
in $\Dscr(P)$, and select the closest one.
This may take $\Theta(n)$ time for a single vertex $p$,
since its degree can be as high as $n-1$. The total number
of edges, however, is at most $3n-6$, so the computation of
all nearest neighbors can be done in only $O(n)$ time.
Another way to look at this result is to observe that,
although the {\it maximum} degree of a point in $\Dscr(P)$
is $n-1$, the {\it average} degree is $O(1)$ (more
precisely, $6\left(1-{2\over n}\right)$).
The computation of the Voronoi, requiring
$O(n\log n)$ time, will be the limiting step in this algorithm.\par

A slight variation of this algorithm also solves the {\it
closest pair problem}, namely finding the pair $p,q\in P$
whose distance is minimal.\par

Another application of the Voronoi diagram is in finding the
largest circle $C$ devoid of $P$.
If no constraints are imposed on $C$,
the problem has no solution: we can always place a circle
of arbitrary radius far enough from $P$ so that
it includes no points of $P$. If the center of $C$ is required
to lie in the covex hull of $P$, the problem becomes
meaningful; using theorem 8.???, we can easily prove
that the center of $C$ is either a Voronoi vertex interior
to the hull of $P$, or the intersection
of a Voronoi edge with an edge of the hull. It turns out that
all those points, and therefore the largest $C$, can be found
in $O(n)$ time {\bf I guess...}.

{\bf ??? Closest $R$-$S$ pair}

\sect{8.6 Data structures for the Voronoi diagram.}

To represent Voronoi diagrams, we can use any of the data structures for planar subdivisions
discussed in section 5.???D. As we observed
there, which data structure
is the best depends obviously on the operations we want to
perform on it; in particular on the the algorithm we are
going to use for point location.\par

Query processing is not the only operation that has to be
considered when selecting the data structure. Some applications may require the occasional insertion, deletion
or displacement of data points, or the merge of
two separate Voronoi diagrams. Even if such capabilities
are not explicitly called for by the application, they may be
required by the algorithm that computes the Voronoi in the
first place. The best data structures for query processing
(e.g, Kirpatrick’s hierarchy) may be quite inadequate for
such operations; so, unless the Voronoi-building algorithm
itself performs point location, it may be convenient
to use a simpler data structure to represent the diagram while
it is being constructed.\par

We can use the basic data structure described in section
5.???D. However, the special nature of the Voronoi diagram
allows us to introduce some simplifying modifications. First,
notice that
all internal angles are less than
$\pi$, so (as discussed in section 5.???D) we can compute
the vertices by the intersection of two adjacent edges,
and don’t have to store their coordinates in the structure. Furthermore,
the supporting line of each edge is the bissector of the data
points corresponding to the adjacent faces; if we store the coordinates of each data point
in its corresponding face node, we don’t need extra
fields in the edge nodes for their parameters.

We can also consider representing the Delaunay
diagram, instead of the Voronoi, by the data
structure of section 5.???D. Actually, the two approaches
give essentially the same data structure! The only
difference will be in the names of nodes and fields:
the vertex nodes of one are the face nodes of the other,
the $\prg{FNext}$ and $\prg{VPrev}$ fields are interchanged,
$\prg{Dest}$ replaces $\prg{RFace}$, and so forth.

To reduce confusion, we will stick to the
Voronoi interpretation of the data structure
(i.e., faces=voronoi regions, vertices=Voronoi vertices).
Since the vertex nodes contain no useful information,
we will do away with them. We will have
therefore one node per data point (with fields $\prg{Edge}$,
$\prg{X}$ and $\prg{Y}$),
and two records per edge (with fields $\prg{RFace}$,
$\prg{FNext}$, $\prg{VPrev}$, and $\prg{Sym}$).\par

As we have seen, an internal face of the Delaunay
diagram for $P$
will have more than three sides only
if four points of $P$ lie on the same circle.
When the data points are reasonably ``random’’,
the second case will hold with probability zero.
Furthermore, when the computations are carried out using
floating-point arithmetic, rounding errors
make it pointless to distinguish between the two cases.\par

In practice, therefore, it may be convenient
to reduce all faces of $\Dscr(P)$ to triangles,
by adding three fictitious data points
$\pi1̌, \pi2̌, \pi3̌$ at infinity, 120\deg apart,
and triangulating all internal faces with more
than three edges. This latter
step is equivalent to replacing every Voronoi vertex
of degree $m>3$ by $m-3$ vertices connected by dummy
edges of length zero:\fig4
The resulting diagram is called the {\it Delaunay
triangulation} for $P$, $\Dscr*(P)$; the Voronoi
diagram corresponding to it will be denoted by $\Vscr*(P)$.
Although this ``patching up’’ of the Voronoi diagram
is slightly ambiguous and mathematically ``ugly’’,
it allow us to eliminate the $\prg{Vprev}$
field of edge records, as we observed in section 5.???D.

\digress{It was once thought that $\Dscr*(P)$ was the
triangulation minimizing the total edge length, and
indeed several proofs of this ``fact’’ were published,
until a counter-example was found.}

\sect {8.7. Computing the Voronoi diagram.}

To compute the edges and vertices of a Voronoi region
$Vp̌$, we can start with $Vp̌\gets \Sscr$ and then perform
$Vp̌\gets Vp̌\inte C(p,q)$ for all $q$ in $P-\set{p}$. After
each step, $Vp̌$ will be essentially a convex polygon, except
that some of its vertices may be at infinity. We can represent
such a region by the list of its vertices, in order of polar
angle around the point $p$, and use the algorithms of chapter
4 to compute each intersection.\par

With this data structure, each intersection can be found
in $O(\log n)$ time, so in $O(n\log n)$ we can compute a
single region, and $O(n↑2\log n)$ the whole diagram. In the
worst case, the running time may be actually $\Theta(n↑2\log
n)$. We clearly need better algorithms.\par

Another straightforward algorithm is an incremental one: starting with the Voronoi
diagram of a single point, we add the remaining points one
by one, updating the diagram after each insertion.
At each step, we have a Voronoi diagram $\Vscr=\Vscr(P)$
and a new data point $p\not\in P$; the goal is to
construct from them
the updated Voronoi $\Vscr\pr$ corresponding
to the augmented set $P\pr=P\un \set{p}$.
For each data point $q$, let’s denote by $Vq̌$ its region
in $\Vscr$, and by $Vq̌\pr$ that in $\Vscr\pr$.\par

The diagrams $\Vscr\pr$ and $\Vscr=$ will differ
only in the neighborhood of $p$; more precisely,
the region $Vq̌\pr$ will be equal to $Vq̌$
minus the interior of $Vp̌\pr$. From this, it is not
hard to see that $Vq̌\pr=Vq̌$, except for the points $q$ that
will be neighbors of $p$ in the dual $\Dscr\pr$
of $\Vscr\pr$..\fig4\par
Therefore, the addition of $p$ requires that we compute
the new region $Vp̌\pr$, and erase that
portion of $\Vscr$ that lies inside it.\par

The hardest part of this computation is figuring out which
points of $P$ will be adjacent to $p$ in $\Dscr\pr$.
It is not hard to get one of them: if we locate the
region $Vq̌$ of the original diagram that contains the
new point $p$, then $q$ will be the nearest neighbor of
$p$ in $P\pr$, and, by lemma 8.11, $p$ and $q$ will be
adjacent in $\Dscr(P\pr)$. More than that, the
midpoint $o$ of $pq$ will be a point of the edge
separating $Vp̌\pr$ from $Vq̌\pr$.\par

It is easy to see that, in the original diagram,
$o$ belongs to $Vq̌$. If we proceed along the
bissector $B(p,q)$, starting from $o$ and moving
counterclockwise as seen from $p$, we will eventually meet
the boundary of $Vq̌$ at some point $v$.
Let $Vř$ be the region of $\Vscr$ that is
on the other side of the boundary at that point, i.e.
let $e{̌qr}$ be the edge of $Vq̌$ containing $v$.\fig4
The point $v$ is equidistant from $p$, $q$ and $r$, and no
other data point is closer to it than those three; therefore,
$v$ is a Voronoi vertex of the new diagram.
The edge $e{̌qr}$ is divided in two pieces by $v$:
the one closest to $p$ has to be discarded,
but the other will be an edge of
$\Vscr\pr$.\par

The third edge of $\Vscr\pr$ incident at $v$ is on
the bissector $B(p,r)$; if we start at $v$ and proceed
counterclockwise (as seen from $p$)
along that bissector, we will eventually leave
$Vř$ at some point $v\pr$, wich will be the
next vertex of $Vp̌\pr$. If we repeat this process, we will
obtain all edges and vertices of $Vp̌\pr$, and eventually
we will get back to the point $o$.\par

At some iteration of this algorithm,
it may happen that the edge we are following
leaves the current region $Vq̌$ through some
vertex $v$ of the current Voronoi diagram,
rather than across an edge.\fig4
In that case, there will be two or more regions
of $\Vscr$ abutting
at $v$, besides the one from which we came. To continue
the traversal, we must move into the region $Vř$ that
is immediately adjacent to $Vq̌↑P$,
in {\it clockwise} order as seen from $v$.\par

One way to avoid this situation is to ensure that
no four data points lie on the same circle. Another
way is to allow edges of zero length (see the
discussion of the `` patched’’ Voronoi diagram,
in section 8.???): if we hit the boundary
of $Vq̌$ at a vertex
$v$, we take either of the edges incident at
$v$, and proceed as if we had hit only that edge.\par

{\bf This is not as simple as it looks. How do
we make sure the resulting graph (abstracting from edge
lengths) is still planar and properly embedded?
I.e., how do we keep the clockwise orderings
of edges around adjacent vertices
consistent if the edges have length zero?
If you throw in numerical errors, the problem
becomes frightening; a satisfactory (but not
easy) compromise would be an algorihm that
works in all cases if exact rational arithmetic
is used throughout. Maybe such a disclaimer
should be added to the introduction.}\par

The algorithm below describes these ideas in more detail.
We are assuming the Voronoi diagram $\Vscr$
is represented by the data structure discussed in
section 8.???.\par

\algorithm{8.1}{Incremental computation of the
Delaunay diagram.}{
\step 1. [Locate nearest neighbor.] {Find the point $q\in P$
that is closest to $p$.}
\step 2. [Initialize.] {Let $a$ be the null edge
record $\Lambda$. Search the list of edges incident
to $q$ in counterclockwise order,
and let $b$ point to the one that would
follow immediately the edge $(q,p)$ (if $q$ has degree 0, i.e.
it is the only point in $P$, set $b\gets \lambda$).}
\step 3. [Insert new edge.] {Create a new record $e$
for the directed edge $(p,q)$, and insert it before $a$, and
its symmetrical $(q,p)$ after $b$. Otherwise,
let $e$ be the first edge incident to $q$, in clockwise order,
that .}}

{\bf ??? proof (using the quadrilateral lemma?}

{\bf ??? Analysis: average $O(n\log n)$?}

{\bf ??? Conjecture: if points are sorted by $x$, the WORST case is $O(n\log n)$.}

{\bf ??? Voronoi by random split}

{\bf ??? $k$-point Voronoi (sorted and unsorted)}

{\bf ??? $R$-$S$ Voronoi (include in Voronoi algorithms?)}

{\bf ??? Line voronoi}

\sect {8.???. Closest point by ``Bucket sort’’.}

An algorithm for the closest-point problem that is entirely
different from the ones based in the\par

The algorithm assumes the points to be uniformly distributed
in some rectangular box of $\Sscr$

\vfill\eject\end