Approximating minimum bending energy path in a simple corridor

Computational Geometry 47 (2014) 349–366

Contents lists available at ScienceDirect

Computational Geometry: Theory and Applications www.elsevier.com/locate/comgeo

Approximating minimum bending energy path in a simple corridor ✩ Lei Xu ∗ , Jinhui Xu Department of Computer Science and Engineering, State University of New York at Buffalo, Buffalo, NY 14260, USA

a r t i c l e

i n f o

Article history: Received 6 August 2012 Accepted 4 September 2013 Available online 18 September 2013 Communicated by J. Mitchell Keywords: Curvature Bending energy Approximation algorithm Corridor Path

a b s t r a c t In this paper, we consider the problem of computing a minimum bending energy path (or MinBEP) in a simple corridor. Given a simple 2D corridor C bounded by straight line segments and arcs of radius 2r, the MinBEP problem is to compute a path P inside C and crossing two pre-specified points s and t located at each end of C so that the bending energy of P is minimized. For this problem, we first show how to lower bound the bending energy of an optimal curve with bounded curvature, and then use this lower bound to design a (1 +  )-approximation algorithm for this restricted version of the MinBEP problem. Our algorithm is based on a number of interesting geometric observations and approximation techniques on smooth curves, and can be easily implemented for practical purpose. It is the first algorithm with a guaranteed performance ratio for the MinBEP problem. © 2013 Elsevier B.V. All rights reserved.

1. Introduction In this paper, we consider the following minimum bending energy path (MinBEP) problem: Given a simple corridor C bounded by the straight line segments and arcs of radius 2r (a precise definition of simple corridor is given later) and two points s and t at each end of C , compute a minimum energy path P crossing s and t and traversing the corridor C . The energy of a smooth path P is measured by its bending energy E b . The bending energy was suggested by Bernoulli in 1738 and studied intensively for the case of planar curves by Euler [10]. Given a smooth parameterized planar  curve s, its bending energy is defined as the integral of the squared curvature over the length of the curve, i.e., E b = s κ 2 ds, where κ is the curvature of the curve at point s. Let s(t ) = (x(t ), y (t )) be the parametrization of an arbitrary point on a smooth curve. Then the curvature [26] of the curve at s(t ) is given by κ (t ) = x˙(˙(xt()ty¨)2(t+)−y˙ (y˙t()t2)¨)x3(/t2) , where x˙ (t ) is the first derivative of x(t ) and y¨ (t ) is the secondary derivative of y (t ). Intuitively, the

curvature of a smooth curve at point s(t ) is defined as the inverse of its minimum curvature radius (i.e., the radius of the osculating circle at s(t )). The MinBEP problem considered in this paper is motivated by applications in interventional procedures for cardiovascular surgery [17,24,25,27]. In such applications, an elastic wire (called guidewire) is often used to guide a catheter/stent through the blood vessels to the narrowed (stenosed) spot to open it up. To optimize the interventional procedure (e.g., selecting the best fitting stent and catheter), it is often desirable to compute the path of the guidewire in advance. To solve this problem, one way is to view the blood vessel as a narrow corridor and determine the path of a stabilized guidewire by computing

✩ The research of this work was supported in part by National Science Foundation (NSF) through a CAREER award CCF-0546509, two grants IIS-0713489 and IIS-1115220. Corresponding author. Tel.: +1 (858) 886 6496; fax: +1 (716) 645 3464. E-mail addresses: [email protected] (L. Xu), [email protected] (J. Xu).

*

0925-7721/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.comgeo.2013.09.001

350

L. Xu, J. Xu / Computational Geometry 47 (2014) 349–366

the minimum bending energy path inside the corridor. Thus the guidewire problem can be formulated as a MinBEP problem in 3D. In this paper we consider the 2D case. The framework of our approach is potentially applicable to the 3D case. Due to the wide use of elastic curves in many fields, the minimum energy path problem has been extensively studied in the past. Most of the studies have focused on computing closed or knotted elastic curves with prescribed boundary conditions and fixed length [12,16,20,21]. Almost all of them are based on numerical techniques (e.g., partial differential equations). For unknotted curves, several results are known for solving the guidewire problem [17,24,25,27]. Among them, the works in [17,27] are based on finite element methods (FEM) and thus converge rather slowly. The work in [24,25] introduces a heuristic algorithm to approximate the path of a guidewire. It demonstrates much improved time efficiency (e.g., seconds vs. hours or even days needed by the FEM approaches). However, due to its heuristic nature, it does not have a provable quality guarantee and the resulting path is in general not smooth. Bending energy has also been used in computational biology for understanding the cellular processes [15]. In addition, in shape completion which is an intriguing problem in geometry processing with applications in CAD and graphics, a 3D Euler spiral — the curve having both its curvature and torsion evolve linearly along the curve is proposed and proved to be existent and unique up to a rigid transformation [14]. A problem closely related to the MinBEP problem is the curvature-bounded shortest path (also called minimum cost path) problem. It has been frequently used for motion planning in robotics and artificial intelligence and extensively studied [1–3,5,6,9,11,19]. In [22], it has been shown that in arbitrary polygons, the shortest path problem with constrained curvature is NP-Hard. An exponential time and space algorithm for generating feasible solutions is given in [11]. More efficient approximation algorithms were also obtained in [4,18,28]. In this paper, we consider the problem of approximating the 2D MinBEP problem. To fix the problems associated with existing solutions (e.g., [17,24,25,27,29]), we expect that the computed path (1) has a bending energy close to optimal, (2) is smooth at every point, and (3) can be computed efficiently. To achieve these goals, our strategy is to first solve a constrained version of the problem. In this constrained version, we require that the to-be-computed path P have a bounded curvature κmax depending on the shape of the corridor. Once we obtain solutions to the curvature bounded MinBEP problem, we can then perform a binary search on the maximum curvature to find the best solution to the MinBEP problem. Thus our focus is on the curvature bounded minimum bending energy path problem. We note that even this constrained MinBEP problem is extremely challenging to solve optimally for several reasons. Firstly, the curvature of an optimal curve could change continuously between 0 and κmax along the corridor, which seems to suggest that the problem may not even belong to the class of NP. Secondly, the optimal path could appear at any place inside the continuous free space of the corridor, making it rather difficult to predict its shape and position. Thirdly, due to the quadratic nature of the objective function, the solution to any subpath depends on the entire unknown optimal path and cannot be determined locally. Thus our goal for the MinBEP problem is to obtain a good approximation solution. To deal with the aforementioned difficulties, we have to be able to estimate the bending energy of the optimal path, predict how the energy is distributed along the corridor, and ensure the quality of the whole computed path. For these purposes, we first use a set of recursive shortest paths and local features of the corridor to establish a lower bound on the energy of an optimal solution as well as the possible energy distribution. With this lower bound, we then place various types of grid points to discretize the continuous space and guess the positions of the optimal path. With the help of a number of interesting observations on the minimum bending energy path, sophisticated approximation techniques, and involving analysis on the quality of the path, we show that it is possible to obtain a (1 +  )-approximation from the discretized space. Our algorithm runs in nearly quadratic time and can be easily implemented for practical purpose. To our best knowledge, this is the first quality guaranteed solution to the MinBEP problem. 2. Preliminaries Let D be a disk of radius r. A pipe is a rectangle with a fixed width of 2r and a non-fixed length e. An elbow is a sector with a fixed radius of 2r and a non-fixed center angle of π − θ for some 0  θ  π . A corridor is a set of alternatively appearing pipes and elbows concatenated in a way that the radius of each elbow completely overlaps with the side of a neighboring pipe with length 2r. A corridor C is simple if every pipe or elbow does not intersect the interior of any other pipe or elbow, and non-simple otherwise (see Fig. 1). Clearly, there is no hole in a such defined simple corridor. Each elbow together with its two adjacent pipes forms an elbow region ER (see Fig. 2). Let η = { v 1 , v 2 , . . . , v n } be the centerline of C . It is easy to see that η is a simple 2D curve, where the edge between every pair of consecutive vertices v i and v i +1 is either a straight line segment or an arc of radius r which alternatively appears along η . Thus C can also be viewed as the Minkowski sum η ⊕ D of η and D. In geometry, the Minkowski sum (also known as dilation) of two sets of position vectors A and B in Euclidean space is formed by adding each vector in A to each vector in B, i.e., the set A + B = {a + b | a ∈ A , b ∈ B }. The corridor C is bounded by straight line segments and arcs. Let C u and C l be the upper and lower boundary of C , separated from the two ends of C . (Assume that the segment v 1 v n is horizontally oriented.) Let {u 1 , u 2 , . . . , un } and { w 1 , w 2 , . . . , w m } be the vertices of C u and C l respectively. Thus the edges between consecutive vertices in C u and C l are either segments or arcs. It is also easy to see that for each vertex in the boundary of C , if it is adjacent to two segments, then it is a reflex vertex (i.e., the inner angle is larger than π ).

L. Xu, J. Xu / Computational Geometry 47 (2014) 349–366

351

