- Email: [email protected]
Computing a maxian point of a simple rectilinear polygon
Operations Research Letters
Operations Research Letters 35 (2007) 54 – 60
www.elsevier.com/locate/orl
Computing a maxian point of a simple rectilinear polygon D.P. Wang∗ Department of International Trade, National Pingtung Institute of Commerce, Pingtung, Taiwan, ROC Received 10 March 2005; accepted 19 December 2005 Available online 7 February 2006
Abstract Let P be a simple rectilinear polygon with n vertices. There are k points in P . The maxian problem is to locate a single facility in P so as to maximize the sum of its distance from it to the k points. We present an O((n × k) log n) time algorithm for this problem. © 2006 Elsevier B.V. All rights reserved. Keywords: Computational geometry; Facility location; Divide and conquer; Analysis of algorithms
1. Introduction The concept of 1-maxian (1-antimedian, 1maxisum) problem has been studied in different structures by different authors, for example, an O(n) time algorithm on a n nodes tree [13], an O(m × n) time algorithm on a n nodes and m edges network [6]. In this paper, we shall consider this problem on a simple rectilinear polygon. A simple polygon P is rectilinear if all edges of P are axis parallel. The inner angles of the vertices of P are either /2 or 3/2. The vertices of P with inner angle /2 (3/2) are called the convex (concave) vertices. Let x1 and x2 be the two points in a simple rectilinear polygon P . The inner L1 -distance between x1 and x2 , denoted as d(x1 , x2 ), is defined to be the ∗ Tel.: +886 8 7238700 6104; fax: +886 8 7238720.
E-mail address: [email protected]. 0167-6377/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.orl.2005.12.006
length of a shortest rectilinear path connecting them inside of P . From now on, for the sake of clarity, whenever we talk of a polygon, we mean a simple rectilinear polygon and whenever we talk of a distance, we mean an inner L1 -distance. Assume that X = {x1 , x2 , . . . , xk } is a k points set in a polygon P having positive weights t1 , t2 , . . . , tk , respectively. Let x be an arbitrary point of P . The total weighted distance function of x is given by F (x) =
k
ti d(x, xi ).
i=1
A point x of polygon P minimizing (maximizing) this function is called a median (maxian) point of P . Chepoi and Dragan [5] can find the median point of P in O(n + k log n) time. In this paper, we shall find the maxian point of P in O((n × k) log n) time.
D.P. Wang / Operations Research Letters 35 (2007) 54 – 60
2. Candidate points
55
we have
A line segment c is said to be a cut of P if c is an axis-parallel line segment in P that touches the boundary of P in exactly two points. From Ref. [11, p. 67], we have the following lemma. Lemma 1. Any cut of P does not touch any convex vertex of P . P
of P , we use the notation For any subpolygon T (P ) to denote the summation of the weights of the points of X in P . That is, T (P )= xi ∈P ti . Suppose e is an edge of P , we use P − e to denote the region of P but exclusive of e. Lemma 2. There exists a maxian point of P which is a convex vertex of P . Proof. Assume that x is a point of P and not a convex vertex of P . It is clear that x is on the interior of at least one cut. Without loss of generality, assume x is on a vertical cut vx . Let Pl and Pr be the two subpolygons of P defined by the vertical cut vx , where Pl is to the left of Pr . Let xl and xr be the two points of Pl and Pr , respectively, and at a horizontal distance from x, as shown in Fig. 1. We can assume that T (Pl ) T (Pr − vx ); otherwise, we can horizontal flip the polygon P to satisfy it. For any point xi ∈ X and in Pl , we have d(xr , x) = d(xr , xi ) − d(x, xi ). And for any point xi ∈ X and in Pr − vx , we have d(xr , x)d(x, xi ) − d(xr , xi ). Thus,
F (xr ) − F (x) =
k
ti (d(xr , xi ) − d(x, xi ))
i=1
⎡ =⎣
⎤ ti (d(xr , xi ) − d(x, xi ))⎦
xi ∈Pl
⎡
+⎣ ⎡ =⎣
xi ∈Pr −vx
⎤ ti (d(xr , xi ) − d(x, xi ))⎦ ⎤
ti d(xr , x)⎦
xi ∈Pl
⎡
−⎣ ⎡ ⎣
xi ∈Pr −vx
xi ∈Pl
⎤ ti (d(x, xi ) − d(xr , xi ))⎦ ⎤
⎡
ti d(xr , x)⎦ − ⎣
⎤ ti d(xr , x)⎦
xi ∈Pr −vx
= (T (Pl ) − T (Pr − vx ))d(xr , x) 0. Therefore, we can conclude that, for any point x of P , if x is not a convex vertex of P , then there always exists a point y of P such that F (y) F (x).
3. Geometric analysis and the algorithm Our divide and conquer algorithm is based on the rectilinear polygon cutting theorem [1]. The theorem can be stated as follows. Lemma 3. Let P be a polygon having n vertices. Then there exists a cut that splits P into two subpolygons Pl and Pr of almost equal size, that is, the number of edges in both Pl and Pr is less than or equal to 3 4 (n + 2).
Fig. 1.
Let vi be a vertex of P and X be a subset of X. We use F (vi , X ) to denote the summation of the weighed distances from the points of X to vi . That is, F (vi , X ) = xj ∈X tj d(vi , xj ). The idea of our algorithm is to compute F (vi ) for each vertex vi of P . And then choose the one which is a convex vertex of P and has a maximum F -value.
56
D.P. Wang / Operations Research Letters 35 (2007) 54 – 60
We compute F (vi ) for each vertex vi of P as follows: We set F (vi ) as a global variable for each vertex vi of P . If P is a rectangle, then we can compute each F (vi ) by a bruteforce method; otherwise, based upon Lemma 3, we can choose a cut c that splits P into two subpolygons Pl and Pr of almost equal size. We partition X into two subset Xl and Xr where Xl = {xi | xi ∈ X and xi is in Pl } and Xr = {xi | xi ∈ X and xi is in Pr − c}. Assume that vi is a vertex of Pl , by definition, F (vi ) = F (vi , X) = F (vi , Xl ) + F (vi , Xr ). We can design a data structure, which is called a histogram partition of a polygon P , to compute F (vi , Xr ) efficiently. And we can use recursive algorithm to compute F (vi , Xl ). If vi is a vertex of Pr , we can use the similar strategy to find F (vi ). Assume that vi is a vertex of Pl . Let us discuss how to find F (vi , Xr ). From Lemma 1 of [12], we know that for each point x of P and a cut c, there is one unique point on c that minimizes the distance from x to c. We call the unique point on c that is closest to x the projection point of x on c and denote it by c (x). The projection point has the following property. Lemma 4 (Schuierer [12]). Let Pl and Pr be the two subpolygons of P defined by a cut c. Assume that x is in Pl and y in Pr , then d(x, y) = d(x, c (x)) + d(c (x), c (y)) + d(c (y), y). Assume that vi is a vertex of Pl and xj is a point of Xr . Based upon Lemma 4, we have d(vi , xj ) = d(vi , c (vi )) + d(c (vi ), c (xj )) + d(c (xj ), xj ). Schuierer [12] has the following lemma. Lemma 5. Let P be a polygon with n vertices and e an edge of P , we can preprocess P in linear time such that, for any point x in P , the projection point of x on e, e (x), can be found in O(log n) time and the distance between x and e (x), i.e., d(x, e (x)), can also be found in O(log n) time. If x is a vertex of P , then both of e (x) and d(x, e (x)) can be found in constant time. Let vi be a vertex of Pl and xj be a point of Xr . From Lemmas 4 and 5, we can find d(vi , xj ) in O(log n) time if P has been preprocessed in O(n) time. We shall use a data structure to improve its time complexity to constant time while keeping the preprocessing time linear.
Fig. 2. The histogram partition of P w.r.t. an edge e.
The structure upon which we can find d(vi , xj ) in O(1) time is the histogram partition of a polygon P . It is the decomposition of P into simpler parts, called histograms. A histogram is a rectilinear polygon H that has one distinguished edge, called the base of H , whose length is equal to the sum of the lengths of all other edges parallel to this base. We say H is a maximal histogram of e if e is the base of H and there is no other histogram in P with the same base that properly contains H . A window of the maximal histogram H is a maximal line segment on the boundary of H that is not part of the boundary of P . The histogram partition of P w.r.t. an edge e, denoted by H (P , e), can be defined as follows: Let He be the maximal histogram of e in P . If P = He , then H (P , e) is defined to be {He }. Otherwise, let wi , 1 i k, be the windows of He . Each window wi splits P into two subpolygons one of which does not contain He . We define this subpolygon as Pwi . The histogram partition of Pw.r.t. an edge e is given recursively by H (P , e) = 1 i k H (Pwi , wi ) ∪ {He }. Please see Fig. 2 for an illustration. There is a natural directed tree structure associated to H (P , e). We call it the histogram tree of H (P , e) and denote it T (P , e). The nodes of T (P , e) are the histograms of H (P , e). He is the root of T (P , e). There is a directed edge from Ha to Hb in T (P , e) if the base of Hb is a window of Ha . We shall use the terminology of trees to the histogram partition. For example, for any two histograms Ha and Hb of H (P , e), we say Ha is the father of Hb if Ha is the father of Hb in T (P , e), and Ha and Hb are brothers if they have a common father in T (P , e). Levcopoulos [8] can find H (P , e) and T (P , e) in O(n) time.
D.P. Wang / Operations Research Letters 35 (2007) 54 – 60
57
exists a shortest path from x to e (x) such that it passes through v. Therefore, e (x) = e (v).
Fig. 3.
Let Hwi be a histogram of the histogram partition H (P , e). Assume that point x is in Hwi and wi = e. We call the end point of wi that is closer to the edge e the gate point of x. Lemma 6. Let P be a polygon and e be an edge of P . Let x be a point in H (P , e). If v is the gate point of x, then v is a vertex of P and e (x) = e (v). Proof. Let Hwi and Hwj be the two histograms in H (P , e) where Hwj is the father of Hwi . Suppose that x is in Hwi , as shown in Fig. 3. (1) v is a vertex of P : Please refer to Lemma 2.1 of [3]. (2) e (x)=e (v): We prove this by showing that there exists a shortest path from x to e (x) that passes through v. Let sp(x, e (x)) be a shortest path from x to e (x) that does not pass through v, as shown in Fig. 3. It is obvious that sp(x, e (x)) must intersect wi and wj . Let a and b denote the first intersection points of sp(x, e (x)) with wi and wj , respectively. From Lemma 5 of [10], the subpath of sp(x, e (x)) from a to b is a shortest path from a to b. From Lemma 1 of [9], the shortest path from a to b must be fully inside of Hwj . A rectilinear path is called a staircase if any axis-parallel line does not intersect the path or intersects the path in exactly one point or a line segment. From Lemma 1 of [4], any shortest path from a to b is a staircase, as shown in Fig. 3. We can transform this staircase by sliding such that it passes through v without increasing its length. This implies that there
Assume that vj is the gate point of x. From Lemma 6, we have d(x, e (x)) = d(x, e (vj )) = d(x, vj ) + d(vj , e (vj )). Based upon Lemma 5, since vj is a vertex of P , d(vj , e (vj )) can be found in O(1) time after P has been preprocessed in O(n) time. Now, what remains to be explained is how to find the gate point of x in constant time. In order to do this, we must determine which histogram of H (P , e) contains x in O(1) time. We can solve it by using a horizontal segment hx which passes through x. Let x ∈ X, we define hx as follows: If x is on a horizontal edge vj vj +1 of P , then hx = vj vj +1 ; otherwise, hx is defined to be the horizontal cut of P that passes through x. For the sake of clarity, assume that hx = xl xr where xl is the left endpoint of hx and xr is the right of it. Goodrich and Tamassia [7] can preprocess P in O(n) time, such that hx can be found in O(log n) time. We use Hx , Hxl and Hxr to denote the histograms of H (P , e) that contain x, xl and xr , respectively. There is a relationship between these histograms. First of all, From Lemma 2.1 of [2], hx can intersect at most three histograms of H (P , e). And we have the following lemma about the relationship between Hxl and Hxr . Lemma 7. The relationship between Hxl and Hxr must be falling in one of the following three cases: (1) Hxl = Hxr . (2) Hxl is the father of Hxr or Hxr is the father of Hxl . (3) Hxl and Hxr are brothers. Proof. There are three possibilities: Case 1: hx intersects one histogram of H (P , e). In this case, hx is fully inside of one histogram, so Hxl = Hxr . Case 2: hx intersects two histograms of H (P , e). In this case, the two histograms must have the father–son relationship. Therefore, Hxl is the father of Hxr or Hxr is the father of Hxl . Case 3: hx intersects three histograms of H (P , e). In this case, the three histograms must have the father–2sons relationship. Therefore, Hxl and Hxr are brothers.
58
D.P. Wang / Operations Research Letters 35 (2007) 54 – 60
Lemma 8. Suppose that we have found Hxl and Hxr , then we can find Hx in O(1) time. Proof. Depending upon the cases of Lemma 7, we have the following three cases: (1) Hxl = Hxr . In this case, Hx = Hxl = Hxr . (2) Hxl is the father of Hxr or Hxr is the father of Hxl . Let w be the window intersected by hx . Note that w is a vertical cut and also a base of Hxl or Hxr . We can find Hx as follows: if xx l does not intersect w, then Hx = Hxl , else Hx = Hxr . (3) Hxl and Hxr are brothers. Let wl and wr be the windows intersected by hx . Note that wl and wr are vertical cuts and are bases of Hxl and Hxr , respectively. We can find Hx as follows: If xx l does not intersect wl , Hx = Hxl ; Otherwise if it intersects wr , Hx = Hxr ; Otherwise Hx is their common father. It takes O(1) time to do the above three cases. The remaining question is: How to find Hxl and Hxr in O(1) time. Let us discuss how to find Hxl . We can find Hxr in a similar way. When xl is a vertex of P , Schuierer [12] can find Hxl in O(1) time. Therefore, we can assume that xl is not a vertex of P and on the edge vi vi+1 . Since hx is horizontal, vi vi+1 is a vertical edge of P . Let Hvi and Hvi+1 be the histograms which contain vi and vi+1 , respectively. We have the following lemma about the relationship between Hvi and Hvi+1 . Lemma 9. The relationship between Hvi and Hvi+1 must be falling in one of the following three cases: 1. Hvi = Hvi+1 . 2. Hvi is the father of Hvi+1 or Hvi+1 is the father of Hvi . 3. Hvi and Hvi+1 are brothers. Proof. vi vi+1 is a vertical edge of P . We can rotate P to change vi vi+1 as a horizontal edge of P . It is a special case of Lemma 7. Using the proof similar to that of Lemma 8, we can have the following lemma. Lemma 10. Suppose that we have found Hvi and Hvi+1 , then we can find Hxl in O(1) time.
Since vi and vi+1 are vertices of P , Schuierer [12] can find Hvi and Hvi+1 in O(1) time. And from Lemmas 8 and 10, we can find Hx in O(1) time. Based upon this result, we can have the following lemma. Lemma 11. Let Pl and Pr be the two subpolygons of P defined by a cut c. Assume that vi is a vertex of Pl and x ∈ Xr . After preprocessing P in O(n) time, we can find d(vi , x) in O(1) time. Proof. From Lemma 4, d(vi , x) = d(vi , c (vi )) + d(c (vi ), c (x)) + d(c (x), x). Assume that v is the gate point of x in H (Pr , c). From Lemma 6, we have d(vi , x) = d(vi , c (vi )) + d(c (vi ), c (v)) + d(c (v), x). Since d(c (v), x) = d(c (v), v) + d(v, x), d(vi , x) = d(vi , c (vi )) + d(c (vi ), c (v)) + d(c (v), v) + d(v, x). From Lemmas 8 and 10, we can find Hx in O(1) time, therefore, v can also be found in O(1) time. From Lemma 5, vi and v are vertices of P , the first three items can be found in O(1) time. By definition, the last item can also be computed in O(1) time. Thus, we can find d(vi , x) in O(1) time. Assume that x ∈ X and hx = xl xr where xl is the left-endpoint of hx and xr is the right of it. Since hx is defined relative to P , when P is split into subpolygons, the definition of hx should be updated relative to the part x lies in. This may be done as follows: (1) If hx does not intersect c or hx = c, then it is not necessary to change hx . (2) If hx intersects c at y, then we have the following two subcases: (2.1) If x is in Pl , then change hx into xl y. (2.2) If x is in Pr − c, then change hx into yx r . Based upon the above, we can have the following algorithms to solve our problem. Algorithm A Input: A polygon P ={v1 →v2 →· · ·→vn } of n vertices and a point set X = {x1 , x2 , . . ., xk } each has a positive weight ti . Output: The maxian point of P with respect to a point set X.
D.P. Wang / Operations Research Letters 35 (2007) 54 – 60
Step 1. Set F (vi ) = 0 for all vertex vi of P . /*F (vi ) is a global variable*/ Step 2. For each point x of X, find hx . Let hX = {hx | x ∈ X}. Step 3. Call Algorithm B(P , X, hX ). Step 4. Return vi to be a maxian point of P if vi is a convex vertex of P with a maximum F -value; i.e., F (vi ) F (vj ), for any other convex vertex vj of P . Algorithm B(P, X, hX ) Input: A polygon P , a points set X and hX which is a set of horizontal segments hx (x ∈ X). Output: Compute F (vi ) for each vertex vi of P . Step 1. If P is a rectangle, then for each vertex vi of P , compute F (vi , X) = t xj ∈X j d(vi , xj ), set F (vi ) = F (vi ) + F (vi , X) and return. Step 2. Compute a cut c which splits P into two subpolygons Pl and Pr such that the number of edges in both Pl and Pr is less than or equal to (3/4)n + 2. Step 3. Construct H (Pl , c) and H (Pr , c). Step 4. Partition X into Xl and Xr where Xl =X ∩ Pl and Xr = X ∩ (Pr − c). Step 5. For each point x of Xl (Xr , resp.), assume hx = xl xr . If hx intersects c at a point y, then change hx into xl y (yx r , resp.). Let hXl ={hx | x ∈ Xl } and hXr ={hx | x∈Xr }. Step 6. For each vertex vi of Pl , compute F (vi , Xr ) = xj ∈Xr tj d(vi , xj ) and set F (vi ) = F (vi ) + F (vi , Xr ). Step 7. For each vertex vi of Pr , compute F (vi , Xl ) = xj ∈Xl tj d(vi , xj ) and set F (vi ) = F (vi ) + F (vi , Xl ). Step 8. Call B(Pl , Xl , hXl ). Step 9. Call B(Pr , Xr , hXr ). Based on Lemma 2, our maxian point is a convex vertex of P . And from Lemma 1, any cut of P does not touch any convex vertex of P . Therefore, our candidate points do not disappear while cutting the polygon. For each vertex vi of P , in Algorithm B, we sum the weighted distances from vi to all points of X. Therefore, Algorithm B correctly computes F (vi ) for each vertex vi of P .
59
Let us analyze the time complexity of Algorithm B. Step 1 takes O(k) time since there are at most k points of X in a rectangle. Step 2 can be done in O(n) time [1]. Step 3 can be done in O(n) time [8]. For each point x of X, the work we have to do in Step 4 is to determine whether x is in Pl or in Pr − c. This can be done by scanning the boundary of P to test whether both of the two end points of hx are in Pl , in Pr or one is in Pl and the other in Pr . Let hx = xl xr , we have the following three cases to test whether x is in Pl or in Pr − c: (1) xl and xr are in Pl : x is in Pl . (2) xl and xr are in Pr − c: x is in Pr − c. (3) xl is in Pl and xr is in Pr − c: If the line segment xx r intersects c, then x is in Pl , else x is in Pr − c. Therefore, it takes O(n + k) time to do Step 4. It needs O(k) time to execute Step 5 since it is only a line segments intersection test. For each vertex vi of Pl and any point xj in Xr , based upon Lemma 11, we can compute d(vi , xj ) in O(1) time. There are at most k points of Xr , therefore F (vi , Xr ) can be computed in O(k) time. Since there are at most n vertices of Pl , Step 6 can be executed in O(n × k) time. The time complexity of Step 7 is identical with that of Step 6, therefore, it also takes O(n × k) time. From the above discussion, we conclude that the divide, conquer and merge time of Algorithm B is O(n × k) and the termination time of it is O(k). Let T (n) denote the time complexity of Algorithm B, we have the following recurrence relation for T (n): T (n)2 × T ( 43 n) + O(n × k) for n > 4, T (n) = O(k) for n = 4.
(1)
From the recurrence relation, we have T (n) = O((n × k) log n). It is easy to see that Step 3 is the bottleneck step of Algorithm A. Therefore, the time complexity of Algorithm A is O((n × k) log n). Combining the above discussion, we have the main result of this paper. Theorem 1. It needs O((n×k) log n) time to compute a maxian point of a polygon P with a point set X, where n is the number of vertices of P and k is the number of points of X.
60
D.P. Wang / Operations Research Letters 35 (2007) 54 – 60
4. Conclusion Our algorithm is simple. We think this algorithm may have some applications in the VLSI layout design and the obnoxious facility location problems. A natural generalization of the considered problem is to improve the time complexity of it. We can do this by improving the time complexity of Steps 6 and 7 of Algorithm B. It is always possible if we can find some properties about the vertices of P . For example, assume vi and vj are two vertices of Pl . If F (vi , Xr ) − F (vj , Xr ) can be found in constant time. Then, the time complexity of Step 6 can be improved to O(n + k). Another research direction is to modify the objective function, for example, using rectilinear link distance metric [1]. That is, find a point x in P , such that the summation of the weighted link distances from x to all points of X is maximized. Acknowledgments The authors wish to thank Prof. B.C. Liaw for his helpful discussion about this problem. We also thank the anonymous referees for the careful reading of the paper and give our constructive comments and suggestions which not only greatly improve our presentation of the results but also lead to the correction of some mistakes. References [1] M. de Berg, On rectilinear link distance, Comput. Geom.: Theory Appl. 1 (1991) 13–34.
[2] M. de Berg, M. Van Kreveld, Rectilinear decompositions with low stabbing number, Inform. Process. Lett. 52 (1994) 215–221. [3] M. de Berg, M. Van Kreveld, B.J. Nilsson, M.H. Overmars, Finding shortest paths in the presence of orthogonal obstacles using a combined L1 and link metric, Lecture Notes Comput. Sci. 447 (1990) 211–224. [4] D.Z. Chen, K.S. Klenk, H.T. Tu, Shortest path queries among weighted obstacles in the rectilinear plane, SIAM J. Comput. 29 (2000) 1223–1246. [5] V. Chepoi, F. Dragan, Computing a median point of a simple rectilinear polygon, Inform. Process. Lett. 49 (1994) 281–285. [6] M. Colebrook, J. Guttierez, J. Sicilia, A new bound and an O(mn) algorithm for the undesirable 1-median problem (maxian) on networks, Comput. Oper. Res. 32 (2005) 309–325. [7] M.T. Goodrich, R. Tamassia, Dynamic ray shooting and shortest path via balanced geodesic triangulation, in: Proceedings of the Ninth Annual ACM Symposium on Computational Geometry, 1993, pp. 318–327. [8] C. Levcopoulos, Heuristics for minimum decompositions of polygons, Ph.D. Thesis, University of Linköping, Linköping, Sweden, 1987. [9] A. Maheshwari, J. Sack, Simple optimal algorithms for rectilinear link path and polygon separation problem, Parallel Process. Lett. 9 (1999) 31–42. [10] K.M. McDonald, J.G. Peters, Smallest paths in simple rectilinear polygons, IEEE Trans. Comput.-Aided Des. 11 (1992) 864–875. [11] J. O’Rourke, Art Gallery Theorems and Algorithm, Oxford University Press, Oxford, 1987. [12] S. Schuierer, An optimal data structure for shortest rectilinear path queries in a simple rectilinear polygon, Int. J. Computat. Geom. Appl. 6 (1996) 205–225. [13] S.S. Ting, A linear time algorithm for maxisum facility location on tree networks, Transport. Sci. 18 (1984) 76–84.
Operations Research Letters 35 (2007) 54 – 60
www.elsevier.com/locate/orl
Computing a maxian point of a simple rectilinear polygon D.P. Wang∗ Department of International Trade, National Pingtung Institute of Commerce, Pingtung, Taiwan, ROC Received 10 March 2005; accepted 19 December 2005 Available online 7 February 2006
Abstract Let P be a simple rectilinear polygon with n vertices. There are k points in P . The maxian problem is to locate a single facility in P so as to maximize the sum of its distance from it to the k points. We present an O((n × k) log n) time algorithm for this problem. © 2006 Elsevier B.V. All rights reserved. Keywords: Computational geometry; Facility location; Divide and conquer; Analysis of algorithms
1. Introduction The concept of 1-maxian (1-antimedian, 1maxisum) problem has been studied in different structures by different authors, for example, an O(n) time algorithm on a n nodes tree [13], an O(m × n) time algorithm on a n nodes and m edges network [6]. In this paper, we shall consider this problem on a simple rectilinear polygon. A simple polygon P is rectilinear if all edges of P are axis parallel. The inner angles of the vertices of P are either /2 or 3/2. The vertices of P with inner angle /2 (3/2) are called the convex (concave) vertices. Let x1 and x2 be the two points in a simple rectilinear polygon P . The inner L1 -distance between x1 and x2 , denoted as d(x1 , x2 ), is defined to be the ∗ Tel.: +886 8 7238700 6104; fax: +886 8 7238720.
E-mail address: [email protected]. 0167-6377/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.orl.2005.12.006
length of a shortest rectilinear path connecting them inside of P . From now on, for the sake of clarity, whenever we talk of a polygon, we mean a simple rectilinear polygon and whenever we talk of a distance, we mean an inner L1 -distance. Assume that X = {x1 , x2 , . . . , xk } is a k points set in a polygon P having positive weights t1 , t2 , . . . , tk , respectively. Let x be an arbitrary point of P . The total weighted distance function of x is given by F (x) =
k
ti d(x, xi ).
i=1
A point x of polygon P minimizing (maximizing) this function is called a median (maxian) point of P . Chepoi and Dragan [5] can find the median point of P in O(n + k log n) time. In this paper, we shall find the maxian point of P in O((n × k) log n) time.
D.P. Wang / Operations Research Letters 35 (2007) 54 – 60
2. Candidate points
55
we have
A line segment c is said to be a cut of P if c is an axis-parallel line segment in P that touches the boundary of P in exactly two points. From Ref. [11, p. 67], we have the following lemma. Lemma 1. Any cut of P does not touch any convex vertex of P . P
of P , we use the notation For any subpolygon T (P ) to denote the summation of the weights of the points of X in P . That is, T (P )= xi ∈P ti . Suppose e is an edge of P , we use P − e to denote the region of P but exclusive of e. Lemma 2. There exists a maxian point of P which is a convex vertex of P . Proof. Assume that x is a point of P and not a convex vertex of P . It is clear that x is on the interior of at least one cut. Without loss of generality, assume x is on a vertical cut vx . Let Pl and Pr be the two subpolygons of P defined by the vertical cut vx , where Pl is to the left of Pr . Let xl and xr be the two points of Pl and Pr , respectively, and at a horizontal distance from x, as shown in Fig. 1. We can assume that T (Pl ) T (Pr − vx ); otherwise, we can horizontal flip the polygon P to satisfy it. For any point xi ∈ X and in Pl , we have d(xr , x) = d(xr , xi ) − d(x, xi ). And for any point xi ∈ X and in Pr − vx , we have d(xr , x)d(x, xi ) − d(xr , xi ). Thus,
F (xr ) − F (x) =
k
ti (d(xr , xi ) − d(x, xi ))
i=1
⎡ =⎣
⎤ ti (d(xr , xi ) − d(x, xi ))⎦
xi ∈Pl
⎡
+⎣ ⎡ =⎣
xi ∈Pr −vx
⎤ ti (d(xr , xi ) − d(x, xi ))⎦ ⎤
ti d(xr , x)⎦
xi ∈Pl
⎡
−⎣ ⎡ ⎣
xi ∈Pr −vx
xi ∈Pl
⎤ ti (d(x, xi ) − d(xr , xi ))⎦ ⎤
⎡
ti d(xr , x)⎦ − ⎣
⎤ ti d(xr , x)⎦
xi ∈Pr −vx
= (T (Pl ) − T (Pr − vx ))d(xr , x) 0. Therefore, we can conclude that, for any point x of P , if x is not a convex vertex of P , then there always exists a point y of P such that F (y) F (x).
3. Geometric analysis and the algorithm Our divide and conquer algorithm is based on the rectilinear polygon cutting theorem [1]. The theorem can be stated as follows. Lemma 3. Let P be a polygon having n vertices. Then there exists a cut that splits P into two subpolygons Pl and Pr of almost equal size, that is, the number of edges in both Pl and Pr is less than or equal to 3 4 (n + 2).
Fig. 1.
Let vi be a vertex of P and X be a subset of X. We use F (vi , X ) to denote the summation of the weighed distances from the points of X to vi . That is, F (vi , X ) = xj ∈X tj d(vi , xj ). The idea of our algorithm is to compute F (vi ) for each vertex vi of P . And then choose the one which is a convex vertex of P and has a maximum F -value.
56
D.P. Wang / Operations Research Letters 35 (2007) 54 – 60
We compute F (vi ) for each vertex vi of P as follows: We set F (vi ) as a global variable for each vertex vi of P . If P is a rectangle, then we can compute each F (vi ) by a bruteforce method; otherwise, based upon Lemma 3, we can choose a cut c that splits P into two subpolygons Pl and Pr of almost equal size. We partition X into two subset Xl and Xr where Xl = {xi | xi ∈ X and xi is in Pl } and Xr = {xi | xi ∈ X and xi is in Pr − c}. Assume that vi is a vertex of Pl , by definition, F (vi ) = F (vi , X) = F (vi , Xl ) + F (vi , Xr ). We can design a data structure, which is called a histogram partition of a polygon P , to compute F (vi , Xr ) efficiently. And we can use recursive algorithm to compute F (vi , Xl ). If vi is a vertex of Pr , we can use the similar strategy to find F (vi ). Assume that vi is a vertex of Pl . Let us discuss how to find F (vi , Xr ). From Lemma 1 of [12], we know that for each point x of P and a cut c, there is one unique point on c that minimizes the distance from x to c. We call the unique point on c that is closest to x the projection point of x on c and denote it by c (x). The projection point has the following property. Lemma 4 (Schuierer [12]). Let Pl and Pr be the two subpolygons of P defined by a cut c. Assume that x is in Pl and y in Pr , then d(x, y) = d(x, c (x)) + d(c (x), c (y)) + d(c (y), y). Assume that vi is a vertex of Pl and xj is a point of Xr . Based upon Lemma 4, we have d(vi , xj ) = d(vi , c (vi )) + d(c (vi ), c (xj )) + d(c (xj ), xj ). Schuierer [12] has the following lemma. Lemma 5. Let P be a polygon with n vertices and e an edge of P , we can preprocess P in linear time such that, for any point x in P , the projection point of x on e, e (x), can be found in O(log n) time and the distance between x and e (x), i.e., d(x, e (x)), can also be found in O(log n) time. If x is a vertex of P , then both of e (x) and d(x, e (x)) can be found in constant time. Let vi be a vertex of Pl and xj be a point of Xr . From Lemmas 4 and 5, we can find d(vi , xj ) in O(log n) time if P has been preprocessed in O(n) time. We shall use a data structure to improve its time complexity to constant time while keeping the preprocessing time linear.
Fig. 2. The histogram partition of P w.r.t. an edge e.
The structure upon which we can find d(vi , xj ) in O(1) time is the histogram partition of a polygon P . It is the decomposition of P into simpler parts, called histograms. A histogram is a rectilinear polygon H that has one distinguished edge, called the base of H , whose length is equal to the sum of the lengths of all other edges parallel to this base. We say H is a maximal histogram of e if e is the base of H and there is no other histogram in P with the same base that properly contains H . A window of the maximal histogram H is a maximal line segment on the boundary of H that is not part of the boundary of P . The histogram partition of P w.r.t. an edge e, denoted by H (P , e), can be defined as follows: Let He be the maximal histogram of e in P . If P = He , then H (P , e) is defined to be {He }. Otherwise, let wi , 1 i k, be the windows of He . Each window wi splits P into two subpolygons one of which does not contain He . We define this subpolygon as Pwi . The histogram partition of Pw.r.t. an edge e is given recursively by H (P , e) = 1 i k H (Pwi , wi ) ∪ {He }. Please see Fig. 2 for an illustration. There is a natural directed tree structure associated to H (P , e). We call it the histogram tree of H (P , e) and denote it T (P , e). The nodes of T (P , e) are the histograms of H (P , e). He is the root of T (P , e). There is a directed edge from Ha to Hb in T (P , e) if the base of Hb is a window of Ha . We shall use the terminology of trees to the histogram partition. For example, for any two histograms Ha and Hb of H (P , e), we say Ha is the father of Hb if Ha is the father of Hb in T (P , e), and Ha and Hb are brothers if they have a common father in T (P , e). Levcopoulos [8] can find H (P , e) and T (P , e) in O(n) time.
D.P. Wang / Operations Research Letters 35 (2007) 54 – 60
57
exists a shortest path from x to e (x) such that it passes through v. Therefore, e (x) = e (v).
Fig. 3.
Let Hwi be a histogram of the histogram partition H (P , e). Assume that point x is in Hwi and wi = e. We call the end point of wi that is closer to the edge e the gate point of x. Lemma 6. Let P be a polygon and e be an edge of P . Let x be a point in H (P , e). If v is the gate point of x, then v is a vertex of P and e (x) = e (v). Proof. Let Hwi and Hwj be the two histograms in H (P , e) where Hwj is the father of Hwi . Suppose that x is in Hwi , as shown in Fig. 3. (1) v is a vertex of P : Please refer to Lemma 2.1 of [3]. (2) e (x)=e (v): We prove this by showing that there exists a shortest path from x to e (x) that passes through v. Let sp(x, e (x)) be a shortest path from x to e (x) that does not pass through v, as shown in Fig. 3. It is obvious that sp(x, e (x)) must intersect wi and wj . Let a and b denote the first intersection points of sp(x, e (x)) with wi and wj , respectively. From Lemma 5 of [10], the subpath of sp(x, e (x)) from a to b is a shortest path from a to b. From Lemma 1 of [9], the shortest path from a to b must be fully inside of Hwj . A rectilinear path is called a staircase if any axis-parallel line does not intersect the path or intersects the path in exactly one point or a line segment. From Lemma 1 of [4], any shortest path from a to b is a staircase, as shown in Fig. 3. We can transform this staircase by sliding such that it passes through v without increasing its length. This implies that there
Assume that vj is the gate point of x. From Lemma 6, we have d(x, e (x)) = d(x, e (vj )) = d(x, vj ) + d(vj , e (vj )). Based upon Lemma 5, since vj is a vertex of P , d(vj , e (vj )) can be found in O(1) time after P has been preprocessed in O(n) time. Now, what remains to be explained is how to find the gate point of x in constant time. In order to do this, we must determine which histogram of H (P , e) contains x in O(1) time. We can solve it by using a horizontal segment hx which passes through x. Let x ∈ X, we define hx as follows: If x is on a horizontal edge vj vj +1 of P , then hx = vj vj +1 ; otherwise, hx is defined to be the horizontal cut of P that passes through x. For the sake of clarity, assume that hx = xl xr where xl is the left endpoint of hx and xr is the right of it. Goodrich and Tamassia [7] can preprocess P in O(n) time, such that hx can be found in O(log n) time. We use Hx , Hxl and Hxr to denote the histograms of H (P , e) that contain x, xl and xr , respectively. There is a relationship between these histograms. First of all, From Lemma 2.1 of [2], hx can intersect at most three histograms of H (P , e). And we have the following lemma about the relationship between Hxl and Hxr . Lemma 7. The relationship between Hxl and Hxr must be falling in one of the following three cases: (1) Hxl = Hxr . (2) Hxl is the father of Hxr or Hxr is the father of Hxl . (3) Hxl and Hxr are brothers. Proof. There are three possibilities: Case 1: hx intersects one histogram of H (P , e). In this case, hx is fully inside of one histogram, so Hxl = Hxr . Case 2: hx intersects two histograms of H (P , e). In this case, the two histograms must have the father–son relationship. Therefore, Hxl is the father of Hxr or Hxr is the father of Hxl . Case 3: hx intersects three histograms of H (P , e). In this case, the three histograms must have the father–2sons relationship. Therefore, Hxl and Hxr are brothers.
58
D.P. Wang / Operations Research Letters 35 (2007) 54 – 60
Lemma 8. Suppose that we have found Hxl and Hxr , then we can find Hx in O(1) time. Proof. Depending upon the cases of Lemma 7, we have the following three cases: (1) Hxl = Hxr . In this case, Hx = Hxl = Hxr . (2) Hxl is the father of Hxr or Hxr is the father of Hxl . Let w be the window intersected by hx . Note that w is a vertical cut and also a base of Hxl or Hxr . We can find Hx as follows: if xx l does not intersect w, then Hx = Hxl , else Hx = Hxr . (3) Hxl and Hxr are brothers. Let wl and wr be the windows intersected by hx . Note that wl and wr are vertical cuts and are bases of Hxl and Hxr , respectively. We can find Hx as follows: If xx l does not intersect wl , Hx = Hxl ; Otherwise if it intersects wr , Hx = Hxr ; Otherwise Hx is their common father. It takes O(1) time to do the above three cases. The remaining question is: How to find Hxl and Hxr in O(1) time. Let us discuss how to find Hxl . We can find Hxr in a similar way. When xl is a vertex of P , Schuierer [12] can find Hxl in O(1) time. Therefore, we can assume that xl is not a vertex of P and on the edge vi vi+1 . Since hx is horizontal, vi vi+1 is a vertical edge of P . Let Hvi and Hvi+1 be the histograms which contain vi and vi+1 , respectively. We have the following lemma about the relationship between Hvi and Hvi+1 . Lemma 9. The relationship between Hvi and Hvi+1 must be falling in one of the following three cases: 1. Hvi = Hvi+1 . 2. Hvi is the father of Hvi+1 or Hvi+1 is the father of Hvi . 3. Hvi and Hvi+1 are brothers. Proof. vi vi+1 is a vertical edge of P . We can rotate P to change vi vi+1 as a horizontal edge of P . It is a special case of Lemma 7. Using the proof similar to that of Lemma 8, we can have the following lemma. Lemma 10. Suppose that we have found Hvi and Hvi+1 , then we can find Hxl in O(1) time.
Since vi and vi+1 are vertices of P , Schuierer [12] can find Hvi and Hvi+1 in O(1) time. And from Lemmas 8 and 10, we can find Hx in O(1) time. Based upon this result, we can have the following lemma. Lemma 11. Let Pl and Pr be the two subpolygons of P defined by a cut c. Assume that vi is a vertex of Pl and x ∈ Xr . After preprocessing P in O(n) time, we can find d(vi , x) in O(1) time. Proof. From Lemma 4, d(vi , x) = d(vi , c (vi )) + d(c (vi ), c (x)) + d(c (x), x). Assume that v is the gate point of x in H (Pr , c). From Lemma 6, we have d(vi , x) = d(vi , c (vi )) + d(c (vi ), c (v)) + d(c (v), x). Since d(c (v), x) = d(c (v), v) + d(v, x), d(vi , x) = d(vi , c (vi )) + d(c (vi ), c (v)) + d(c (v), v) + d(v, x). From Lemmas 8 and 10, we can find Hx in O(1) time, therefore, v can also be found in O(1) time. From Lemma 5, vi and v are vertices of P , the first three items can be found in O(1) time. By definition, the last item can also be computed in O(1) time. Thus, we can find d(vi , x) in O(1) time. Assume that x ∈ X and hx = xl xr where xl is the left-endpoint of hx and xr is the right of it. Since hx is defined relative to P , when P is split into subpolygons, the definition of hx should be updated relative to the part x lies in. This may be done as follows: (1) If hx does not intersect c or hx = c, then it is not necessary to change hx . (2) If hx intersects c at y, then we have the following two subcases: (2.1) If x is in Pl , then change hx into xl y. (2.2) If x is in Pr − c, then change hx into yx r . Based upon the above, we can have the following algorithms to solve our problem. Algorithm A Input: A polygon P ={v1 →v2 →· · ·→vn } of n vertices and a point set X = {x1 , x2 , . . ., xk } each has a positive weight ti . Output: The maxian point of P with respect to a point set X.
D.P. Wang / Operations Research Letters 35 (2007) 54 – 60
Step 1. Set F (vi ) = 0 for all vertex vi of P . /*F (vi ) is a global variable*/ Step 2. For each point x of X, find hx . Let hX = {hx | x ∈ X}. Step 3. Call Algorithm B(P , X, hX ). Step 4. Return vi to be a maxian point of P if vi is a convex vertex of P with a maximum F -value; i.e., F (vi ) F (vj ), for any other convex vertex vj of P . Algorithm B(P, X, hX ) Input: A polygon P , a points set X and hX which is a set of horizontal segments hx (x ∈ X). Output: Compute F (vi ) for each vertex vi of P . Step 1. If P is a rectangle, then for each vertex vi of P , compute F (vi , X) = t xj ∈X j d(vi , xj ), set F (vi ) = F (vi ) + F (vi , X) and return. Step 2. Compute a cut c which splits P into two subpolygons Pl and Pr such that the number of edges in both Pl and Pr is less than or equal to (3/4)n + 2. Step 3. Construct H (Pl , c) and H (Pr , c). Step 4. Partition X into Xl and Xr where Xl =X ∩ Pl and Xr = X ∩ (Pr − c). Step 5. For each point x of Xl (Xr , resp.), assume hx = xl xr . If hx intersects c at a point y, then change hx into xl y (yx r , resp.). Let hXl ={hx | x ∈ Xl } and hXr ={hx | x∈Xr }. Step 6. For each vertex vi of Pl , compute F (vi , Xr ) = xj ∈Xr tj d(vi , xj ) and set F (vi ) = F (vi ) + F (vi , Xr ). Step 7. For each vertex vi of Pr , compute F (vi , Xl ) = xj ∈Xl tj d(vi , xj ) and set F (vi ) = F (vi ) + F (vi , Xl ). Step 8. Call B(Pl , Xl , hXl ). Step 9. Call B(Pr , Xr , hXr ). Based on Lemma 2, our maxian point is a convex vertex of P . And from Lemma 1, any cut of P does not touch any convex vertex of P . Therefore, our candidate points do not disappear while cutting the polygon. For each vertex vi of P , in Algorithm B, we sum the weighted distances from vi to all points of X. Therefore, Algorithm B correctly computes F (vi ) for each vertex vi of P .
59
Let us analyze the time complexity of Algorithm B. Step 1 takes O(k) time since there are at most k points of X in a rectangle. Step 2 can be done in O(n) time [1]. Step 3 can be done in O(n) time [8]. For each point x of X, the work we have to do in Step 4 is to determine whether x is in Pl or in Pr − c. This can be done by scanning the boundary of P to test whether both of the two end points of hx are in Pl , in Pr or one is in Pl and the other in Pr . Let hx = xl xr , we have the following three cases to test whether x is in Pl or in Pr − c: (1) xl and xr are in Pl : x is in Pl . (2) xl and xr are in Pr − c: x is in Pr − c. (3) xl is in Pl and xr is in Pr − c: If the line segment xx r intersects c, then x is in Pl , else x is in Pr − c. Therefore, it takes O(n + k) time to do Step 4. It needs O(k) time to execute Step 5 since it is only a line segments intersection test. For each vertex vi of Pl and any point xj in Xr , based upon Lemma 11, we can compute d(vi , xj ) in O(1) time. There are at most k points of Xr , therefore F (vi , Xr ) can be computed in O(k) time. Since there are at most n vertices of Pl , Step 6 can be executed in O(n × k) time. The time complexity of Step 7 is identical with that of Step 6, therefore, it also takes O(n × k) time. From the above discussion, we conclude that the divide, conquer and merge time of Algorithm B is O(n × k) and the termination time of it is O(k). Let T (n) denote the time complexity of Algorithm B, we have the following recurrence relation for T (n): T (n)2 × T ( 43 n) + O(n × k) for n > 4, T (n) = O(k) for n = 4.
(1)
From the recurrence relation, we have T (n) = O((n × k) log n). It is easy to see that Step 3 is the bottleneck step of Algorithm A. Therefore, the time complexity of Algorithm A is O((n × k) log n). Combining the above discussion, we have the main result of this paper. Theorem 1. It needs O((n×k) log n) time to compute a maxian point of a polygon P with a point set X, where n is the number of vertices of P and k is the number of points of X.
60
D.P. Wang / Operations Research Letters 35 (2007) 54 – 60
4. Conclusion Our algorithm is simple. We think this algorithm may have some applications in the VLSI layout design and the obnoxious facility location problems. A natural generalization of the considered problem is to improve the time complexity of it. We can do this by improving the time complexity of Steps 6 and 7 of Algorithm B. It is always possible if we can find some properties about the vertices of P . For example, assume vi and vj are two vertices of Pl . If F (vi , Xr ) − F (vj , Xr ) can be found in constant time. Then, the time complexity of Step 6 can be improved to O(n + k). Another research direction is to modify the objective function, for example, using rectilinear link distance metric [1]. That is, find a point x in P , such that the summation of the weighted link distances from x to all points of X is maximized. Acknowledgments The authors wish to thank Prof. B.C. Liaw for his helpful discussion about this problem. We also thank the anonymous referees for the careful reading of the paper and give our constructive comments and suggestions which not only greatly improve our presentation of the results but also lead to the correction of some mistakes. References [1] M. de Berg, On rectilinear link distance, Comput. Geom.: Theory Appl. 1 (1991) 13–34.
[2] M. de Berg, M. Van Kreveld, Rectilinear decompositions with low stabbing number, Inform. Process. Lett. 52 (1994) 215–221. [3] M. de Berg, M. Van Kreveld, B.J. Nilsson, M.H. Overmars, Finding shortest paths in the presence of orthogonal obstacles using a combined L1 and link metric, Lecture Notes Comput. Sci. 447 (1990) 211–224. [4] D.Z. Chen, K.S. Klenk, H.T. Tu, Shortest path queries among weighted obstacles in the rectilinear plane, SIAM J. Comput. 29 (2000) 1223–1246. [5] V. Chepoi, F. Dragan, Computing a median point of a simple rectilinear polygon, Inform. Process. Lett. 49 (1994) 281–285. [6] M. Colebrook, J. Guttierez, J. Sicilia, A new bound and an O(mn) algorithm for the undesirable 1-median problem (maxian) on networks, Comput. Oper. Res. 32 (2005) 309–325. [7] M.T. Goodrich, R. Tamassia, Dynamic ray shooting and shortest path via balanced geodesic triangulation, in: Proceedings of the Ninth Annual ACM Symposium on Computational Geometry, 1993, pp. 318–327. [8] C. Levcopoulos, Heuristics for minimum decompositions of polygons, Ph.D. Thesis, University of Linköping, Linköping, Sweden, 1987. [9] A. Maheshwari, J. Sack, Simple optimal algorithms for rectilinear link path and polygon separation problem, Parallel Process. Lett. 9 (1999) 31–42. [10] K.M. McDonald, J.G. Peters, Smallest paths in simple rectilinear polygons, IEEE Trans. Comput.-Aided Des. 11 (1992) 864–875. [11] J. O’Rourke, Art Gallery Theorems and Algorithm, Oxford University Press, Oxford, 1987. [12] S. Schuierer, An optimal data structure for shortest rectilinear path queries in a simple rectilinear polygon, Int. J. Computat. Geom. Appl. 6 (1996) 205–225. [13] S.S. Ting, A linear time algorithm for maxisum facility location on tree networks, Transport. Sci. 18 (1984) 76–84.