Optimal Parallel Algorithms for Finding Proximate Points, with Applications.
AbstractConsider a set P of points in the plane sorted by x-coordinate. A point p in P is said to be a proximate point if there exists a point q on the x-axis such that p is the closest point to q over all points in P. The proximate point problem is to determine all the proximate points in P. Our main contribution is to propose optimal parallel algorithms for solving instances of size n of the proximate points problem. We begin by developing a work-time optimal algorithm running in O(log log n) time and using ${{n \over {\log \log n}}}$ Common-CRCW processors. We then go on to show that this algorithm can be implemented to run in O(log n) time using ${{n \over {\log n}}}$ EREW processors. In addition to being work-time optimal, our EREW algorithm turns out to also be time-optimal. Our second main contribution is to show that the proximate points problem finds interesting, and quite unexpected, applications to digital geometry and image processing. As a first application, we present a work-time optimal parallel algorithm for finding the convex hull of a set of n points in the plane sorted by x-coordinate; this algorithm runs in O(log log n) time using ${{n \over {\log \log n}}}$ Common-CRCW processors. We then show that this algorithm can be implemented to run in O(log n) time using ${{n \over {\log n}}}$ EREW processors. Next, we show that the proximate points algorithms afford us work-time optimal (resp. time-optimal) parallel algorithms for various fundamental digital geometry and image processing problems. Specifically, we show that the Voronoi map, the Euclidean distance map, the maximal empty circles, the largest empty circles, and other related problems involving a binary image of size nn can be solved in O(log log n) time using$${{{n^2} \over {\log \log n}}}$$Common-CRCW processors or in O(log n) time using ${{{n^2} \over {\log n}}}$ EREW processors.
Introduction
Consider a parallel algorithm that solves an instance of size n of some problem in T p (n) time
using p processors. Traditionally, the main complexity measure for assessing the performance of
Work supported in part by NSF grant CCR-9522093, by ONR grant N00014-97-1-0526, and by Grant-in-Aid for
Encouragement of Young Scientists (08780265) from Ministry of Education, Science, Sports, and Culture of Japan
y Dept. of Electrical and Computer Engineering, Nagoya Institute of Technology, Showa-ku, Nagoya 466, JAPAN,
z Department of Computer Science, Old Dominion University, Norfolk, Virginia 23529, USA, olariu@cs.odu.edu
the algorithm is the amount W (n) of work performed by the algorithm, defined as the product
(n). The algorithm is termed work-optimal if W (n) 2 \Theta(T   (n)), where T   (n) is the running
time of the fastest sequential algorithm for the problem. The algorithm is work-time optimal [20] if
it is work-optimal and, in addition, its running time T p (n) is best possible among the work-optimal
algorithms in that model. Needless to say that one of the challenges of parallel algorithm design is
to produce not only work-optimal but, indeed, whenever possible, work-time optimal algorithms.
Occasionally, an even stronger complexity metric is being used - the so-called time-optimality.
Specifically, an algorithm is time-optimal in a given model, if the problem cannot be solved faster
in that model, even if an unbounded number of processors were available.
In this paper we assume the Parallel Random Access Machine (PRAM, for short) which consists
of synchronous processors, each having access to a common memory. We refer the interested reader
to [20] for an excellent discussion of the PRAM model.
Let P be a set of points in the plane sorted by x-coordinate. A point p is a proximate point
of P if there exists a point on the x-axis closer to p than to any other point in P . The proximate
points problem asks to determine all proximate points in P . Clearly, the proximate points problem
can be solved, using an algorithm for finding the Voronoi diagram. However, as argued in [16],
the computation of the Voronoi diagram
log n) time even if the n points are sorted
by x-coordinate. Thus, this naive approach does not yield an optimal solution to the proximate
points problem. Recently, Breu et al. [6] proposed a linear-time algorithm for the proximate points
problem. In spite of its optimality, the algorithm of Breu et al. [6] relies in crucial ways on stack
operations, notoriously hard to parallelize.
Our first main contribution is to propose parallel algorithms for solving instances of size n of
the proximate points problem. Specifically, we first exhibit an algorithm running in O(log log n)
time using n
log log n Common-CRCW processors. We then go on to show that this algorithm can be
implemented to run in O(log n) time using n
log n EREW processors. Our Common-CRCW algorithm
is work-time optimal; the EREW algorithm turns out to also be time-optimal. We establish the
work-time optimality of our Common-CRCW algorithm by a reduction from the minimum finding
problem; the time-optimality of our EREW algorithm follows by a reduction from the OR problem.
Our second main contribution is to show that the proximate points problem has interesting,
and quite unexpected, applications to digital geometry and image processing. To begin, we present
a work-time optimal parallel algorithm for computing the convex hull of a set of n points in the
plane sorted by x-coordinate. This algorithm runs in O(log log n) time using n
log log n Common-
CRCW processors or in O(log n) time using n
log n EREW processors. We show that this algorithm
is work-time optimal in the CRCW model and, in addition, time-optimal in the EREW.
Numerous parallel algorithms have been proposed for computing the convex hull of sorted points
in the plane [4, 7, 12, 13, 21]. Recently, Chen [7] presented an O(log n)-time algorithm using n
log n
processors. Chen et al. [12] presented work-optimal algorithms running in O(log n)-time
algorithm and using n
log n EREW processors, and in an O(log log n)-time algorithm using n
log log n
Common-CRCW processors. Quite recently, Berkman et al. [3] presented an O(log log n)-time
algorithm using n
log log n Common-CRCW processors. Our algorithm features the same performance
as those in [3, 7, 12]. However, our algorithm is much simpler and more intuitive. Further, to the
best of our knowledge, the work-time optimality of the CRCW version and the time-optimality of
the EREW version algorithm has not been solved yet.
Given a binary image the Voronoi map assigns to each pixel in the image the position of the
nearest black pixel. The Euclidean distance map assigns to each pixel the Euclidean distance to the
nearest black pixel. An empty circle of the image is a circle whose interior contains only white pixels.
A maximal empty circle is an empty circle contained in no other empty circle. A largest empty circle
is an empty circle of the largest radius. We refer the reader to Figure 1 for an illustration. The
largest square, diamond, n-gon, etc. are defined similarly. These computations are known to have
numerous applications ranging from clustering and shape analysis [2, 17] to handoff management
in cellular systems [26] to image compression, decomposition, and reconstruction [5, 23, 27, 28, 31].
As further applications, we propose algorithms for computing the Voronoi map, the Euclidean
distance map, the maximal empty circles, and the largest empty circles of a binary image of size
n \Theta n. We begin by presenting a work-time optimal algorithm that computes the Voronoi map
and the Euclidean distance map of a binary image of size n \Theta n in O(log log n) time using n 2
log log n
Common-CRCW processors or in O(log n) time using n 2
log n EREW processors. We also show that
the distance map for various metrics including the well known L k metrics, (k - 1), can also be
computed in the same manner. We then go on to show that all the maximal empty circles and
a largest empty circle of an n \Theta n binary image can be found in O(log log n) time using n 2
log log n
Common-CRCW processors or in O(log n) time using n 2
log n EREW processors. As it turns out,
with minimal changes, this algorithm is applicable to various other kinds of empty figures including
squares, diamonds, n-gon etc.
Recently, Chen et al. [8, 11] and Breu et al. [6] presented O(n 2 )-time sequential algorithms for
computing the Euclidean distance map. Roughly at the same time, Hirata [19] presented a simpler
sequential algorithm to compute the distance map for various distance metrics including
Euclidean, 4-neighbor, 8-neighbor, chamfer, and octagonal. A number of parallel algorithms for
computing the Euclidean distance map have been developed for various parallel models [1,9,10,14].
In particular, the following results have been reported in the recent literature. Lee et al. [22]
Figure