Fig. 1. (a) Simple corridor C ; (b) A non-simple corridor C  .

Fig. 2. An elbow region ER.

Fig. 3. Illustration of convert corridor to simple polygon.

3. Lower bound of the MinBEP To design a quality guaranteed algorithm for the MinBEP problem, we need to know the minimum energy of an optimal solution, i.e., the lower bound of the MinBEP path. For this purpose, we observe that the bending energy of a smooth curve is the integral of the squared curvature along the curve. This suggests us to lower bound (1) the arc length of the portions of an optimal curve where the curvature is non-zero, (2) the curvatures at those portions. We first consider the problem of lower bounding (1). To estimate the curve length of (1), we consider the locations of the corridor where any curve has to make turns, since the curvature at those places cannot be zero. To find out those places, we consider the shortest path inside the corridor C . From [13], we know that the shortest path inside a simple polygon can be computed in linear time. (Note that although the corridor C is not a simple polygon, we can easily convert it into a simple polygon by replacing each arc with a convex polygonal curve around the arc (see dashed line in Fig. 3).) Let SP be the shortest path inside C with starting point s and ending point t. We assume that s and t are not visible to each other in C (since otherwise the segment st will be the trivial shortest path and also the minimum bending energy path in C ). Clearly, SP consists of a set of straight line segments and has unbounded curvature. Let V sp = {s1 , s2 , . . . , sk } be the set of interior vertices of SP (i.e., excluding s which is s0 and t which is sk+1 ) with the indices following their orders along SP. Let V´ be the set of reflex vertices in the boundary ∂ C of the corridor C . Clearly the boundary curve ∂ C is nondifferentiable at any point in V´ . The following lemma shows that the shortest path SP makes turns only at reflex vertices. Lemma 1. V sp ⊆ V´ . Proof. We prove this lemma by contradiction. Suppose the lemma is not true. Then there must exist at least one point in V sp but not in V´ . Let si be such a point (i.e., si ∈ V sp and si ∈ / V´ ). Let si −1 and si +1 be the predecessor and successor of si in SP respectively. Consider the triangle si −1 si si +1 . If ∂ C touches si −1 si from the inside of si −1 si si +1 , then let si −1 be the touching point in si −1 si closest to si . Otherwise, let si −1 be si −1 . Similarly, we define si +1 as the closest touching point (from the inside of si −1 si si +1 ) in si si +1 to si or si +1 if ∂ C does not touch si si +1 from the inside of si −1 si si +1 (see Fig. 4). Since si is not a reflex vertex of ∂ C and si −1 , si and si +1 are not collinear, neither si −1 nor si +1 can be coincident with si . This is because the inner boundary of the corridor consists of straight line segments and arcs. Thus if si −1 or si +1 (say si −1 ) is coincident with si , then si −1 must be an interior point of a segment or an arc due to the fact that si is not a reflex vertex of ∂ C . Since si −1 touches si −1 si from the inside of si −1 si si +1 , if si −1 is an interior point of a segment e, then e must go

352

L. Xu, J. Xu / Computational Geometry 47 (2014) 349–366

Fig. 4. Illustration of Lemma 1.

Fig. 5. Illustration of Lemma 2.

outside of si −1 si si +1 . This implies that si si +1 is outside of the corridor C , which contradicts the fact that si si +1 is inside of C . If si −1 is a point of an arc e  , since si −1 touches si −1 si from the inside of si −1 si si +1 , part of e  must be outside of si −1 si si +1 . This also implies that si si +1 is outside of the corridor C . Thus a contradiction. Since neither si −1 nor si +1 is coincident with si , there must exist a shortest path between si −1 and si +1 which is fully contained inside si −1 si si +1 and has shorter length than the subpath of SP between si −1 and si +1 . Thus we can replace the subpath of SP by the shorter path to obtain a shorter path between s and t, which is a contradiction to the fact that SP is the shortest path inside C . 2 Given any pair of consecutive vertices si and si +1 on the shortest path SP, let si and si +1 be the two closest (to si si +1 ) intersection points of ∂ C and the supporting line of segment si si +1 , where si (or si +1 ) is the intersection point closer to si (or si +1 ) than si +1 (or si ) (see Fig. 5). Let M be an arbitrary curvature-bounded smooth curve traversing the corridor C with the same starting and ending points s and t as SP, and κmax be its maximum curvature. Let M (t ) be the parametrization of M with time t as the parameter, where t 1  t  t 2 and s = M (t 1 ) and t = M (t 2 ) (i.e., we view the curve M as the loci of a point moving with unit speed). The following two lemmas show how M intersects with si si +1 . Lemma 2.  M ∩ si si +1   2 for any 1 < i < k (i.e., M intersects si si +1 at least twice and each intersection could be a consecutive portion of M overlapping with si si +1 ).

Proof. For this lemma, we consider two cases, (a) si and si +1 are on the same side of ∂ C (see Fig. 5) and (b) si and si +1 are on opposite sides of ∂ C . Below we will focus on case (a), as case (b) is an easier case and can be proved similarly. From Lemma 1, we know that both si and si +1 are vertices in V´ . Thus segment si si +1 partitions the corridor C into several (at

least 3 and at most 4) connected regions. Let C i be the region containing si −1 si and bounded by si si , and C i+1 be the region containing si +1 si +2 and bounded by si +1 si +1 . Segment si si +1 partitions the remaining region into possibly two connected

regions. Let C i,i +1 be the region bounded by si si +1 and the portion of ∂ C between si and si +1 (see Fig. 5), and C i ,i +1 be the possibly empty region bounded by si si +1 and the portion of ∂ C between si and si +1 . Clearly, C = C i ∪ C i+1 ∪ C i,i +1 ∪ C i ,i +1 , where C i and C i,i +1 are separated by segment si si , C i,i +1 and C i+1 are separated by segment si +1 si +1 , and C i,i +1 and C i ,i +1 are separated by segment si si +1 . Since SP is a shortest path inside C and SP makes turns at si and si +1 respectively, s and si +1 are not visible to each other (otherwise SP can be shortened by replacing the subpath of SP from s to si +1 with segment ssi +1 ). Similarly, we know that si and t are not visible to each other. Thus segments ssi +1 and si t are not entirely contained inside C . Also by the fact that SP is the shortest path in C , s must be in C i and t must be in C i+1 (otherwise SP does not need to make turns at si and si +1 ). To see this, consider the case of t. Suppose t is not in C i+1 . Then t cannot be in C i ,i +1 , since otherwise SP has a

L. Xu, J. Xu / Computational Geometry 47 (2014) 349–366

353

self-intersection at si si +1 . Also t cannot be in C i (since otherwise SP does not need to cross si si and visit si +1 ). Thus t can only be in C i,i +1 . Since si +2 ∈ C i+1 , there must exist a vertex s j ∈ V sp with j > i + 2 such that s j , s j +1 , . . . , t are all in C i,i +1 .

This implies that the subpath of SP between si +2 and s j intersects segment si +1 si +1 at a point, say v. Thus, if we replace the subpath of SP between si +1 and v by the segment si +1 v, we obtain a shorter path. This contradicts the fact that SP is the shortest path. Thus t must be in C i+1 . Similarly, we can show that s ∈ C i . Since M is a smooth curve in C with starting point s and ending point t, it must enter C i,i +1 (from C i ) through si si and

enter C i+1 (from C i,i +1 ) through si +1 si +1 . Thus  M ∩ si si +1   2.

2

Similar results hold for ss1 and sk t (i.e., for the cases of i = 1 and i = k), where s1 and sk are the first and last interior vertices of SP. Lemma 3.  M ∩ ss1   2 and  M ∩ sk t   2. Proof. Since the roles of s and t are symmetric, we only need to prove  M ∩ ss1   2.

First since s ∈ M,  M ∩ ss1   1. It is sufficient to show that M intersects the segment of ss1 at least once at a position

other than s. Clearly, ss1 partitions C into three connected regions C 1 , C 2 and C 3 , where C 1 and C 2 are separated by

segment ss1 , and C 2 and C 3 are separated by segment s1 s1 . By an argument similar to the one in the proof of Lemma 2,

we know that t must be in C 3 . Thus M must cross s1 s1 at least once, in order to enter from C 2 to C 3 . Hence the lemma follows. 2 − −−−− → − −−−− →

− −−−− →

− −−−− →

Let si ∈ V sp be any interior vertex of SP, and αi = si −1 si , si si +1 be the minor angle between vectors si −1 si and si si +1 . − −−−− → − −−−− → Clearly, si −1 si , si si +1 = π − si −1 si si +1 . Let r (si , si −1 ) be the ray emitting from si and crossing si −1 , and si −1 be the first intersection point of r (si , si −1 ) ∩ ∂ C after passing si −1 . Similarly, we have r (si , si +1 ) and si +1 . The two rays r (si , si −1 ) and r (si , si +1 ), together with the portion of ∂ C on the opposite side of si form a closed region in C . We denote it as C i , and call it the turning region of si . Let e i ,i −1 be the segment of si si −1 if si and si −1 are on the same side of ∂ C , and si si −1 otherwise. Similarly, we have e i ,i +1 . Clearly, e i ,i −1 and e i ,i +1 separate the turning region C i from the remaining part of C . ˙ (t x ) be the unit tangent vector of M at point x (i.e.,  T x  = ˙x = 1), Let x be any point on M with x = M (t x ), T x = x˙ = M and κx be the curvature of M at x. To simplify our discussion, from now on we assume every vector is a unit vector (unless we specify otherwise). Below we show that for any smooth curve M, its accumulated total angle change around each turn of SP is at least αi . The following two key lemmas are for the case in which the shortest path has at least three interior vertices. Lemma 4. Let M be any planar smooth curve with non-zero curve length and x1 x2 x3 be any triangle with x1 ∈ M. Then there exists a point y in M but not in x1 x2 x3 . Proof. Suppose that there is no such a point in M. Let B x1 , be a ball centered at x1 and with an arbitrarily small radius  > 0. By the assumption, we know that M ∩ B x1 , ⊆ x1 x2 x3 . Let M (t ) be the parametrization of M with parameter t ˙ (t ), and T (x1 )+ = limt →t x ;t >t x M ˙ (t ). Since M ∩ B x1 , is inside and x1 = M (t x1 ) for some t x1 . Let T (x1 )− = limt →t x ;t
1

1

1

−−− → − −−− → x1 x2 x3 , we have  T x−1 , T x+1  − x2 x1 , x1 x3 . This implies that T x−1 =

T x+1 , i.e., T x1 does not exist. This contradicts the fact that M (t ) is a smooth curve and differentiable at any interior point. 2

Lemma 5. Let si ∈ V sp be any interior vertex of SP and M be any smooth curve in C (between s and t). There exist two points m1 and m2 in M ∩ C i so that  T m1 , T m2 > αi , where T m1 and T m2 are the unit tangent vectors of m1 and m2 respectively. Proof. Let v j and v j +1 be the two vertices in the centerline η closest to si , and v i be the intersection point of the two lines supporting segments v j −1 v j and v j +1 v j +2 respectively. Let l(si , v i ) be the line crossing points si and v i . It is easy to see that l(si , v i ) partitions C into two regions at si , with one containing s, and the other containing t. Thus any smooth curve M ∈ C between s and t intersects l(si , v i ) (see Fig. 6). Since M ∩ l(si , v i ) = Φ , we consider the following two cases: A. si ∈ / M ∩ l(si , v i ); B. si ∈ M ∩ l(si , v i ). For the two cases, first let x = si be any point in si −1 si and y = si be any point in si si +1 . Since si −1 si and si si +1 are edges of the shortest path SP and si is a reflex point of ∂ C , it is easy to see that segment xy is not fully inside C and si v i ∩ {si −1 si si +1 } = {si }. Then we analyze each of the two cases.

354

L. Xu, J. Xu / Computational Geometry 47 (2014) 349–366

Fig. 6. Illustration of Lemma 5.

A: For case A, since l(si , v i ) ∩ M = φ , there must exist a point ξ = si and ξ ∈ l(si , v i ) ∩ M. From Lemma 2, we know that

M intersects l(si −1 , si ) twice with one at each side of si (i.e., at si −1 si −1 and si si respectively). Similarly, M intersects l(si , si +1 ) at both sides of si . Thus there must exist two points x1 ∈ M ∩ l(si −1 , si ) and x2 ∈ M ∩ l(si , si +1 ) so that M passes through x1 , ξ, x2 in this order. Without loss of generality, we assume that x1 is the last intersection point of M and l(si −1 , si ) before M reaches ξ , and x2 is the first intersection point of M and l(si , si +1 ) after M leaves ξ . Then we can build two lines l(x1 , ξ ) and l(ξ, x2 ). Since M is smooth and second order differentiable between x1 and ξ , by Mean Value Theorem,1 we know that there exists a point m1 ∈ M between x1 and ξ such that T m1 =

M (t ξ )− M (t x1 ) t ξ −t x1

− −− →

= x1 ξ , where M (t x1 ) = x1 , M (t ξ ) = ξ , and − −− → M (tm1 ) = m1 . Similarly, there exists an m2 ∈ M between ξ and x2 such that T m2 = ξ x2 . − −− → − −− → − −−−− → − −−−− → Now consider the triangles x1 si x2 and x1 ξ x2 . Since x1 ∈ l(si −1 , si ) and x2 ∈ l(si , si +1 ), x1 si , si x2 = si −1 si , si si +1 − −− → − −− → − −− → − −− → = αi . Also since x1 ξ x2 < x1 si x2 , x1 ξ , ξ x2 = π − x1 ξ x2 , and αi = π − x1 si x2 , we have αi < x1 ξ , ξ x2 . This means  T m1 , T m2 > αi . B: For case B, we know si ∈ M ∩ l(si , v i ). There must exist an arbitrarily small constant  > 0 so that in the  -neighborhood of 1. 2. 3. 1:

si (i.e., the ball B si , centered at si and with radius  ) exactly one of the following cases occurs. M ∩ B si , ∩ si −1 si = {si }, and M ∩ B si , ∩ si si +1 = {si }; M ∩ B si , ∩ si −1 si ⊃ {si }; M ∩ B si , ∩ si si +1 ⊃ {si }. This case means that when M enters the  -neighborhood of si (i.e., B si , ), M intersects si −1 si and si si +1 only at si . Let M  = M ∩ B si , . Removing si from M  will result in two smooth sub-curves, M − and M + , where M − is the portion of M  before si and M + is the portion of M  after si . Based on the position of the two sub-curves, we have three cases. 1.1. Both M − and M + are outside of si −1 si si +1 ; 1.2. M − is outside of si −1 si si +1 and M + is inside of si −1 si si +1 ; 1.3. M − is inside of si −1 si si +1 and M + is outside of si −1 si si +1 . Note that by Lemma 4, we know that M − and M + cannot be both inside of si −1 si si +1 . Below we first consider case 1.1. For case 1.1, by a similar argument as in case A, we know that there exist x1 ∈ l(si −1 , si ) ∩ M and x2 ∈ l(si , si +1 ) ∩ M so that M passes through x1 , si , x2 in this order. Without loss of generality, we assume that x1 and x2 are such intersection points closest to si . Since both M − and M + are outside of si −1 si si +1 , there must exist an ξ ∈ l(si , v i ) such that the two lines l(ξ, x1 ) and l(ξ, x2 ) both intersect M \ {si }. To see this, we can imagine there is a point ξ continuously moving from si towards v i . Initially, l(ξ, x1 ) and l(ξ, x2 ) both intersect M at si . When ξ moves away from si by an arbitrarily small distance, l(ξ, x1 ) and l(ξ, x2 ) will still intersect M − and M + respectively. This is because both M − and M + are continuous curves and are outside of si −1 si si +1 . Let y 1 and y 2 be the two − −− → intersection points. By using Mean Value Theorem, we know that there exist m1 ∈ M and m2 ∈ M such that T m1 = x1 ξ − −− → − −− → − −− → − − − − − → − − − − − → and T m2 = ξ x2 . Considering the triangles x1 si x2 and x1 ξ x2 , we know that x1 ξ , ξ x2 > si −1 si , si si +1 . Thus,

 T m , T m > αi . 1 2

The Lagrange Derivative Mean Value Theorem or Mean Value Theorem says: If f (x) is a continuous function in interval [a, b] and has the first order f (b)− f (a) . derivative ˙f (x) for a < x < b, then ∃ξ , a < ξ < b, so that ˙f (ξ ) = b−a 1

L. Xu, J. Xu / Computational Geometry 47 (2014) 349–366

355

For cases 1.2 and 1.3, since they are symmetric, we only consider case 1.2. For this case, we first consider M + and

the line l(si , si +1 ). By Lemma 2, we know that M intersects segment si +1 si +1 . Also it is possible for M to intersect

segment si si +1 . Let x2 be the first intersection point of M and si si +1 after M leaves si . Since si , x2 ∈ M, by Mean − −− → Value Theorem we know that there exists m2 ∈ M such that T m2 = si x2 . − For M  , since it is outside of si −1 si si +1 , by the same argument given in case 1.1, we know that there ex− −− → ist ξ ∈ l(si , v i ) and m1 ∈ M such that T m1 = x1 ξ . Let ξ  be the intersection point of l(x1 , ξ ) and l(si , si +1 ) (note −−−→ −−−→

− −−−− → − −−−− →

that ξ  ∈ si si ). Considering the triangles x1 si x2 and x1 ξ  x2 , we know that x1 ξ  , ξ  x2 > si −1 si , si si +1 . Thus,

 T m , T m > αi . 1 2 2: Let M  = B si , ∩ M, M − be the portion of M  before si , and M + be the portion of M  after si . Since  is arbitrarily small, case 2 means that after M enters the  -neighborhood of si (i.e., B si , ), it “stays” on segment si −1 si until it reaches si (i.e., M − ⊆ si −1 si ). − −−−− → Let m1 be any point in M − . Clearly, its tangent vector T m1 = si −1 si . For M + , since it is outside of si −1 si si +1 , by the same argument given in case 1.1, we know that there exist ξ ∈ l(si , v i ) and x2 ∈ si si +1 such that l(ξ, x2 ) intersects − −− → M + . By Mean Value Theorem, we know that there exists an m2 ∈ M such that T m2 = ξ x2 . It is easy to see that − − − − − → − − − − − →