1: Illustrating the Euclidean distance map and the largest empty circle
presented an O(log 2 n)-time algorithm using n 2 EREW processors. Pavel and Akl [24] presented an
algorithm running in O(log n) time and using n 2 EREW processors. Clearly, these two algorithms
are not work-optimal. Chen [8] presented a work-optimal O( n 2
)-time algorithm using p, (p log p -
n), EREW processors. This yields an O(n log n)-time algorithm using n
log n EREW processors.
Fujiwara et al. [18] presented a work-optimal algorithm running in O(log n) time and using n 2
log n
EREW processors and in O( log n
log log n ) time using n 2 log log n
log n Common-CRCW processors. Although
Fujiwara et al. [18] claim that their algorithm is applicable to various distance maps, a closer
analysis reveals that it only applies to a few distance metrics. The main problem seems to be
that their algorithm uses a geometric transform that depends in a crucial way on properties of the
Euclidean distance and, therefore, does not seem to generalize. As we see it, our Euclidean distance
map algorithm has three major advantages over Fujiwara's algorithm. First, the performance of
our algorithm for the CRCW is superior; second, our algorithm applies to a large array of distance
finally, our algorithm is much simpler and more intuitive.
The remainder of this paper is organized as follows: Section 2 introduces the proximate points
problem for the Euclidean distance metric and discusses a number of technicalities that will be
crucial ingredients in our subsequent algorithms. Section 3 presents our parallel algorithms for the
Common-CRCW and the EREW. Section 4 proves that these algorithms are work-time, respectively
time-optimal. Section 5 presents a work-time optimal parallel algorithm for computing the convex
hull of sorted points in the plane. Section 6 uses the proximate points algorithm to computing the
Voronoi map, the Euclidean distance map, the maximal empty circles, and the largest empty circles
of a binary image. Section 7 offers concluding remarks and open problems. Finally, the Appendix
discusses other distance metrics to which the algorithms presented in Section 3 apply.
2 The proximate points problem: a first look
In this section we introduce the proximate points problem along with a number of geometric results
that will lay the foundation of our subsequent algorithms. Throughout, we assume that a point p
is represented by its Cartesian coordinates (x(p); y(p)). As usual, we denote the Euclidean distance
between the planar points p and q by d(p;
Consider a collection of n points sorted by x-coordinate, that is, such that
We assume, without loss of generality that all the points in P have
distinct x-coordinates and that all of them lie above the x-axis. The reader should have no difficulty
confirming that these assumptions are made for convenience only and do not impact the complexity
of our algorithms.
Recall that for every point p i of P the locus of all the points in the plane that are closer to
than to any other point in P is referred to as the Voronoi polygon associated with p i and is
denoted by V (i). The collection of all the Voronoi polygons of points in P partitions the plane into
the Voronoi diagram of P (see [25] p. 204). Let I i , (1 - i - n), be the locus of all the points q
on the x-axis for which d(q; In other words, q 2 I i if and
only if q belongs to the intersection of the x-axis with V (i), as illustrated in Figure 2. In turn, this
implies that I i must be an interval on the x-axis and that some of the intervals I i ,
may be empty. A point p i of P is termed a proximate point whenever the interval I i is nonempty.
Thus, the Voronoi diagram of P partitions the x-axis into proximate intervals. Since the points of
are sorted by x-coordinate, the corresponding proximate intervals are ordered, left to right, as
point q on the x-axis is said to be a boundary point between p i and p j if q is
equidistant to p i and p j , that is, d(p It should be clear that p is a boundary point
between proximate points p i and p j if and only if the q is the intersection of the (closed) intervals
I i and I j . To summarize the previous discussion we state the following result.
Proposition 2.1 The following statements are satisfied:
ffl Each I i is an interval on the x-axis;
ffl The intervals I 1 ; I lie on x-axis in this order, that is, for any non-empty I i and I j
lies to the left of I
I 1 I 2 I 4 I 6 I 7

Figure

2: Illustrating proximate intervals
I 1 I 2 I 3 I 4
I 0

Figure

3: Illustrating the addition of p to g.
ffl If the non-empty proximate intervals I i and I j are adjacent, then the boundary point between
Referring again to Figure 2, among the seven points, five points are proximate
points, while the others are not. Note that the leftmost point p 1 and the rightmost point p n are
always proximate points.
Given three points we say that dominated by p i and p k whenever
fails to be a proximate point of the set consisting of these three points. Clearly, p j is dominated
by p i and p k if the boundary of p i and p j is to the right of that of p j and p k . Since the boundary
of any two points can be computed in O(1) time, the task of deciding for every triple (p
whether p j is dominated by p i and p k takes O(1) time using a single processor.
Consider a collection of points in the plane sorted by x-coordinate, and a
point p to the right of P , that is, such that x(p 1 x(p). We are interested
in updating the proximate intervals of P to reflect the addition of p to P as illustrated in Figure 3.
We assume, without loss of generality, that all points in P are proximate points and let
I n be the corresponding proximate intervals. Further, let I 0
p be the up-dated
proximate intervals of P [ fpg. Let p i be a point such that I 0
i and I 0
are adjacent. By (iii) in
Proposition 2.1, the boundary point between p i and p separates I 0
i and I 0
. As a consequence, (ii)
implies that all the proximate intervals I 0
n must be empty. Furthermore, the addition of p
to P does not affect any of the proximate intervals I j , 1 In other words, for all 1
I 0
are empty, the points p are dominated by p i and p. Thus, every
point n), is dominated by otherwise, the boundary between
would be to the left of that of that between p j and p. This would imply that the non-empty
interval between these two boundaries corresponds to I 0
j , a contradiction. To summarize, we have
the following result.
Lemma 2.2 There exists a unique point p i of P such that:
ffl The only proximate points of P [ fpg are
ffl For the point p j is not dominated by
I 0
ffl For dominated by and the interval I 0
j is empty.
i and I 0
are consecutive on the x-axis and are separated by the boundary point between
and p,
be a collection of proximate points sorted by x-coordinate and let p
be a point to the left of P , that is, such that x(p) ! x(p 1
reference we now take note of the following companion result to Lemma 2.2. The proof is identical
and, thus, omitted.
Lemma 2.3 There exists a unique point p i of P such that:
ffl The only proximate points of P [ fpg are
ffl For not dominated by p and p j+1 . Moreover, for
I 0
ffl For the point p j is dominated by p and p j+1 and the interval I 0
j is empty.
p and I 0
are consecutive on the x-axis and are separated by the boundary point between p and
The unique point p i whose existence is guaranteed by Lemma 2.2 is termed the contact point
between P and p. The second statement of Lemma 2.2 suggests that the task of determining the
unique contact point between P and a point p to the right or left of P reduces, essentially, to binary
search.
Now, suppose that the set
into two subsets g. We are interested
in updating the proximate intervals in the process or merging PL and PR . For this purpose, let
I 2n be the proximate intervals of PL and PR , respectively. We as-
sume, without loss of generality, that all these proximate intervals are nonempty. Let I 0
be the proximate intervals of . We are now in a position to state and prove the next
result which turns out to be a key ingredient in our algorithms.
Lemma 2.4 There exists a unique pair of proximate points PR such that
ffl The only proximate points in PL [ PR are
are empty, and I 0
ffl The proximity intervals I 0
i and I 0
are consecutive and are separated by the boundary point
between
Proof. Let i be the smallest subscript for which p i 2 PL is the contact point between PL and a
point in PR . Similarly, let j be the largest subscript for which the point PR is the contact
point between PR and some point in PL . Clearly, no point in PL to the left of p i can be a proximate
point of P . Likewise, no point in PR to the left of p j can be a proximate point of P .
Finally, by Lemma 2.2 every point in PL to the left of p i must be a proximate point of P .
Similarly, by Lemma 2.3 every point in PR to the right of p i must be a proximate point of P , and
the proof of the lemma is complete.
The points p i and p j whose existence is guaranteed by Theorem 2.4 are termed the contact
points between PL and PR . We refer the reader to Figure 4 for an illustration. Here, the contact
points between PL and PR are p 4 and p 8 .
Next, we discuss a geometric property that enables the computation of the contact points p i
and p j between PL and PR . For each point p k of PL , let q k denote the contact point between p k
and PR as specified by Lemma 2.3. We have the following result.
Lemma 2.5 The point p k is not dominated by p k\Gamma1 and q k if 2 - k - i, and dominated otherwise.
I 1 I 2 I 3 I 4 I 5
I 6 I 7 I 8 I 9 I 10
I 0
9 I 0
Figure

4: Illustrating the contact points between two sets of points
Proof. If dominated by p k\Gamma1 and q k , then I 0
k must be empty. Thus, Lemma 2.4
guarantees that p k , (2 - k - i), is not dominated by p k\Gamma1 and q k . Suppose that p k , (i
is not dominated by p k\Gamma1 and q k . Then, the boundary point between p k and q k is to the right of
that between p k\Gamma1 and p k . Thus, the non-empty interval between these two boundaries corresponds
to I 0
k , a contradiction. Therefore, p k , (i n), is dominated by p k\Gamma1 and q k , completing the
proof.
Lemma 2.5 suggests a simple, binary search-like, approach to finding the contact points p i
and between two sets PL and PR . In fact, using a similar idea, Breu et al. [6] proposed a
sequential algorithm that computes the proximate points of an n-point planar set in O(n) time.
The algorithm in [6] uses a stack to store the proximate points found and, consequently, seems very
hard to parallelize.
3 Parallel algorithms for the proximate points problem
We begin by discussing a parallel algorithm for solving the proximate points problem on the
Common-CRCW. The algorithm will then be converted to run on the EREW. We rely, in part, on
the solution to the well-known LEFTMOST-ONE problem: given a sequence b 1
determine the smallest i, (1 - i - n), such that b
Lemma 3.1 [20] An instance of size n of the LEFTMOST-ONE problem can be solved in O(1)
time using n Common-CRCW processors.
Consider a set of points such that x(p 1 To capture
the neighboring proximate points of each point use three indices c i , l i and r i
defined as follows:
non-proximate points
proximate points

Figure

5: Illustrating indices l i , c i , and r i for a point p i
1. is an proximate pointg;
2. l is an proximate pointg;
3. r is an proximate pointg.
We refer the reader to Figure 5 for an illustration. Note that we must have l
there is no proximate point p j such that l i
then c
Next, we are interested in finding the contact point between the set and a
new point p with x(p n We assume that for every i, (1 - i - n), c i , l i , and r i are available,
and that m, (m - n), processors are at our disposal. The algorithm is essentially performing m-ary
search using Lemma 2.2.
Algorithm Find-Contact-Point(P; p)
Extract a sample S(P ) of size m consisting of the points p c
in P . For
every k, (k - 0), check whether the point p c k n
is dominated by p l k n
and p, and whether
is dominated by p c k n
and p. If p c k n
is not dominated but p r k n
is dominated,
then
is the desired contact point.
such that the point p r k n
is not dominated by p c k n
and p, and p c (k+1) n
is dominated by p l (k+1) n
and p.
Step 3 Execute recursively this algorithm for the set of points P
g to find the contact point.
1, the set P 0 contains at most l (k+1) n
points. Hence, the depth of the recursion is O( log n
log m ). Notice, further, that
algorithm Find-Contact-Point does not perform concurrent reading or writing. Thus, we have
the following result.
Lemma 3.2 Given a set of n points in the plane sorted by x-coordinate and a
point p, with x(p n the task of finding the contact point between P and p can be performed
in O( log n
log using m EREW processors.
Next, consider two sets of points in the
plane such that x(p 1 Assume that for every i the indices c i , l i and r i are
given and that m processors are available to us. The following algorithm finds the contact points
of PL and PR by
m-ary search using Lemma 2.5.
Algorithm Find-Contact-Points-Between-Sets(PL ; PR )
sample points S(PL
from PL . By using the algorithm
Find-Contact-Point and p m of the processors available each, determine for each
sample point
1), the corresponding contact point q c k n
in PR .
Step 2 For each k, (0 - k -
check whether the point p c k n
is dominated by p l k n
and q c k n
, and whether the point p r k n
is dominated p c k n
and q c k n
. If p c k n
is not dominated, yet p r k n
is, output
and q c k n
as the desired contact points.
Step 3 Find k such that the point p r k n
is not dominated by p c k n
and q c k n
is dominated by p l (k+1) n
and q c (k+1) n
Step 4 Execute recursively algorithm Find-Contact-Points-Between-Sets for the sets P 0
and PR and return the desired contact points.
It is not hard to see that algorithm Find-Contact-Points-Between-Sets involves concurrent
reads (because several processors may access a point concurrently), but does not involve concurrent
operations. By Lemma 3.2, Step 1 can be takes O( log n
log m ) time on the CREW model. Steps 2
and 3 run, clearly, in O(1) time. Since P 0
L contains at most n
the depth of recursion
is O( log n
log m ). Thus, altogether, algorithm Find-Contact-Points-Between-Sets runs in O( log 2 n
using m CREW processors.
Lemma 3.3 Given the sets of points in the
plane such that x(p 1 the task of finding the contact points between PL and
PR can be performed in O( log 2 n
using m CREW processors.
Next, we are interested in designing an algorithm to compute the proximate points of a set P
on n points in the plane sorted by x-coordinate in O(log log n) time on the Common-CRCW. We
assume that n processors are available to us. We begin by determining for every i, the indices c i ,
l i , and r i . With this information available, all that remains to be done is to retain all the points p i
for which c i. The details follow.
Algorithm Find-Proximate-Points(P)
Partition the set P into n 1=3 subsets such that for every k, (0 - k -
g. For every point p i in P k , (0 - k - n 1=3 \Gamma 1),
determine the indices c i , l i , and r i local to P k .
Compute the contact points of each pair of sets P i and using
n 1=3 of the processors available. Let q i;j 2 P i denote the contact point between P i and P j .
Step 3 For every P i , find the rightmost contact point p rc i
among all the points q i;j with
and find the leftmost contact point p lc i
over all points q i;j with j ? i. Clearly, x(p rc i
ig.
Step 4 For each set P i , the proximate points lying between rc i and lc i (inclusive) are proximate
points of P . Update each c
It is clear that Step 2 can be performed in O( log n 2=3
log runs in O(1)
time using Lemma 3.1. At this moment, the reader may wonder how the updating of the indices
can be performed efficiently. In fact, as it turns out, this update can be done
in O(1) time. Since the task of updating l i and r i is, essentially, the same as that of updating c i ,
we will only focus on c i . In each P i , for all the points p j , (rc the value of c j is not
changed. For all the points p j with lc i ! j, the value of c j must be changed to lc i . For all points p j
with the value of c j is changed to lc has an proximate point. However, if P
contains no proximate points, we have to find the nearest subset that contains a proximate point.
To do this, first check whether each P i has a proximate point using n 2=3 processors each. Thus,
totally, processors are used for this task. Next, using Lemma 3.1 we determine
k such that contains a proximate pointg for each P i . Since P has n 1=3
groups, this task can be done in O(1) time and n 1=3 processors each and, totally, n 1=3 \Delta n
processors are used. Thus, Step 4 can be done in O(1) time using n processors.
Let TCRCW (n) be the running time of this algorithm. To find the recurrence describing the worst
case running time of algorithm Find-Proximate-Points, we note that Step 1 executes recursively
this algorithm for n 2=3 points, while Steps 2, 3, and 4 run in O(1) time. Thus, we have
confirming that TCRCW (n) 2 O(log log n). Thus, we have:
Lemma 3.4 An instance of size n of the proximate points problem can be solved in O(log log n)
time using n Common-CRCW processors.
Next, we show that the number of processors can be reduced by a factor of log log n without increasing
the running time. The idea is as follows: begin by partitioning the set P into n
log log n subsets
log log n
each of size log log n. Next, using algorithm Sequential-Proximate-Points
find the proximate points within each subset in O(log log n) sequential time and, in the process,
remove from P all the points that are not proximate points. For every i, (1 - i - n
log log n ), let
proximate points in the set P i .
At this moment, execute algorithm Find-Proximate-Point on P 1
log log n
. Since n
processors are required in order to update the indices c i , l i , and r i in O(1), we will proceed slightly
differently. The idea is the following: while executing the algorithm, some of the (currently)
proximate points will cease to be proximate points. To maintain this information efficiently, we
use ranges
log log n
log log n
such that for each P i , fp i;L
are the current proximate points. While executing the algorithm, P i may contain no proximate
points. To find the neighboring proximate points, we use the pointers L 0
log log n
and
log log n
such that
and the set P j contains a proximate pointg,
and the set P j contains a proximate pointg.
By using this strategy, we can find the contact point between a point and P in O( log n
log using
processors as discussed in Lemma 3.2. Thus, the contact points between two subsets can be
found in the same manner as in Lemma 3.3. Finally, the algorithm for Lemma 3.1 can update
i in Step 4 in O(1) time by using O( n
log log n ) processors. To summarize, we have
the following result.
Theorem 3.5 An instance of size n of the proximate points problem can be solved in O(log log n)
time using n
log log n Common-CRCW processors.
We close this section by pointing out that algorithm Find-Proximate-Points can be implemented
efficiently on the EREW. For this purpose, we rely, in part, on the following well known
result [20].
Lemma 3.6 A single step execution of the m-processor CRCW can be simulated by an m-processor
EREW in O(log m) time.
By Lemma 3.6, Steps 2, 3, and 4 of the algorithm can be performed in O(log n) time using
processors, as the CRCW performs these steps in O(1) time using n processors. Let
TEREW (n) be the worst-case running time on the EREW. Then, the recurrence describing the
confirming that T (n) 2 O(log n). Consequently, we have:
Lemma 3.7 An instance of size n of the proximate points problem can be solved in O(log n) time
using n EREW processors.
Using, essentially, the same idea as for the Common-CRCW, we can reduce the number of
processors by a factor of log n without increasing the computing time. Specifically, in case of the
EREW, the n points are partitioned into n
log n subsets each of size log n. Thus, we have
Theorem 3.8 An instance of size n of the proximate points problem can be solved in O(log n) time
using n
log n EREW processors.
4 Lower Bounds
The main goal of this section is to show that the running time of the Common-CRCW algorithm
for the proximate points problem developed in Section 3 cannot be improved while retaining work-
optimality. This, in effect, will prove that our Common-CRCW algorithm is work-time optimal.
We then show that our EREW algorithm is time-optimal.
The work-optimality of both algorithms is obvious; in order to solve the proximate points
problem every point must be accessed at least once.
n) work is required of any algorithm
solving the problem.
Our lower bound arguments rely, in part, on the following fundamental result of Valiant [30].
Lemma 4.1 The task of finding the minimum (maximum) of n real numbers
log n)
time on the CRCW provided that n log O(1) n processors are available.
We now show that the lower bound of Lemma 4.1 holds even if all the item are non-negative.
Lemma 4.2 The task of finding the minimum (maximum) of n non-negative (non-positive) real
numbers requires
log n) time on the CRCW provided that n log O(1) n processors are available.
Proof. Assume that the minimum (maximum) of n non-negative numbers can be computed in
o(log log n) time using n log O(1) n CRCW processors.
With this assumption, we can find the minimum of n real numbers in o(log log n) time as follows:
first, in O(1) time, check whether there are negative numbers in the input. If not, the minimum
of input items can be computed in o(log log n) time. If negative numbers exist, replace every non-positive
number by 0 and find the maximum of their absolute values in the resulting sequence in
o(log log n) time. The maximum thus computed corresponds to the minimum of the original input.
Thus, the minimum of n real numbers can be computed in o(log log n) time, contradicting Lemma
4.1.
Further, we rely on the following classic result of Cook et al. [15].
Lemma 4.3 The task of finding the minimum (maximum) of n real numbers
on the CREW (therefore, also on the EREW) even if infinitely many processors are available.
We shall reduce the task of finding the minimum of a collection A of n non-negative a 1 ; a
to the proximate points problem. Our plan is to show that an instance of size n of the problem of
finding the minimum of a collection of non-negative numbers can be converted, in O(1) time, to an
instance of size 2n of the the proximate points problem involving sorted points in the plane.
For this purpose, let be a set of arbitrary non-negative real numbers that
are input to the minimum problem. We construct a set of points in the plane
by setting for every i, (1
Notice that this construction guarantees that the points in P are sorted by x-coordinate and that
for every i, (1 - i - n), the distance between the point p i and the origin is exactly
Intuitively, our construction places the 2n points circles centered at
the origin. More precisely, for every i, (1 - i - n), the points p i and p n+i are placed on such a
circle C i with radius
. It is very important to note that the construction above can be
carried out in O(1) time using n EREW processors.
In our subsequent arguments, we find it convenient to rely on the next technical result.
Lemma 4.4 Both p i and p i+n are proximate points if and only if a i is the minimum of A.
Proof. Let a i be the minimum of A and refer to Figure 6. Clearly, C i is the circle of smallest
radius containing p i and p i+n , while all the other points lie outside C i . Hence, p i and p i+n are the
closest points of P from the origin. Thus, the boundary between
lies to the left of the origin: were this not true, p j would be closer to the origin than p i+n . The
following simple facts are proved in essentially the same way.
a i
a i+n
I i+n
I i
O

Figure

Illustrating P for Lemma 4.4
1. The boundary point between p j and p i+n lies to the left of the origin if 1 and to
the right if i
2. The boundary point between p j and p i lies to the left of the origin if 1 and to the
right if
for each point n), the boundary between p j and p i lies to the left of that
between p j and p n+i , p j is not proximate point. Thus, for j 6= i, either p j or p j+n fails to be a
proximate point. Further, for the point p i the boundary with lies to the left of the
origin, and that with n) lies to the right of the origin (or is the origin itself). Thus, p i
is a proximate point. The fact that p i+n is a proximate point follows by a mirror argument. This
completes the proof.
Lemma 4.4 guarantees that the minimum of A can be determined in O(1) time once the proximate
points of P are known. Now, Lemma 4.1 implies the following important result.
Theorem 4.5 Any algorithm that solves an instance of size n of the proximate points problem on
the CRCW must take \Omega\Gammake/ log n) time in the worst case, provided that n log O(1) n processors are
available.
Using exactly the same construction, in combination with Lemma 4.3 we obtain the following lower
bound for the CREW.
Theorem 4.6 Any algorithm that solves an instance of size n of the proximate points problem on
the CREW (also on the EREW) must take \Omega\Gammake/ n) time, even if an infinite number of processors
are available.
Notice that the EREW algorithm for the proximate points problem presented in Section 3
running in O(log n) time using n
log n processors features the same work and time performance on
the CREW-PRAM. By Theorem 4.6 the corresponding CREW algorithm is also time-optimal.
It is straightforward to extend the previous arguments to handle the case of the L k metric.
Specifically, in this case, for every i, (1 - i - n), the points
allow us to find the minimum of A. Thus, Theorems 4.5 and
4.6 provide lower bounds for solving the proximate points problem for the distance metric L k .
5 Computing the convex hull
The main goal of this section is to show that the proximate points algorithms developed in Section 3
yield a work-time optimal (resp. time-optimal) algorithm for computing the convex hull of a set of
points in the plane sorted by x-coordinate. We begin by discussing the details of this algorithm.
In the second subsection we establish its work-time (resp. time) optimality.
5.1 The convex hull algorithm
be a set of n points in the plane with x(p 1
line segment partitions the convex hull of P into the lower hull, lying below the segment, and
the upper hull, lying above it. We focus on the computation of the lower hull only, the computation
of the upper hull being similar.
For a sequence a 1 , a of items, the prefix maxima is the sequence a 1 , maxfa 1 ; a
g. For later reference, we state the following result [20, 29].
Lemma 5.1 The task of computing the prefix maxima (prefix minima) of an n-item sequence can
be performed in O(log n) time using n
log n EREW processors or in O(log log n) time using n
log log n
Common-CRCW processors.
be a set of n points in the plane sorted by x-coordinate as x(p 1
We define a set let of n points by setting
for every i, (1 - i - n), q
It is important to note that the points in Q

Figure

7: Illustrating the proof of Lemma 5.2
are also sorted by x-coordinate. The following surprising result captures the relationship between
the sets P and Q we just defined.
Lemma 5.2 For every j, (1 - j - n), p j is an extreme point of the lower hull of P if and only if
q j is a proximate point of Q.
Proof. If is an extreme point of P and q j is a proximate point of
Q. Thus, the lemma is correct for Now consider an arbitrary j in the range
be arbitrary subscripts such that 1 -
be the boundaries between q i and q j , and between q j and q k , respectively, and refer to Figure 7.
Clearly,
Thus, we have
Similarly, we obtain
It is easy to see that the point q j is not dominated by q i and q k if and only if x(b
Notice that the slopes of the line segments are 2x(b i ) and 2x(b k ), respectively.
Thus, the point p j lies below the segment p i p k if and only if 2x(b Consequently, the
point lies below the segment p i p k if and only if the point q j is not dominated by q i and q k . In
other words, the point p j is an extreme point of the lower hull of P if and only if q j is a proximate
point of Q.
Lemma 5.2 suggests the following algorithm for determining the extreme points of the lower
hull of g.
Algorithm Find-Lower-Hull(P )
Construct the set by setting for every i, (1 - i - n),
Determine the proximate points of Q and report p i as an extreme point of the lower hull
of P whenever q i is a proximate point of Q.
The preprocessing in Step 0 amounts to translating the set P vertically in such a way that for
every This affine transformation does not affect the convex hull
of P . The correctness of this simple algorithm follows directly from Lemma 5.2. To argue for the
running time, we note that by Lemma 5.1 Step 0 takes O(log log n) time and optimal work on the
Common-CRCW or O(log n) time and optimal work on the EREW. Step 1 runs in O(1) time using
optimal work on either the Common-CRCW or the EREW. By Theorems 3.5 and 3.8, Step 2 takes
O(log log n) time and optimal work on the Common-CRCW or O(log n) time and optimal work on
the EREW. Thus, we have proved the following result.
Theorem 5.3 The task of computing the convex hull of a set of n points sorted by x-coordinate
can be performed in O(log log n) time using n
log log n Common-CRCW processors or in O(log n) time
using n
log n EREW processors.
5.2 The optimality of the convex hull algorithm
The main goal of this subsection is to show that the convex hull algorithm described in the previous
subsection is work-time optimal on the Common-CRCW and, in addition, time-optimal on the
CREW and EREW.
Clearly, every point must be read at least once to solve the proximate points problem. Thus,
O(n)-time is required to solve the problem, and our convex hull algorithms (Common-CRCW or
are work-optimal.
Next, we show that given a set of n non-negative integers their maximum
can be determined by using any algorithm for computing the convex hull of a set of sorted points
in the plane. For this purpose, we exhibit an O(1)-time reduction of the maximum problem to the
convex hull problem. The proof technique is similar to the one employed for the proximate points
problem.
With A given construct a set of points in the plane by setting for every i,
is a set of point in
the plane sorted by x-coordinate. The following result relates the sets A and P .
Lemma 5.4 The item a i is the maximum of A if and only if both p i and p i+n are points on the
upper hull of P .
Proof. Let a i be the maximum of A. By construction, both p i and p i+n are points of the upper
hull of P . Further, none of the points p can belong to the upper hull of P .
Thus, there exist no subscript j, (j 6= i), for which both p j and p n+i belong to the upper hull of P .
This completes the proof.
Consequently, to find the maximum of A all we need do is to find an index i such that both p i
and p n+i are points of the upper hull. Therefore, the problem of finding the upper hull of 2n sorted
points in the plane is at least as hard as the problem of finding the maximum of n non-negative
numbers. Thus, we have the following important result.
Theorem 5.5 The task of finding the convex hull of n points
log n) time on the
CRCW, provided that n log O(1) n processors are available.
Similarly, we have the following companion result.
Theorem 5.6 The task of finding the convex hull of n points
n) time on the CREW,
even if infinitely many processors are available.
By Theorems 5.5 and 5.6 the convex hull algorithms developed in the previous section are work-time
optimal. In addition, the EREW algorithm is both work-time and time-optimal.
6 Applications to image processing
A binary image I of size n \Theta n is maintained in an array b i;j , (1 - n). It is customary to refer
to pixel (i; j) as black if b The rows of the image will be numbered
bottom up starting from 1. Likewise, the columns will be numbered left to right, with column 1
being the leftmost. In this notation pixel b 1;1 is in the south-west corner of the image.
The Voronoi map associates with every pixel in I the closest black pixel to it (in the Euclidean
metric). More formally, the Voronoi map of I is a function I such that for every (i; j),
only if
where
is the Euclidean distance between pixels (i;
The Euclidean distance map of image I associates with every pixel in I the Euclidean distance
to the closest black pixel. Formally, the Euclidean distance map is a function R such that
for every (i; j), (1 -
In our subsequent arguments we find it convenient to rely on the solution to the NEAREST-ONE
problem: given a sequence of 0's and 1's, determine the closest 1 to every item
in A. As a direct corollary of Lemma 5.1 we have
Lemma 6.1 An instance of size n of the NEAREST-ONE problem can be solved in O(log log n)
time using n
log log n Common-CRCW processors or in O(log n) time using n
log n EREW processors.
We assume a binary image I of size n \Theta n as discussed above and the availability of n 2
processors, where T log n) for the Common-CRCW and T n) for the EREW.
We now outline the basic idea of our algorithm for computing the Voronoi map and the Euclidean
distance map of image I. We begin by determining, for every pixel in row j, (1 - j - n), the
nearest black pixel, if any, in the same column of the subimage of I. More precisely, with every
pixel (i; j) we associate the value
Next, we construct an instance of the
proximate points problem for every row j, (1 - j - n), in the image I involving the set P j of points
in the plane defined as P
Having solved, in parallel, all these instances of the proximate points problem, we determine, for
every proximate point p i;j in P j its corresponding proximity interval I i . With j fixed, we determine
for every pixel (i; (that we perceive as a point on the x-axis) the identity of the proximity interval
to which it belongs. This allows each pixel (i; j) to determine the identity of the nearest pixel to
it. The same task is executed for all rows in parallel, to determine for every pixel (i;
in row j the nearest black pixel. The details are spelled out in the following algorithm.
Algorithm Voronoi-and-Euclidean-Distance-Map(I)
Step 1 For each pixel (i; j), compute the distances
ng
to the nearest black pixel in the same column as (i; j) in the subimage of I.
Step 2 For every j, (1 ng. Compute the proximate
points E(P j ) of P j .
Step 3 For every point p in E(P j ) determine its proximity interval of P j .
Step 4 For every i, (1 - i - n), determine the proximate intervals of P j to which the point (i;
(corresponding to pixel (i; j)) belongs.
The correctness of this algorithm being easy to see we turn to the complexity. Step 1 can be
performed in O(T (n)) time using the processors available by using Lemma 6.1. Theorem 3.5
and 3.8 guarantee that Step 2 takes O(T using n
processors. By Lemma 6.1, Steps 3
and 4 can be performed in the same complexity. Thus, we have the following important result.
Theorem 6.2 The task of computing the Voronoi map and the Euclidean distance map of a binary
image of size n \Theta n can be performed in O(log log n) time using n 2
log log n Common-CRCW processors
or in O(log n) time using n 2
log n EREW processors.
Recall that an empty circle in the image I is a circle filled with white pixels. The task of
computing the largest empty circles in an image is a recurring theme in pattern recognition, robotics,
and digital geometry [17]. An empty circle is said to be maximal if it is contained in no other empty
circle. An empty circle is said to me maximum if its radius is as large as possible. It is clear that
a maximum empty circle is also a maximal, but not conversely. We now turn to the task of
determining all maximal (resp. maximum) empty circles in an input image I.
Algorithm All-Maximal-Empty-Circles(I)
Compute the Euclidean distance map m of I.
Step 2 For each pixel (i; j), (1 - I compute the smallest distance u
jg to the border of the image. Then, compute r which is
the largest radius of every empty circle centered at the pixel (i; j).
Step 3 For each pixel (i; check whether there exists a neighboring pixel (i
1), such that the circle with radius r i;j and origin (i; j) is included by the circle with
radius r i 0 ;j 0 and origin (i no such circle exists, label the circle of radius r i;j centered at
as a maximal empty circle.
g. Every pixel (i; in I for which r its empty
circle as the largest empty circle of I.
Clearly, all the steps of this simple algorithm can be performed in O(log log n) time using n 2
log log n
Common-CRCW processors or in O(log n) time using n 2
log n EREW processors. Thus we have
Corollary 6.3 The task of labeling all the maximal empty circles and of reporting a maximum
empty circle of a binary image of size n \Theta n can be performed in O(log log n) time using n 2
log log n
Common-CRCW or in O(log n) time using EREW n 2
log n processors.
Conclusions
Our first main contribution is to propose optimal parallel algorithms for solving instances of size
n of the proximate points problem. Our first algorithm runs in O(log log n) time and uses n
log log n
Common-CRCW processors. This algorithm can, in fact, be implemented to run in O(log n) time
using n
log n EREW processors. The Common-CRCW algorithm is work-time optimal; the EREW
algorithm is, in addition, time-optimal. out to also be time-optimal.
Our second main contribution is to show that the proximate points problem finds interesting,
and quite unexpected, applications to digital geometry and image processing. As a first application
we presented a work-time optimal parallel algorithm for finding the convex hull of a set of n points
in the plane sorted by x-coordinate; this algorithm has the same complexity as the proximate points
algorithm. Next, we showed that the proximate points algorithms afford us work-time optimal (resp.
time-optimal) parallel algorithms for various fundamental digital geometry and image processing
problems. Specifically, we show that the Voronoi map, the Euclidean distance map, the maximal
empty circles, the largest empty circles, and other related problems.
Further, we have proved the work-time, respectively, the time optimality of our proximate points
and convex hull algorithms. However, for the image processing problems discussed, it is not known
whether the algorithms developed are optimal. We conjecture that, for these
problems,\Omega\Gammaobl log n)
is a time lower bound on the CRCW, provided that the algorithms are work-time optimal. For the
CREW and EREW, the logical-OR problem can be reduced to these image processing problems
quite easily.
Therefore,\Omega\Gammaher n) is a time lower bound for both the CREW and the EREW.



--R

Euclidean distance transform on polymorphic processor array.
Computer Vision.
A fast parallel algorithm for finding the convex hull of a sorted point set.

Centres of maximal discs in the 5-7-11 distance transform
Linear time Euclidean distance transform algorithms.
Efficient geometric algorithms on the EREW PRAM.
Optimal algorithm for complete Euclidean distance transform.
Designing systolic architectures for complete Euclidean distance trans- form
An efficient algorithm for complete Euclidean distance transform on mesh-connected SIMD
A fast algorithm for Euclidean distance maps of a 2-d binary image
Optimal parallel algorithms for computing convex hulls.
A parallel method for the prefix convex hulls problem.
SIMD hypercube algorithm for complete Euclidean distance transform.
Upper and lower time bounds for parallel random access machines without simultaneous writes.
On computing Voronoi diagrams for sorted point sets.
Pattern Classification and Scene Analysis
An optimal parallel algorithm for the Euclidean distance maps.
A unified linear-time algorithm for computing distance maps
An Introduction to Parallel Algorithms.
Efficient parallel geometric algorithms on a mesh of trees.
Parallel computation of exact Euclidean distance transform.
Modified distance transform with raster scanning value propagation.
Efficient algorithms for the Euclidean distance transform.
Computational Geometry: An Introduction.


A skeletonization algorithm by maxima tracking on Euclidean distance transform.
Finding the maximum
Parallelism in comparison problem.
On the generation of skeletons from discrete Euclidean distance maps.
--TR

--CTR
Ling Chen , Yi Pan , Xiao-hua Xu, Scalable and Efficient Parallel Algorithms for Euclidean Distance Transform on the LARPBS Model, IEEE Transactions on Parallel and Distributed Systems, v.15 n.11, p.975-982, November 2004
Amitava Datta , Subbiah Soundaralakshmi, Fast and scalable algorithms for the Euclidean distance transform on a linear array with a reconfigurable pipelined bus system, Journal of Parallel and Distributed Computing, v.64 n.3, p.360-369, March 2004
parallel algorithms;proximate points;digital geometry;convex hulls;pattern recognition;largest empty circles;cellular systems;image analysis