s i −1 s i , s i s i +1 <  T m , T m . Thus we have  T m , T m > αi . 1 2 1 2 3: Similar to case 2. Thus for all subcases in case B, we have  T m1 , T m2 > αi . Hence the lemma follows.

2

Next we consider the case that SP has only one or two interior points, i.e., V sp = {s1 } or V sp = {s1 , s2 }. Lemma 6. If V sp = {s1 } or V sp = {s1 , s2 }, then ∃m1 ∈ M, and ∃m2 ∈ M, such that αi <  T m1 , T m2 for i ∈ {1, 2}. Proof. The proof of this lemma is similar to that of Lemma 5. The only difference is that in these two cases, we need to use Lemma 3 to replace Lemma 2 in some part of the proof. 2 Corollary 1. Let si ∈ V sp be any interior vertex of SP, and M be any smooth curve in C i intersecting both e i ,i −1 and e i ,i +1 at a position other than si . Then there exist two points m1 and m2 in M ∩ C i such that  T m1 , T m2 > αi . The above two lemmas and corollary suggest that the angle change by a smooth curve M around any interior point si of the shortest path SP is at least as large as that of SP (i.e., αi ). The following lemma lower bounds the arc length of M for making such an angle change. Lemma 7. Let M be a smooth curve in C between s and t with maximum curvature of α with non-zero curvature in C i is at least κ i .

κmax . Then the arc length of the portion of M

max

Proof. Since the curvature of M is no more than κmax , the radius of the osculating circle at any point of M is at least κ 1 . max By Lemmas 5 and 6, we know that the angle change of M in C i is at least αi . Thus, the minimum arc length for M to make α an angle change of αi is at least κ i . 2 max

The above lemma gives a lower bound on the arc length of an optimal curve P (which is a smooth curve) in the turning region C i with non-zero curvature. Let E i be the bending energy of P in C i . Then to lower bound E i , we need to determine the curvature of this portion of P . Lemma 8. Among all smooth curves with an accumulated angle change of same curvature at every point.

αi , the curve with the minimum bending energy has the

Proof. Consider a discretization of the curve M where each consecutive pair of discrete points in M has length s. Since s is infinitely close to 0, we can assume that the curvature of M between any pair of consecutive points is the same at every point. This means that the curve between consecutive pair of points is an arc. Let κ j be the curvature of M between the j-th and ( j + 1)-th points and θ j be the angle of the corresponding arc. Thus the bending energy of M can be written





θ



j 2 2 as E i = j κ j s = j ( s )  s = equal. Thus the lemma follows. 2

θ 2j

j s ,

where



j

θ j = αi . It is easy to see that E i achieves its minimum when all θ j are

The above lemma implies that locally the minimum bending energy is achieved when the curve is an arc. Thus to find the minimum curvature κi , we only need to determine the radius r i of the maximum osculating circle which has a circular arc B max in C i of angle αi . Below we discuss our ideas on determining r i .

356

L. Xu, J. Xu / Computational Geometry 47 (2014) 349–366

Fig. 7. Four types of C i depending on e i ,i −1 and e i ,i +1 .

Fig. 8. An example of cases A.1 and A.4.

Definition 1. An elbow in C is a short elbow if SP does not touch its reflex vertex u. An elbow region containing a short elbow is called a short elbow region. We have two cases about the turning region C i : (A) C i contains no short elbow, and (B) C i contains at least one short elbow. We first consider case (A). As discussed above, C i is bounded by segment e i ,i −1 and e i ,i +1 . Let si −1 (si +1 ) be the farthest endpoint of e i ,i −1 (e i ,i +1 ) away from si . Clearly, si −1 (si +1 ) is same as si −1 (si +1 ) if si and si −1 (si +1 ) are on the same side of ∂ C and si −1 (si +1 ) otherwise. Thus we have four cases of C i depending on whether si −1 = si −1 and si +1 = si +1 (see Fig. 7). Let ERi be the elbow region in C i with reflex vertex u = si and an outer angle of θi (see Fig. 2). Let uu 1 and uu 2 be the two incident edges of u, e 1 and e 2 be the corresponding edges of uu 1 and uu 2 on the opposite side of ∂ C , and u 1 and u 2 be the corresponding vertices of u 1 and u 2 . Let |e i |, i ∈ {1, 2}, denote the length of e i , and w denote the width of C . Clearly |uu 1 | = |e 1 |, |uu 2 | = |e 2 |, and w = 2r. To determine r i , we need to first find B max in C i . Based on the relationships of B max and ERi , we have four subcases to consider (see Figs. 8 and 9). A.1. B max is tangent to both e 1 and e 2 at their interior points. B max starts and ends at the two tangent points. In this case, si is incident to B max due to the fact that B max is the largest circular arc, the following conditions are satisfied. θi θi θi w w θi cos 2  |e 1 |, and θi cos 2  |e 2 |. Since r i = r i sin 2 + w, we have 1−sin

1−sin

2

ri =

w 1 − sin

θi

.

2

(1)

2

A.2. B max is tangent to e 1 only. In this case, si and si +1 are incident to B max . B max starts at the tangent point of e 1 and ends at si +1 .

L. Xu, J. Xu / Computational Geometry 47 (2014) 349–366

357

Fig. 9. An example of the two subcases of case A.3.

A.3. B max is tangent to e 2 only. This is the symmetric case of case 2. In this case, si and si −1 are incident to B max . B max starts at si −1 and ends at the tangent point of e 2 . A.4. B max is tangent to none of e 1 and e 2 . This implies that si , si −1 and si +1 are all incident to B max , and B max starts at si −1 and ends at si +1 . Thus we have

1

ri =

2 sin θi

 |e i ,i −1 |2 + |e i ,i +1 |2 − 2|e i ,i −1 ||e i ,i +1 | cos θi .

(2)

Since our goal is to lower bound E i and cos θi > −1 (i.e., 0 < θi < π ) in Eq. (2), we can choose

1

ri =

2 sin θi

  |e i ,i −1 | + |e i ,i +1 | .

(3)

In each case of C i , B max is always inside C i and the accumulative angle change of B max is following lemma follows directly from the above analysis.

αi (by Corollary 1). The

Lemma 9. In case (A), B max has only the above 4 cases and r i is a function of the local features (i.e., w , θi , |e 1 |, e 2 |, |e i ,i −1 | and

|e i ,i +1 |) of C i .

Proof. We have 4 cases depending on the local features of ERi in C i . 1: |e 1 |  2: |e 1 | < 3: |e 1 | <

w 1−sin w

θi

1−sin w

θi

cos

2

cos

2

θ

cos

θi 2

θi 2

θi

, and |e 2 |  , and |e 2 | < , and |e 2 | >

w 1−sin w

θi

1−sin w

θi

cos

θi

cos

θi

cos

θi

2

2

θ

2 2

. This is in case A.1 and we use Eq. (1) to compute r i . . This is in case A.4 and we could use Eq. (3) to compute r i . . We have two subcases here.

2 2 1−sin 2i 1−sin 2i 3.1: si −1 and si are on the opposite sides of ∂ C (see Fig. 9(a)). Let |dmax | be the distance between the center of B max and e i ,i −1 and θi = θi + arcsin |e w | . We have i ,i −1



|dmax | =

X 1 − (tan4 θi + 2 tan2 θi )( X 2 + X 3 ) − X 4 tan4 θi + 2 tan2 θi

where X 1 = (tan3 θi |e i ,i −1 |

2

|e i ,i −1 | 2

+

w tan θi (tan2 θi +1) sin θi

(1 + tan2 θi )2 , and X 4 = (tan3 θi

4

Hence,



ri =

|dmax |2 + √

|e i ,i −1 |2 4

|e i ,i −1 | 2

(4)

,

+ tan θi )2 , X 2 = (tan2 θi +

w tan θi (tan2 θi +1) sin θ  i

|e i ,i −1 | 2

+

w tan2 θi 2 ) , sin θi

X 3 = (tan θi +

w tan θi 2 ) sin θi

−

+ tan θi ).

(5)

.

Thus if |e i ,i +1 | > 2r w − w 2 , it is in case A.3. Otherwise it is in case A.4. 3.2: si −1 and si are on the same side of ∂ C (see Fig. 9(b)). In this case, |dmax | has the same form as shown in Eq. (4),

except for replacing θi with θi . Thus r i has the same form as in Eq. (5). Hence if |e i ,i +1 | > case A.3. Otherwise it is in case A.4. θi θi w w 4: |e 1 | > 2 θi cos 2 , and |e 2 | < θi cos 2 . It is symmetric to the case 3. 1−sin

2

1−sin

√

2r w − w 2 , it is in

2

Note that we could compute r i faithfully following the description in cases A.2, A.3. However, since our goal is to obtain a lower bound on E i , it is no harm if we have a larger r i of B max . In this way, the lower bound E li of E i will be smaller,

358

L. Xu, J. Xu / Computational Geometry 47 (2014) 349–366

Fig. 10. C i contains one short elbow.

and consequently our approximation solution (given in the next section) will be closer to the optimal. The following lemma give an alternative way for determining r i in cases A.2 and A.3. Lemma 10. r i in case A.3 (or A.2) can be computed by using either Eq. (1) as in case A.1 or Eq. (3) as in case A.4. Proof. Since case A.2 is the symmetric case of case A.3, we only consider the computation of r i in case A.3. We have two subcases, depending on whether si and si −1 are on the same side of ∂ C (see Fig. 9). If si and si −1 are on the same side of ∂ C , r i can be determined by using Eq. (5) after replacing θi with θi in Eq. (4). Such obtained r i has a smaller value than

using Eq. (1) in case A.1. Thus we can use Eq. (1) to compute r i to obtain a larger r i (i.e., a lower E li ). Otherwise (si and si −1 are on different sides of ∂ C ), B max is a circular arc decided by {si , si , st } where st is the tangent point of e 2 . This is similar 1 to case A.4. Thus we can use Eq. (3) as in case A.4 to compute r i . 2 Note that if C i contains s (or t), r i can still be computed in a similar way. We only need to replace si −1 by s (or si +1 by t). It is easy to see that in each case, B max can always be determined by the local features of C i . From Lemmas 7, 9, and 10, we immediately have the following lemma on the lower bound of the bending energy of the optimal curve P in C i . Lemma 11. The optimal curve P in the turning region of C i in case (A) has a bending energy

E i > E li =

αi 2 α 1 κi = i 2 , κmax κmax ri

(6)

where E li is the lower bound of E i . Now we consider case (B) (i.e., C i contains at least one short elbow). In this case, it is difficult to compute r i since C i could be quite complicated, and the local features of ERi may not be sufficient to determine r i . We consider two subcases: (B.1) C i contains exactly one short elbow, and (B.2) C i contains more than one short elbow. B.1: Consider the case in Fig. 10 (other cases can be handled similarly). Let ERk be the only short elbow region in C i with reflex vertex k and k be the middle point of the arc of ERk . Obviously, to determine a good lower bound of the bending energy of P in C i , we need to take into consideration of all the local features of C i which could be quite complicated. To better solve this problem, we first partition C i , as well as C i +1 , into two parts by introducing the segment kk . Let C i and C i+1 be the adjusted turning regions of si and si +1 respectively. Obviously, kk intersects the optimal curve P . Let pk be the intersection point, and assume that the exact position of pk is known in advance (in the next section, we will show how to guess the positions of pk ). With pk , we can then compute the shortest paths SP si−1 , pk and SP pk ,si+1 between si −1 and pk and between pk and si +1 respectively. This will give us more accurate α angles αi , αi+1 at si and si +1 . Further, we can use a similar approach as in case (A) to derive better lower bounds of the bending energy for the portion of P in turning regions C i and C i+1 respectively by using SP si−1 , pk and SP pk ,si+1 . The only difference is that B max ends (or starts) at pk in C i (or C i+1 ). B.2: Let {ER1 , ER2 , . . . , ERm } be the sequence of short elbow regions in C i between si −1 and si +1 , with reflex vertices {k1 , k2 , . . . , km }. Assume that we know pm (same as pk in case B.1) in advance. We can then compute shortest paths SP si−1 , pm , SP pm ,si+1 , and SP si ,si+1 (if necessary), and adjust the turning regions for si −1 , si and si +1 . If there is some short elbow in any of the adjusted turning regions (based on the newly computed shortest paths), we can recursively apply the above strategy to adjust the turning regions and shortest paths until we eliminate all short elbows.

L. Xu, J. Xu / Computational Geometry 47 (2014) 349–366

359

Lemma 12. The bending energy of an optimal curve P is at least

1 k

El =

2

E li ,

(7)

i =1

where E li is the lower bound on the bending energy of P in C i . Proof. The coefficient 12 follows from the fact that two consecutive regions C i and C i +1 overlap each other and the bending energy in the overlapped part are counted twice in the summation. 2 4. Approximation algorithm for computing MinBEP In this section, we present a (1 +  )-approximation algorithm for the MinBEP problem with bounded curvature. Let E b be the bending energy of an optimal curve P and E a be the bending energy of the path P A computed by our algorithm. As mentioned in Section 1, there is a number of difficulties to compute the MinBEP. To overcome them, our main idea is to make use of the lower bound obtained in the last section to achieve a good approximation. From the lower bound discussion, we know that most of the bending energy of the optimal curve P is distributed in the turning regions of those reflex vertices. Thus, to approximate the optimal path P , P A has to be similar to P in those regions. To achieve this, we first compute the shortest path SP in C and then isolate the task of approximating the optimal curve within each turning region (based on SP). For every turning region, we first sample a set of points as the possible positions for the optimal curve to either pass through or appear in their neighborhoods. Then we connect those sample points by smooth curves to form a network or graph G so that there exists at least one smooth curve in G which is geographically close (i.e., has a small Hausdorff or Fréchet distance) to the optimal curve P and differs in bending energy from P by only a small constant factor. To make this approach work, we have to first solve the following problems: (1) How to sample the points? (2) How to connect the set of sample points? (3) How to ensure the existence of a smooth curve with similar bending energy? (4) How to find such a smooth curve and how to ensure that the time complexity is truly polynomial of the input size? To simplify our discussion, we first perform an isotropic scaling on the corridor C . Isotropic (or uniform) scaling is a linear transformation in Euclidean space that enlarges or shrinks objects with the same scale factor in all dimensions. Lemma 13. For any smooth curve P c with bending energy E c in the corridor C , after isotropic scaling with scale factor γ , the bending energy E c of the scaled curve P c is γ1 E c . Proof. Consider a discretization of P c with each consecutive pair of discrete points separated by a distance s infinitely close to 0. By a similar argument given in the proof of Lemma 8, we know that

Ec =



κi2 s =

 1

i

r i θi r i2

i

=

 θi i

ri

,

where r i is the radius of arc between the i-th and (i + 1)-th discrete points. After isotropic scaling, θi is still preserved. The radius r i , however, has been scaled by a factor of

E c =

 i

κi 2 s =

 i

1

(γ ri )2

γ r i θi =

1

γ

Ec.

γ . Thus we have

2

The above lemma suggests that isotropic scaling does not affect the performance ratio of an approximate path. Thus, to design a (1 +  )-approximation for the MinBEP problem, we can first choose an appropriate factor to scale the corridor C so that the width w becomes a large enough constant (depending on κmax ), and then compute a (1 +  )-approximating path P A in the scaled corridor. After obtaining P A , we can then scale the corridor together with the path back to its original size. The scaled path will still be a (1 +  )-approximating path for the original corridor. Thus from now on, we assume that both w and κmax are constants. Below we discuss our ideas on each of the aforementioned four problems. 4.1. Sampling points (Problem 1) To sample the possible positions of the optimal curve P , we have two sub-problems to consider: (i) How to sample points inside a turning region C i in case (A); (ii) How to sample points inside a turning region C i in case (B). To solve sub-problem (i), we first determine B max which could span over at most three elbow regions. Let a and b be the two endpoints of B max , and la and lb be the two lines which cross a and b respectively and are orthogonal to the center line η . We denote the region of C bounded by la and lb as ϕi (which could be either a subset or superset of C i , see Fig. 11).

360

L. Xu, J. Xu / Computational Geometry 47 (2014) 349–366

Fig. 11.

ϕi .

Fig. 12. Grid for ERi .

α

From Lemma 11, we know that the bending energy of C i is at least E li (i.e., κ i 12 ). The energy E li is estimated based on max r i the curvature of B max and αi (obtained from SP). Since the angle of B max is no smaller than αi , this implies that ϕi allows a smooth curve to complete an angle change of αi with the smallest curvature 1/r i . This indicates that it is sufficient to focus only on ϕi when sampling points for C i . To sample points in ϕi , consider the elbow region ERi in C i with reflex point u (i.e., si = u). Let ERi be the extended elbow region of ERi by enlarging or shortening the length of its two pipes so that ϕi is barely contained inside ERi . To simplify our discussion on sampling points for ϕi , we can first sample points in ERi and then throw away all points outside of ϕi to obtain the points for ϕi . Thus we will focus on ERi only. Our idea for sampling points in ERi is to place a bended grid in ERi so that there always exists a set of grid points close to P (see Fig. 12(b)). In a normal grid, the grid points are the intersections of a set of evenly spaced horizontal and vertical lines (see Fig. 12(a)). In an elbow region, however, since the corridor is bended, our grid bends around the elbow. The “horizontal” lines in our grid are evenly spaced curves “parallel” to the center line η of the corridor C . The “vertical” lines in the two pipes are evenly spaced lines orthogonal to the boundary edges of the pipes, and the “vertical” lines in the elbow are evenly spaced rays emitted from the reflex point u. The intersections of the “horizontal” and “vertical” lines forms the set of grid points. Let s be the distance between two neighboring “horizontal” lines and s be the maximum distance between two neighboring “vertical” lines. We choose s = κ 2 . For s, on one hand it should be small enough so that the computed max path P A is close enough to the optimal curve P . On the other hand, they should not be too small since otherwise the resulting graph would be very large, resulting in a slow algorithm. To determine the value for s, we only consider two types of C i in cases A.1 and A.4, since case A.2 or A.3 can be reduced to one of them. To avoid putting too many grid points in ERi (especially in a long pipe), we treat the pipes and elbow of ERi differently. More specifically, we let s = sp for a pipe of ERi and s = se for the elbow. Depending on the types of C i and the lower bound of the MinBEP problem discussed in the above section, we have the following choices for s. A.1: sp =

w cos

θi 2

κmax αi (1−sin

θi 3 )

A.4: We have two subcases: A.4.1: If |e i ,i −1 | + |e i ,i +1 |  and se =

w (π −θi ) . θ κmax αi (1−sin 2i )2

and se =

2

2w 1−sin

θi

w (π −θi ) κmax αi sin2 θi (1−sin

A.4.2: If |e i ,i −1 | + |e i ,i +1 | >

2w 1−sin

, then sp =

2

θi 2

θi 2 )

w

κmax αi sin2 θi (1−sin

θi 3 )

,

2

.

2

, then sp =

(|e i ,i −1 |+|e i ,i +1 |)3 and 8κmax αi sin2 θi w 2

se =

(π −θi )(|e i ,i −1 |+|e i ,i +1 |)2 . 4w κmax αi sin2 θi

L. Xu, J. Xu / Computational Geometry 47 (2014) 349–366

361

Now we consider sub-problem (ii) (i.e., sampling points for C i in case (B)). Note that for this problem, we cannot faithfully follow the idea for deriving lower bound for case B. This is because it requires us to recursively compute the shortest paths and each higher level shortest path requires the exact position of the optimal curve P in some short elbow which is in general unavailable. To overcome this difficulty, our idea is to guess the position of P by placing grid points in the last short elbow ERk in C i (see case B.2). More specifically, let kk be defined as in case B.1 in the last section (i.e., the segment connecting the reflex vertex and the middle point of the arc). We place a set of grid points evenly spaced in kk with a separating distance of 2 s = κmax . With the set of grid points, we can then recursively compute higher level shortest paths starting from (or ending at) each grid point to obtain the adjusted turning region as well as αi angle for each reflex point in the set of shortest paths. For an adjusted turning region with multiple αi values (from different grid points), we can use the same approach as in case (A) to obtain a set of densities (i.e., s) values of the grid and choose the densest one for its grid. Note that if we compute the shortest path between each pair of grid points, in the worst case, the number of shortest paths we need to compute is exponential in m, where m is the total number of short elbows in C i . Thus the total time for computing the shortest paths in all levels is exponential. On the other hand, since our goal is to decide the density (sp and se) of grid points in C i , there is no need to compute all levels of shortest paths. In case (B.1), each grid point on kk corresponds to a distinct set of three adjusted turning regions C i , C k and C i+1 as well as αi angle for each adjusted turning region. Note that (1) sampling grid points in each of the three turning regions is one subcase of case (A) and (2) αi always appears in the denominator of sp and se in case (A). Thus it is sufficient to find the grid point on kk that maximizes αi . (The maximum αi means the densest grid points for the corresponding turning region.) Thus we choose the extreme point on kk to compute αi . Once we have this information, it can be reduced to case (A). In case (B.2), for each short elbow in C i , we place kk with grid points as discussed above. For each kk , we compute the maximum αi as in case (B.1) and thus choose a set of extreme grid points on each kk . Then sampling grid points in each adjusted turning region is one subcase of case (A). In this way, we could sample points inside a turning region C i in case (B). For a turning region C i containing some short elbows, it is easy to see that the bending energy of P in the region spanned by those short elbow regions could be 0. To ensure that locally the bending energy of our computed path P A is within a (1 +  )-factor away from that of P in such a region, our algorithm has to be able to produce a 0-energy path if P has 0 bending energy. We adopt two strategies to ensure this, (1) adding inter-elbow grid points and (2) adding inter-elbow edges in the graph G (construction of graph G is given later). We discuss our ideas on (1) next and (2) in the next section. Now we discuss how to generate inter-elbow grid points. This type of grid points is necessary. Let si and si +1 be two consecutive interior vertices of the shortest path SP. Let {ER1 , ER2 , . . . , ER j } be the sequence of short elbow regions between si and si +1 , ϕi and ϕi +1 be the corresponding grid region of si and si +1 . If P has zero-curvature in those short elbow regions while the grid points are not dense enough in ϕi and ϕi +1 , it is possible that the grid points in ϕi and ϕi +1 are not visible to each other. Consequently, P A will be forced to use a small amount of energy to pass each short elbow region and the total energy could be relatively large. Thus, it is difficult to guarantee the performance of P A . Actually, this scenario may occur as long as any vertical line segment l v i in ϕi and any vertical line segment l v i+1 in ϕi +1 are weakly visible to each other, but there is no grid point visible to each other. (Two segments l v i and l v i+1 are weakly visible to each other if there exists a pair of visible points on l v i and l v i+1 respectively.) Moreover, increasing the density of ϕi and ϕi +1 does not solve the problem since the visible tunnel between l v i and l v i+1 could be very narrow (i.e., close to a line segment). To overcome this difficulty, our idea is to place inter-elbow grid points. More specifically, for each vertical line l v i in ϕi , we do two main computations, (1) the weakly visible portion between l v i and the farthest weakly visible vertical line lt i close to t and (2) the weakly visible portion between l v i and the farthest weakly visible vertical line l si close to s to check whether we need to put inter-elbow grid points (see Algorithm 1). Algorithm 1 Algorithm of putting inter-elbow grid points. for each vertical line l v i in ϕi do

compute the farthest weakly visible vertical line lt i close to t. The portion of e on l v i is weakly visible from the portion of e  on lt i if no dotted line between e and e  then if one grid point g i on e and another grid point gti on e  then continue else if no grid point on e and one grid point gti on e  then place a grid point g i on e else if no grid point on e  and one grid point g i on e then place a grid point gti on e  else if no grid points on e and e  then place a grid point g i on e and place another grid point gti on e  end if connect g i and gti by a dotted line else place the grid point g i on e as the intersection point between l v i and dotted line which touched e  end if Compute the farthest weakly visible vertical line l si close to s and connect it to l v i by a dotted line ld . If there is no gird point at the intersection between ld and l v i , place a gird point at the intersection end for

362

L. Xu, J. Xu / Computational Geometry 47 (2014) 349–366

Fig. 13. Weakly visible segments.

To implement the algorithm, we first compute each shortest path map (using the algorithm in [23]) of the endpoints of l v i and endpoints of each vertical line lk not in ϕi . For this purpose, we first convert C into a simple polygon C P by replacing each arc with a convex polygonal curve around the arc. It is easy to see that expanding the corridor in this way does not affect the shortest maps of the endpoints of l v i and endpoints of lk , since the shortest path between any pair of vertices will never go out of C . Then we determine the cusp vertices on the shortest path in constant time. Clipping l v i and lk by the common tangent lines defined by these cusp vertices gives the portion e and e  (in constant time). Thus if the total number of vertical lines is O (n), the worst running time of the above algorithm is O (n2 ). Moreover, since s is a constant, each vertical line l v i could be charged at most two inter-elbow grid points by simple amortized analysis even if the weakly visible portion is distinct for distinct pair of vertical lines. This is because each vertical line l v i is connected to at most two vertical lines (l si and lt i respectively) by dotted lines. In Fig. 13, there are four vertical line segments. Segment 1 is weakly visible by segments 2, 3, 4 in different portions specified by dotted, dashed and solid arrow headed lines. We have the following lemma. Lemma 14. The total number of grid points is O (n). Proof. Based on the above discussion, we only have to prove in one turning region the number of vertical lines in case A.1 (A.4) is constant since other cases are similar to the two cases. We denote the constant number as en1 in case A.1.

en1 

=

2r i cos

θi 2

sp 2w cos

θi 2

+

w (π − θi )

se

κmax αi (1 − sin θ2i )3

w (1 − sin

θi

) cos θ2i 2

2 θi

 3κmax αi 1 − sin

+

w (π − θi )κmax αi (1 − sin

θi 2 ) 2

w (π − θi )

2

 6πκmax . For case A.4, we compute it in the similar way and obtain the same result above. Thus the lemma follows.

2

4.2. Connecting the grid points (Problem 2) In this section, we discuss our idea on how to connect the set of grid points. The purpose of sampling grid points is to make sure that the optimal curve P passes through these points or their neighborhoods. To connect them, a naive way is to use straight line segments. But this will lead to non-smooth paths with infinitely large bending energy. To ensure the resulting path is smooth, for each grid point our idea is to first draw a set of circles tangent to the corresponding “horizontal” line (from both sides of the horizontal line) at the grid point with curvatures 0, κ , 2κ , . . . , κmax (the set of circles are said to be associated with the grid point), and then connect every pair of circles generated at different grid points by a tangent segment (i.e., a segment tangent to both circles) which is fully inside C . We now discuss how to choose κ . Notice that on one hand, κ cannot be too large, since otherwise we may not be able to guarantee that P A has a similar bending energy as P . On the other hand, κ cannot be too small, since otherwise it will be computationally costly. Thus κ has to be properly chosen to meet the two contradicting expectations. To determine κ , consider Fig. 14. Let l and l be two consecutive “vertical” lines of the grid, and p be a grid point on the vertical line l associated with two circles c 1 and c 2 with curvature κ1 and κ2 respectively and satisfying the condition |κ1 − κ2 |  κ . Let q be a grid point on l associated with two circles c 1 and c 2 in a similar way. Let  E κ1 ( E κ2 ) be the bending energy for a smooth curve using part of c 1 (c 2 ) and c 1 (c 2 ) to travel through the portion of C between l and l . Let en be the total number of vertical lines of the grid passed by P in ϕi . Our idea is to choose κ such that the following inequality holds.

L. Xu, J. Xu / Computational Geometry 47 (2014) 349–366

363

Fig. 14. Connect circles by tangent line segments.

en | E κ1 −  E κ2 | 

1 2

 E li .

(8)

The idea behind the above inequality is to ensure that the difference of bending energy between any pair of smooth curves that travel through the region ϕi with a curvature difference no more than κ is within an  -factor. To find the exact value of κ , we consider cases A.1 and A.4 since other cases are similar. In case A.1, it is easy to see that en =

(

2w cos

(1−sin

θi 2

θi 2

)sp

+

w (π −θi ) se ).

Since |κ1 − κ2 |  κ , | E κ1 −  E κ2 | satisfies



κ 2π κ 2π κ  2π κ  2π | E κ1 −  E κ2 |  2 1 − 2 + 1  − 2  κ1 κ2 κ1 κ2 Combining (9) and (8), we have κ 

 2 24π κmax w2

 4κπ .

(9)

. Thus we can choose κ =

 2 24π κmax w2

.

We obtain the same result for case A.4. The lemma below follows from Lemma 14 and the fact that κ is a constant. Lemma 15. The total number of circles in C is O ( n ). Now we discuss our idea on how to connect the circles associated with different grid points. Let c 1 and c 2 be two circles associated with different grid points. c 1 and c 2 may be (1) in the same elbow region or (2) in different elbow regions. In (2), we have to add inter-elbow edges into G to ensure that the path we computed has 0 energy in the region if P has 0 energy in the same region. To connect c 1 and c 2 , we first compute the four possible tangent segments between c 1 and c 2 . Since only those tangent segments inside C are generated, we have to determine the visibility of the endpoints of each tangent segment. This implies that we need to check the visibility for the endpoints of every tangent segment. To solve this problem more efficiently, we treat each pair of the endpoints (i.e., tangent points) as a visibility query, and build a data structure to answer those queries. Our idea is to use the optimal data structure developed in [7] in which a data structure was presented for solving the ray shooting problem inside a simple polygon P O with n vertices. In particularly, it was shown that the data structure of size n can be preprocessed in O (n) time, and for each pair (q, u ) of point q and direction u, it finds in O (log n) time the first edge of P O intersected by the ray emitting from q moving in the direction of u. In [8], a much simpler algorithm was presented to solve the visibility query problem. We can use either of them for our purpose. To make use of this data structure, we first convert the corridor into a simple polygon C P and then build the data structure for C P . Thus, for each tangent segment, we can determine whether it is an interior tangent segment of C in O (log n) time. The tangent points on each circle c i partitions c i into a linear number (in terms of the number of tangent points) of arcs. Some arcs may go outside of the corridor. Since each circle is drawn using the local feature of the elbow region, we can determine whether it is entirely inside C in O (1) time. 2

Lemma 16. The total number of interior arcs and interior tangent segments is O ( n 2 ) in the worst case. Proof. The lemma follows from Lemmas 14 and 15 and the facts that each circle has at most O (n/ ) tangent points, thus O (n/ ) tangent segments and arcs. 2 With the set of interior tangent segments and arcs, we can now construct the graph G. In this graph, the vertices are the set of interior arcs. Each arc has an edge to its two neighboring arcs in the same circle. For two arcs in difference different circles, there is an edge between them if and only if the following two conditions are met: (i) There is an interior tangent segment connecting them (i.e., the tangent segment is between a pair of endpoints of the two arcs); (ii) The two arcs and the tangent segment together form a smooth curve. Thus there are two types of edges in G, tangent segments and edges connecting neighboring arcs in the same circle. Since the number of the second types of edges is no more than that of the arcs, the total number of edges is still bounded by O (n2 / 2 ). Also note that in the above visibility query computation, we treat s and t as circles of radius 0. Thus s and t are vertices of G and are connected to other vertices by tangent segments.

364

L. Xu, J. Xu / Computational Geometry 47 (2014) 349–366

Fig. 15. (a) Step 1 in case I. (b) Construct P si . (c) Tangent line (dashed).

The weight of each edge is 0 and the weight of each vertex is the bending energy of the corresponding arc. The lemma below follows directly from Lemma 16 and the above construction of the graph. 2

2

Lemma 17. In graph G, | V | = O ( n 2 ) and | E | = O ( n 2 ) in the worst case. Furthermore, for any path in G, its corresponding path is a smooth curve in C . 4.3. The existence of a smooth curve with similar bending energy (Problem 3) The following lemma shows the existence of a smooth curve in G with similar bending energy to that of P . Lemma 18. There is a curve P s in G close to P with bending energy E s  (1 +  ) E b , where E b is the bending energy of P . Proof. We prove the lemma by constructing a path P s with similar energy to P . Given any P , we first compute the shortest path SP in C . Thus we could obtain a set of turning regions as discussed in the last section. Let C i be any turning region, and P i be the portion of P inside C i . Denote the portion of P s in C i as P si . In order to prove the lemma, we have to show (1) P si has a similar energy to P i in each C i and (2) P s has similar energy to P after we connect each P si . Now we first prove (1). From the above discussion, we know that C i is a turning region in either case (A) or case (B). Thus we have two cases to consider, (I) C i contains no short elbow, and (II) C i contains at least one short elbow. Accordingly we have two cases about P i . I: In this case, the bending energy E P i of P i is positive, since C i contains no short elbow. We have four subcases corresponding to those in case (A) given in the last section. We only prove case (A.1) since (A.2), (A.3) and (A.4) are similar. In (A.1), we construct P si in C i in the following way. (1) On each “vertical line” in C i , we choose the two closest grid points to P i (one above and one below). (2) For each chosen grid point we compute k to generate a series of circles. (3) The circles are connected by tangent lines as discussed above Section 4.2. Now we describe each step and prove the quality of P si inside C i . In step (1), we call the two closest grids point to P i on a vertical line kk g 1 and g 2 . It is easy to see that at least one of them is inside C i . We draw a circle c 12 on kk with radius κ 1 in which the diameter ends at g 1 and g 2 (see max Fig. 15(a)). We draw the similar circle on each “vertical line” in C i . Obviously those circles are entirely inside C . We sort the “vertical line” by its closeness to the starting point s and consider each vertical line as a time t, the more closer to s the earlier time. Thus the set of circles on all “vertical lines” could be considered as a sample of circle’s movement inside C . One circle on a “vertical line” corresponds to the circle’s location at that specific time. Since P is inside C , every circle constructed in this way is entirely in C . Thus if we connect those consecutive circles by a thin pipe bounded by two separate tangent lines (which indicate the movement path of the circle), the region CC R covered by the whole movement is entirely in C . Now we construct P si in C c to approximate the energy of E P i . Consider a portion P p of P i with same curvature κ p , our goal is to construct a corresponding portion P sp of P si satisfying (1) with similar shape and energy and (2) inside C . We first increase κ p to the closest j κ for some integer j. In other words, we shrink P p with radius 1/κ p to P sp with reduced radius j 1 κ . P sp has less length than P p and thus we could move it up or down to pass g 1 or g 2 depending on which one makes P sp entirely inside in CC R (see Fig. 15(b)). Since the CC R is entirely in C , P si is a smooth curve satisfying (1) passing one grid point on the vertical line and (2) P si is inside C . We repeat the above procedure for each portion of P i with same curvature and then connect those consecutive P sp s by tangent lines to form P si . P si is still entirely in C since each tangent line is entirely inside CC R (see Fig. 15(c)). From the analysis of choosing sp, se and k, P si has similar energy as P i . II: In this case, C i contains at least one short elbow. In order to construct P si , we first compute a higher level shortest path in C i (say with recursion level j), based on the optimal point pk on kk of the last short elbow in C i . Let C i be any j

j

adjusted turning region corresponding to that level shortest path, and g 1 and g 2 be the two closest grid points close to

L. Xu, J. Xu / Computational Geometry 47 (2014) 349–366 j

365

j

pk on kk . Note that we choose one of g 1 and g 2 that makes smaller αi angle change for the shortest path computed by pk . II.1: (a) If C i has no short elbow, then we construct P si in the portion of C i as in case I. By the argument in case I, the

lemma holds for the portion of P si in C i . Otherwise, (b) if either g 1 or g 2 are visible to at least one of another pair of g 1x and g 2x for some x < j in the other end of C i , then we can construct P si in C i by connecting circles at the visible grid points. This operation could be guaranteed by the addition of inter-elbow grid points. Thus, the energy of such constructed P s is within a (1 +  )-factor away from that of P . II.2: Otherwise, we reconstruct P s recursively, following the recursive computation of the shortest paths in C i level by level until either (a) or (b) in case II.1 is encountered. Now we prove (2). If a portion of C is overlapped by two C i s, there is no need to prove (2). The lemma holds by Lemma 12. Thus we only need to show how to connect P si and P sj if the corresponding C i and C j doest not overlap. It is easy to see that we could use tangent line to connect the circle generated by the last grid point on P si and the circle generated by the first grid point on P sj . The tangent line is inside C since C c is formed by the whole closest grid points to P inside C i and C i is inside C . j

j

From the above discussion, we know that the bending energy of P s in each turning region is always within a factor of (1 +  ) away from that of P . Thus the lemma follows. 2 4.4. The (1 +  )-approximation algorithm (Problem 4) In this section, we describe how to find the path P A with E a  (1 +  ) E b . Algorithm ((1 +  ) minimum bending energy path approximation algorithm). Input: A corridor C , the maximum curvature κmax and  Output: (1 +  )-approximation path P A 1. Perform an isotropic scaling on C to make w be a large enough constant. 2. Find the turning region by computing the shortest path SP and all maximum αi s for short elbows. 3. For each turning region (or adjusted turning region) C i , determine s and s and generate the gird points. Add inter-elbow gird points if necessary. 4. Determine κ and draw a set of circles for each grid point. 5. Connect the circles by tangent line segments and determine interior tangent segments and arcs. 6. Construct graph G based on interior arcs and tangent segments. 7. Find the shortest path P A between s and t in G. 8. Scale the corridor and P A back to its original size. The theorem below follows from all previous lemmas and the shortest path algorithm (e.g., Dijkstra’s algorithm) running in O ( V log V + E ) time where E is the number of edges and V is the number of vertices. Theorem 1. There exists a (1 +  )-approximation algorithm for the curvature bounded MinBEP problem with worst case running time O (( n )2 log( n )). 5. Conclusion In this paper, we propose a novel approach based on a mix of geometric and approximation algorithm techniques for the MinBEP problem. We first show how to lower bound the bending energy of an optimal curve with bounded curvature, and then use this lower bound to design a (1 +  )-approximation algorithm for this restricted version of the MinBEP problem. It is the first algorithm with a guaranteed performance ratio for the MinBEP problem. There are some other related open problems. For example, how to solve the MinBEP problem in arbitrary simple polygon or polygonal region? Since the MinBEP problem considered in this paper is motivated by interventional applications for cardiovascular surgery, an approach applicable to the 3D case is also highly desirable. References [1] P. Agarwal, T. Biedl, S. Lazard, S. Robbins, S. Suri, S. Whitesides, Curvature-constrained shortest paths in a convex polygon, SIAM J. Comput. 31 (6) (2002) 1814–1851. [2] P. Agarwal, P. Raghavan, H. Tamaki, Motion planning for a steer-constrained robot through moderate obstacles, in: Proceedings of the Twenty-seventh Annual ACM Symposium Theory on Computing, 1995. [3] J. Backer, D. Kirkpatrick, Finding curvature-constrained paths that avoid polygonal obstacles, in: Proceedings of the Twenty-third Annual ACM-SIAM Symposium on Discrete Algorithms, 2007. [4] J. Backer, D. Kirkpatrick, A polynomial-time algorithm for finding a bounded-curvature path in a narrow corridor, in: 16th Fall Workshop on Computational and Combinatorial Geometry, 2008.

366

L. Xu, J. Xu / Computational Geometry 47 (2014) 349–366

[5] S. Bereg, P. Bose, A. Dumitrescu, F. Hurtado, P. Valtr, Traversing a set of points with a minimum number of turns, in: Proceedings of the Twenty-third Annual ACM-SIAM Symposium on Computational Geometry, 2007. [6] S. Bereg, D. Kirkpatrick, Curvature-bounded traversals of narrow corridors, in: Proceedings of the Twenty-first Annual ACM-SIAM Symposium on Computational Geometry, 2005. [7] B. Chazelle, L.J. Guibas, Visibility and Intersection Problems in Plane Geometry, Discrete Comput. Geom. 4 (6) (1989) 551–581. [8] D. Chen, H. Wang, Weak visibility queries of line segments in simple polygons, in: Proceedings of the Twenty-third International Symposium on Algorithms and Computation, 2012, pp. 609–618. [9] R. Drysdale, G. Rote, A. Sturm, Approximation of an open polygonal curve with a minimum number of circular arcs and biarcs, Comput. Geom. 41 (2008) 31–47. [10] L. Euler, Methodus Inveniendi Lineas Curvas Maximi Minimive Propriatate Gaudentes, Sive Solutio Problematis Isoperimetrici Lattisimo Sensu Accepti, Bousquet, Lausannae et Genevae, E65A. O.O. ser. I, vol. 24, 1744. [11] S. Fortune, G. Wilfong, Planning constrained motion, Ann. Math. Artif. Intell. 3 (1) (1991) 21–82. [12] M. Freedman, Z. He, Z. Wang, On the Mobius energy of knots and unknots, Ann. Math. 139 (1) (1994) 1–50. [13] L. Guibas, J. Hershberger, D. Leven, M. Sharir, R. Tarjan, Linear time algorithm for visibility and shortest path problems inside simple polygons, Algorithmica 1 (2) (1987) 133–162. [14] G. Harary, A. Tal, 3D Euler spirals for 3D curve completion, in: Proceedings of the Fifteenth Annual ACM Symposium on Computational Geometry, 2010. [15] R. Heekeren, F. Faas, L. Vliet, Finding the minimum-cost path without cutting corners, Lect. Notes Comput. Sci. 4522 (2007). [16] T. Ivey, D. Singer, Knot types, homotopies and stability of closed elastic rods, Proc. Lond. Math. Soc. 79 (1999) 429–450. [17] M.K. Konings, E.B. van de Kraats, T. Alderliesten, W.J. Niessen, Analytical guidewire motion algorithm for simulation of endovascular interventions, Med. Biol. Eng. Comput. 41 (6) (2003) 689–700. [18] P. Jacobs, J. Canny, Planning smooth paths for mobile robots, in: Proceedings of the 1989 IEEE International Conference on Robotics and Automation, 1989. [19] D. Lutterkort, J. Peters, Smooth paths in a polygonal channel, in: Proceedings of the Fifteenth Annual ACM Symposium on Computational Geometry, 1999. [20] H. von der Mosel, Minimizing the elastic energy of knots, J. Asymptot. Anal. 18 (1–2) (1998) 49–65. [21] J. O’Hara, Energy of a knot, Topology 30 (1991) 241–247. [22] J. Reif, H. Wang, The complexity of the two dimensional curvature-constrained shortest-path problem, in: Proceedings of the 3rd Workshop on the Algorithmic Foundations of Robotics, 1998. [23] J.R. Sack, J. Urrutia, Handbook of Computational Geometry, North-Holland Publishing Co., 2000, pp. 854–855. [24] S. Schafer, V. Singh, K. Hoffmann, P. Noël, J. Xu, Planning image-guided endovascular interventions: Guidewire simulation using shortest path algorithms, Proc. SPIE 6509 (2007). [25] S. Schafer, V. Singh, P.B. Noël, J. Xu, A. Walczak, K.R. Hoffmann, Real time endovascular guidewire position simulation using shortest path algorithms, Int. J. Comput. Aided Radiol. Surg. (IJCARS) 4 (2009) 597–608. [26] J. Thorpe, Elementary Topic in Differential Geometry, Springer-Verlag, 1979. [27] T. van Walsum, S.A.M. Baert, W.J. Niessen, Guide wire reconstruction and visualization in 3D using monoplane fluoroscopic imaging, IEEE Trans. Med. Imaging 24 (5) (2005) 612–623. [28] H. Wang, P. Agarwal, Approximation algorithms for curvature-constrained shortest paths, in: Proceedings of the Seventh Annual ACM-SIAM Symposium on Discrete Algorithms, 1996. [29] L. Xu, Y. Tian, X. Jin, J. Chen, S. Schafer, K. Hoffmann, J. Xu, An improved endovascular guidewire position simulation algorithm, in: Proceedings of the Ninth IEEE International Symposium on Biomedical Imaging, 2012, pp. 1196–1199.