Corner Detection and Arc Approximation of Planar Data Set
CCorner Detection and Arc Approximation of Planar Data Set
Maurizio Scarparo
Abstract
The aim of this paper is to present a new method of approximation of planar data set usingonly arcs or segments. The first problem we are trying to solve is the following: the CNCmachines can work only with simple curves (arcs or segments, indeed). So, if we have acontour of an object and we want to cut, for example, glass or metal following that profile,we need a continuous curve composed by linear or circular paths only. Moreover, we wantto minimize the number of these paths. The second problem is the following: the contour ofan object is detected by a specific laser, which collects a discrete data set. The laser detects agood approximation of a point if this point is not a corner of the object. So in many cases ourdata set will not contain the real corners. Our task is to present a method on how to find theseparticular points. A third purpose is to study the regularity of the approximating curve, whichcan be G or C continuous. The entire method is developed using single arcs, segments andin certain cases biarcs, in order to ensure the smoothness of the final path. Keywords : arcs, biarcs, least-square fitting, point approximation, corner detection, digitalcurves
Contacts : [email protected]
The problem of approximation of planar dataset has been already studied by several au-thors, such as Les A. Piegl, Wayne Tiller,Hyungjun Park, D. Meek, D. Walton (see page17). In their articles these authors explainthe biarc model and how to apply it in CNCmachining. One of the most powerful meth-ods, introduced by Piegl and Tiller, consistsin finding the least square B-spline and ap-proximating this one with biarcs (see [7], [8],[9]). The problem of corner detection has alsobeen studied by authors such as C. Harris, M.Stephens, H. Freeman, L. Davis. However, wedecided to analyze these problems in a dif-ferent way. The first task was to elaborate a new method of approximating a planar dataset without using B-splines. Indeed, we foundout that the B-spline method is very powerfulwhen the data set has no corner points, butit is quite difficult to get a good approximationwhen the data set is not smooth. Moreover, wefound out that the approximation heavily de-pends on the data set and on the choice of thetangent vectors, so we wanted to cancel thedependence on the vectors. We also discoveredthat single arcs are better than biarcs if wewant to minimize the number of paths. Sec-ondly, we wanted to develop a straightforwardmethod for corner detection in a planar dataset captured by a laser. As already explainedin the abstract, this data set may not include1 a r X i v : . [ m a t h . NA ] N ov he exact corner points, because of the laser er-rors. This is the main difference between ouralgorithm and the methods explained in [1],[2], [6], [11], [12]: all methods developedthere need the corner points to be included inthe data set. We noticed that this problem hasbeen carefully studied for digital images usingpixel properties, but as far as we know thereare not many articles concerning this topic out-side image analysis. The first time we try to approximate a planardata set with only biarcs we can find out an in-convenient: if the distance between two con-secutive points is small and the tangent vectorsare quite different, the biarc between the twopoints will be serpentine , that means uselessand insignificant for a good approximation. Sothe first step is to forget the biarc model fora while and to focus on the least square con-strained arcs . Let P , P , . . . , P N +1 be N + 2 distinct planar points. We want to find oneparticular arc with initial point P , final point P N +1 , approximating P , . . . , P N and mini-mizing an objective error function. Let R , C ,be respectively the radius and the center of ourarc. The function we want to minimize is thefollowing: f ( t ) = N (cid:88) k =1 (cid:12)(cid:12)(cid:12)(cid:12) R − (cid:107) P k − C ( t ) (cid:107) (cid:12)(cid:12)(cid:12)(cid:12) . (1)The variable t represents the distance betweenthe centre C and the segment P P N +1 . In par-ticular, let d = (cid:107) P N +1 − P (cid:107) be the distancebetween P and P N +1 , let u be the normalizedvector u = P N +1 − P d (2) Notation: P i = [ x i , y i ] and let u (cid:48) be the π counterclockwise rotationof u . Then, the center C can be representedby the following: C ( t ) = P + d u + t u (cid:48) . (3)Let us discuss about the objective function (1).This is a non negative continous function on R .One can think that f depends on two unknown R and C ; however, since the initial point of ourarc is P and the final one is P N +1 , the radiusdepends on the center and R ( t ) = (cid:107) P − C ( t ) (cid:107) = (cid:107) P N +1 − C ( t ) (cid:107) . (4)So, according to (4), the function f in (1) be-comes f ( t ) = N (cid:107) P − C ( t ) (cid:107) + N (cid:88) k =1 (cid:107) P k − C ( t ) (cid:107) − (cid:107) P − C ( t ) (cid:107) N (cid:88) k =1 (cid:107) P k − C ( t ) (cid:107) . (5)According to Weierstrass theorem, f admitsminimum on each compact set [ a, b ] ⊂ R .Moreover, one can easily see that the firstderivative of f is always continuous on R if thecenter C ( t ) and the points P k are distinct forevery k = 0 , , . . . , N + 1 . We will prove that f admits finite limits at ±∞ and these limits arethe same. Theorem 2.1. f has finite limit at ±∞ and lim t →±∞ f ( t ) = N (cid:88) i =1 (cid:20) u (cid:48) · ( P − P i ) (cid:21) , (6) where · means the standard scalar product in R . roof. We have that (cid:107) P − C ( t ) (cid:107) = (cid:13)(cid:13)(cid:13)(cid:13) d u + t u (cid:48) (cid:13)(cid:13)(cid:13)(cid:13) = t + d , (cid:107) P i − C ( t ) (cid:107) = t + 2 t (cid:20) u (cid:48) · (cid:18) P − P i + d u (cid:19)(cid:21) + (cid:13)(cid:13)(cid:13)(cid:13) P − P i + d u (cid:13)(cid:13)(cid:13)(cid:13) . (7)Let D i = P − P i + d u . We notice that f is thesum of N functions f i ( t ) such that f i ( t ) = ( (cid:107) P − C ( t ) (cid:107) − (cid:107) P i − C ( t ) (cid:107) ) == ( g ( t ) − g i ( t )) , where g ( t ) = (cid:114) t + d ,g i ( t ) = (cid:113) t + 2 t (cid:0) u (cid:48) · D i (cid:1) + (cid:107) D i (cid:107) . (8)Multiplying f i ( t ) by ( g ( t )+ g i ( t )) ( g ( t )+ g i ( t )) one obtains f i ( t ) = ( g ( t ) − g i ( t )) ( g ( t ) + g i ( t )) == (cid:18) d − t (cid:0) u (cid:48) · D i (cid:1) − (cid:107) D i (cid:107) (cid:19) (cid:18)(cid:113) t + d + (cid:113) t + 2 t (cid:0) u (cid:48) · D i (cid:1) + (cid:107) D i (cid:107) (cid:19) (9)By equation (9) one can easily see that lim t →±∞ f i ( t ) = (cid:0) u (cid:48) · D i (cid:1) = (cid:0) u (cid:48) · ( P − P i ) (cid:1) , (10)and (6) is proved using limit properties on thesum. Corollary 2.1.
The limit of f at ±∞ is zeroif and only if all the points P , . . . , P N arecollinear on the segment P P N +1 . Theorem 2.1 is quite important because wecan say that, if L is the limit of f at ±∞ , ∀ ε > ∃ t ∈ R : | t | > | t | = ⇒ | f ( t ) − L | < ε. So the set B (0 , | t | ) = { t ∈ R : | t | ≤ | t |} is com-pact, f is continuous on B (0 , | t | ) and admitsminimum there. In order to find a minimumof f one can apply one of the optimizationtechniques, such as Newton or Quasi-Newtonmethods, BFGS methods and others. In liter-ature, one can find many methods of build-ing a least square circle given some planarpoints. One of the most efficient methods isTaubin’s method, explained in detail in [15].The method explained int this paper has onemain difference with Taubin’s: the least squarearc is constrained, that is it passes through twogiven points. Taubin’s arc is better than ours,because it minimizes the error on all the pos-sible centers and radii, but in CNC approxima-tion it’s very difficult to find two consecutiveTaubin’s arcs connecting each other in a con-tinuous way.In the following page, there is a simple ex-ample of the function f and the related leastsquare arc. In this section we want to explain how to ap-ply least square constrained arcs in the approx-imation of a planar data set without corners(we suppose that the profile of our object islocally smooth). Suppose that P , P , . . . , P N are N consecutive points approximating onesection of the object contour.The first step of the method consists in try-ing to approximate the first four points witha constrained arc (it starts from four becausewe already know that there always exists onearc passing through three points, eventuallywith infinite radius). We decide an error limitand we compute the distance of each point P k f (on the right).from the arc. If there exists one point whosedistance is greater than our limit, we decreasethe number of approximated points and, inthis case, we build the arc passing through thefirst three points. Conversely, if all the pointsare within the error limit, we increase thenumber of approximated points and we try tobuild a longer arc. The method repeats thesesteps for each point P k , k = 1 , , . . . , N − .At the end we obtain a N × matrix M :in the k th place there will be the number ofpoints that can be approximated with a singleleast square arc starting from P k — includedthe first and the last points —. We will set M ( N −
1) = M ( N ) = 0 . Here follows thealgorithm (which needs at least four points inorder to be applied):1. Set i = 1 , k = 3 , f lag = 1 and choose ε > ;2. While i ≤ N − do
3. Set l = k − ;4. Find the center C and the radius R of the least square arc passingthrough P i , P i + k ; 5. For j = 1 , . . . , l d = | R − (cid:107) P i + j − C (cid:107)| ;7. If d > ε set f lag = 0 ; break ;8. If f lag = 0
9. build the least square arcpassing through P i , P i + k − M ( i ) = k , f lag = 1 , i = i + 1 ;11. If i = N − set k = 2 Else set k = 3 ; Else k = k + 1 ;13. If i + k > N
14. build the least squarearc passing through P i , P i + k − M ( i ) = k , i = i + 1 ;16. If i = N − set k = 2 Else set k = 3 ;17. M ( N −
1) = 0 , M ( N ) = 0 .4ne can notice that M ( k ) ≥ for all k =1 , . . . , N − , so that the worst case occurswhen all the elements of M are equal tothree. Next step is to decide how to choosethe longest arcs. We will explain the choice of longest arcs withtwo simple examples. Fix k ≥ and supposethat M ( k ) = n , for a certain n ≥ ∈ N .Suppose also that M ( i ) < M ( k ) for all i = k + 1 , k + 2 , . . . , k + n − . This is the sim-plest case: the arc starting from P k is approx-imating an higher number of points than the n − consecutive arcs, so it is good for an op-timized approximation. Now suppose that anindex j exists such that ≤ j ≤ n − and M ( k + j ) > M ( k ) . In this situation we can’tsay that the arc starting from P k is the best forour task. Let S = (cid:22) M ( k )2 (cid:23) and let I k be the set I k = { j ∈ N | k + 1 ≤ j ≤ k + S } . Then, we decide to apply the following test: if M ( j ) ≤ M ( k ) for all j ∈ I k , the arc startingfrom P k is good for our approximation and wepass from P k to P k + M ( k ) − with a single arc.We repeat the process on P k + M ( k ) − (the newstarting point). Conversely, if an index j ∈ I k exists such that M ( j ) > M ( k ) we pass directlyfrom the point P k to P j without any segmentsor arcs — for the moment — and we applyagain the test on the new starting point P j .This kind of selection is justified by the follow-ing fact: suppose there exists an index j / ∈ I k such that M ( j ) > M ( k ) and M ( k ) is a highnumber ( M ( k ) ≥ for example ); then, thegap between P k and P j can be so wide that we are not able to obtain a good approxima-tion for all those points among them, so thetest is applied only on the set I k .Here follows the algorithm ( N is, as usual, thenumber of points we are approximating):1. Set i = 1 , G = M ;2. While i < N do S = (cid:22) M ( i )2 (cid:23) , j = i + 1 , f lag = 1 , e = i + M ( i ) − ;4. While j ≤ i + S and f lag = 1 do If M ( j ) ≤ M ( i ) set j = j + 1 Else
6. Set f lag = 0 ;7. for k = i, . . . , j − set G ( k ) = 0 ;8. Set i = j ; e = i + M ( i ) − ; break ;9. If f lag = 1 For k = i + 1 , . . . , e − G ( k ) = 0 ;11. Set i = e ;In this algorithm we introduce a new matrix G .At the beginning, G is equal to M , but duringthe steps the element G ( i ) may become null.This happens in two cases, if P k is a startingpoint: • M ( k ) ≥ M ( j ) for all j ∈ I k . Then we set G ( j ) = 0 for all j such that k + 1 ≤ j ≤ k + M ( k ) − ;5 there exists some index j ∈ I k such that M ( k ) < M ( j ) . Then we set G ( i ) = 0 forall i ∈ I k satisfying k + 1 ≤ i ≤ j − andfor i = k .At the end of the process we obtain a N × matrix G with elements G ( k ) satisfying: (cid:5) G ( k ) > , if there exists a single arc pass-ing through the points P k , P k + G ( k ) − andapproximating all the other points amongthem, i.e. if the point P k is a starting pointfor one arc; (cid:5) G ( k ) = 0 , if the point P k is not a startingpoint, so G ( N −
1) = G ( N ) = 0 .Suppose now G ( k ) > and G ( k + G ( k ) − > : then there exist two consecutive arcs (con-tinous arcs) with common point P k + G ( k ) − .This is the best case: all points between P k and P k + G ( k ) − G ( k + G ( k ) − − are approxi-mated within the tolerance. But what happenswhen two consecutive arcs don’t have a com-mon point? This is the next step: how to fillthe gaps between the built arcs. Before proceeding to the next step, we want todiscuss about biarcs and their applications.
Definition 5.1.
We say that two circulararcs a , a form a constrained biarc passingthrough given points P s , P e with unit tangentvectors T s , T e , if they have a common point P satisfying the following properties: ∗ a starts from P s , ends in P and T s is tangent to a in P , with orientationcorresponding to a parametrization of a from P s to P ; ∗ a starts from P , ends in P e and T e is tangent to a in P e , with orientationcorresponding to a parametrization of a from P to P e ; ∗ a and a have proportional tangent vec-tors in P (i.e. the biarc is G continousin P ).In this paper we will use the standard notationused by Piegl and Tiller in [7], [8]. The controlpoints of the biarc are determined by P = P s + α T s , P = P e − β T e , with both α and β positive. Using this conven-tion, the sweep angle of each arc of the biarcwill always be in [0 , π ) . Now, let V = P e − P s .We will prove that there always exists a biarcpassing through P s , P e when T s · T e (cid:54) = 1 (11)and V · V (cid:54) = 2( V · T s )( V · T e ) T s · T e − . (12) Theorem 5.1.
Given two distinct planar points P s , P e with associated unit vectors T s , T e satis-fying (11) and (12), there always exists a biarcpassing through P s , P e with tangents T s , T e .Proof. Let P = P s + α T s , with α > . Thenit’s sufficient to prove that there exists a point P such that P = P e − β T e , β > , and (cid:107) P − P (cid:107) = α + β . This is equivalent to thecondition ( P − P ) · ( P − P ) = ( α + β ) . Substituting the values of P , P , we find thecondition ( V − α T s − β T e ) · ( V − α T s − β T e ) = ( α + β ) . Expliciting the products, we find V · V − α V · T s − β V · T e +2 αβ ( T s · T e −
1) = 0 . So, expliciting β : β = 2 α V · T s − V · V α ( T s · T e − − V · T e . (13)6ow, assume that exists α > such that β > .Then, one can define the junction point P as P = βα + β P + αα + β P . Since (cid:107) P − P (cid:107) = α and (cid:107) P − P (cid:107) = β , wecan build two arcs c , c satisfying the condi-tions of the Definition 5.1.The proof of the theorem uses the following Lemma 5.1.
Given two distinct planar points P s , P e with associated unit vectors T s , T e sat-isfying [11] and [12], there always exists a pos-itive value of α such that β > , where β is de-fined in equation [13].Proof. Let α = V · V V · T s , α = V · T e T s · T e − ,K = V · T s T s · T e − . Then, by hypothesis, α (cid:54) = α . From equation(13) we can write β = K α − α α − α . There are three cases to be examined. • K > 0. We have V · T s < , so α < and β is positive for α > max { α , } . Noticethat if α = α one obtains β ≡ K > regardless α . In this case one arc of thebiarc is degenerate (a segment), but thetheorem doesn’t concern this situation. • K < 0. We have V · T s > , so α > .There are three subcases: – α ≤ < α . We choose α ∈ (0 , α ) ; – < α < α . We choose α ∈ ( α , α ) ; – ≤ α < α . We choose α ∈ ( α , α ) . • K = 0. We have V · T s = 0 , so β = − V · V T s · T e − α − α ) . We choose α > max { α , } .The condition of the lemma V · V (cid:54) = 2( V · T s )( V · T e ) T s · T e − (14)can be read in a geometrical way. Let θ , θ , θ be respectively the angles between ( V , T s ) , ( V , T e ) , ( T s , T e ) . Then, from (14), we find (cid:107) V (cid:107) (cid:54) = 2 (cid:107) V (cid:107) cos θ cos θ cos θ − , or, simplifying, cos θ (cid:54) = 2 cos θ cos θ + 1 . (15)It’s easy to see that a necessary condition tohave cos θ = 2 cos θ cos θ + 1 is that − ≤ cos θ cos θ ≤ . The condition (15) includesmany cases, but in all of these situations onecan solve the connection problem using twoconsecutive biarcs. Clearly, things get simplerwhen α = β . In this situation, there are onlytwo cases in which we can’t build a single arc: T s · T e = 1 and V · ( T s + T e ) = 0 . In bothcases we have to build two consecutive biarcs,i.e. four arcs, matching with G continuity. Ingeneral, one can define the biarc ratio r = αβ instead of choosing a value of α to have β > .The ratio r is very useful when we want tobuild a constrained biarc minimizing the dis-tance from given planar points. An efficient al-gorithm has been studied by Hyungjun Park in[5]: it’s called optimal single biarc fitting , andits approach is to search an optimized valueof r decreasing the width of the range of r ineach step. We will use this approach in thenext steps.7igure 2: Example of constrained biarc passing through P s = (0 , , P e = (5 , with tangentvectors T s = (1 , , T e = (0 , − (on the left) and construction of a biarc (on the right). Let us consider the biarc in Figure 2 (on theright). Let ab = B − A (cid:107) B − A (cid:107) , bc = C − B (cid:107) C − B (cid:107) . Thenwe have D = ( ab , bc ) = cos 2 θ, where θ is the half sweep angle of the biarc,and we suppose that θ ∈ (cid:2) , π (cid:1) . By bisectionformulas one can find cos θ = (cid:114) D , sin θ = (cid:114) − D , tan θ = (cid:114) − D D . (16)So the explicit formula for θ becomes θ = arctan (cid:114) − D D .
Since R = h cos θ and (cid:107) B − A (cid:107) = h sin θ wecan find the radius and the centre of the firstarc: R = (cid:114) D − D (cid:107) B − A (cid:107) , O = B + R cos θ ab − bc (cid:107) ab − bc (cid:107) . (17) Suppose P = A , P = B , P = C and { w , w , w } = { , cos θ, } . Then we easily ob-tain the parametrization of the first arc: G ( t ) = (cid:80) i =0 w i P i B i ( t ) (cid:80) i =0 w i B i ( t ) , t ∈ [0 , , (18)where B ni ( t ) , i = 0 , . . . , n , are the ( n + 1) Bern-stein polynomials defined by B ni ( t ) = (cid:18) ni (cid:19) t i (1 − t ) ( n − i ) , t ∈ [0 , . Now remember that a biarc has five controlpoints P , . . . , P , where the first and the lastones are the initial and final point of the biarc, P is the junction point. Let us consider thecharacteristic functions χ [0 , , χ (1 , . Then wecan express the biarc in Bézier form: B ( t ) = χ [0 , (2 t ) (cid:80) i =0 w i P i B i (2 t ) (cid:80) i =0 w i B i (2 t ) ++ χ (1 , (2 t ) (cid:80) i =2 w i P i B i − (2 t − (cid:80) i =2 w i B i − (2 t − , t ∈ [0 , , (19)8here the weight vector is { w i } i =0 = { , cos θ , , cos θ , } and θ , θ are the half sweep angles of the twoarcs.In general, we can write a PCC continous curveas a finite sum of arcs using Bézier form. Let T i f ( t ) = f ( t − i ) be the translation of f . Thena PCC curve C composed by n + 1 arcs has theform C ( t ) = n (cid:80) i =0 T i χ [0 , ( t ) T i G i ( t ) , if t ∈ [0 , n ) G n (1) , if t = n , (20)where G i ( t ) = (cid:80) j =0 w ij P ij B j ( t ) (cid:80) j =0 w ij B j ( t ) , t ∈ [0 , , (21)is a classical Bézier rational curve, w ij are de-termined in order to obtain arcs (cosine of thehalf sweep angles) and P ij are control pointsof the arc i (i.e. G i ( t ) is an arc with startingpoint P i and ending point P i ). Moreover, toachieve continuity, the control points must sat-isfy P i, ≡ P i − , , i = 1 , , . . . , n . Recall that in the last step we found the ma-trix G : G ( k ) > if P k is a starting pointfor some arc, G ( k ) = 0 otherwise. Now wewant to build the longest arcs using G and tounderstand the regularity of our final curve.This is quite simple to do when G ( k ) > and G ( k + G ( k ) − > . As we explained be-fore, this is the case of two consecutive arcsmatching with C continuity. So we build thearc starting from P k , ending in P k + G ( k ) − andminimizing the error function f mentioned insection 2 on page 2. Then, we can progress to the next point with the same method, be-cause we already know that the arc startingfrom P k + G ( k ) − exists. Generally, it’s quite dif-ficult that two consecutive longest arcs con-nect each others with C continuity. So, theidea is to proceed in a different way accordingto the number of not approximated points be-tween two consecutive arcs. We will explainbetter what we mean with two examples. Sup-pose that one of the longest arcs has its end-ing point in P k and that the consecutive arcstarts from P k +1 . In this case we can builda segment between the two points P k , P k +1 in order to have C continuity. Suppose nowthat an ending point is P k but the next startingpoint is P k +2 . In this case, to achieve continu-ity, we build the arc passing through P k , P k +1 and P k +2 , which always exists. We find theworst case when the number of not approxi-mated points is greater than . We can solvethis problem in two different ways: • using constrained biarcs and optimal sin-gle biarc fitting; • reapplying the algorithms to this smallersection in order to find a new G and newapproximating arcs.Generally, the first approach is not ideal for anoptimized curve. In fact, as we underlined atthe beginning, the biarc approximation failssometimes when the section to be approxi-mated is quite small. However, using biarcsensures G continuity, so the first method cre-ates more regular curves. Here follows the al-gorithm. The matrix G is supposed to be N × ,with the last two elements equal to zero.1. Set i = 1 , k = 0 ;2. While i ≤ N do If G ( i ) = 0 do
4, 54. i = i + 1 ; k = k + 1 ;9. If i > N i = i − ; k = k − ;7. If k = 1 build the segment P i − P i ; k = 0 ;8. If k = 2 build the circle passingthrough P i − , P i − , P i k = 0 ;9. If k ≥ reapply the algorithmson the set { P i − k +1 , . . . , P i } or try to approximatewith biarcs; k = 0 ;10. i = i + 1 ;11. Else do
12, 13, 14, 1512. If k = 1 build the segment P i − P i ; k = 0 ;13. If k = 2 build the circle passingthrough P i − , P i − , P i k = 0 ;14. If k ≥ reapply the algorithms onthe set { P i − k +1 , . . . , P i } or try to approximatewith biarcs; k = 0 ; 15. i = i + G ( i ) − ; Before studying the regularity of the finalcurve, we need a brief summary. We startfrom a set of planar data points and we tryto build all the possible arcs approximatingany number of points (construction of matrix M ). Then, we make a selection: we take inconsideration only the longest arcs (construc-tion of matrix G ). Finally, if there are somepoints which haven’t been approximated yet,we reapply all the previous steps in the miss-ing sections, creating a continuous curve com-posed only by arcs or segments.Clearly, it’s quite difficult that all the section ofthe final curve join each other with C continu-ity, because the approximation uses both arcsand segments. So the next problem we wantto discuss is about how to smooth our junctionpoints (j.p.). To do this we first need to recog-nize good j.p. from bad j.p. Let us fix a positivetolerance ε and let P k be a junction point forsome k . Suppose that P k is the j.p. of twoconsecutive arcs A k − , A k and let T e,k − , T s,k be respectively the final unit tangent vector of A k − and the initial unit tangent vector of A k .We say that P k is a bad j.p. if T s,k · T e,k − = cos ϑ < ε, (22)where ϑ ∈ [0 , π ] is the angle between the twotangent vectors. For example, if ε = 0 . ,each angle greater than 6 degrees will be a bad angle. Each angle that doesn’t satisfy (22) iscalled a good angle. Our task is to smooth out bad angles: we can only accept piecewise cir-cular curves (PCC) whose junction points are good . This is a practical motivation: if the finalobject, whose profile is defined by our curve,10igure 3: Example of good (on the left) and bad (on the right) junction points (green points).On the left the curve will have C continuity. On the right, the red and blue curve is an exampleof smoothing biarc. Notice that the biarc is within the tolerance and that the final curve will be G continuous in the junction point.Figure 4: Example of smoothing arc in a neighbourhood of a bad point.11as good junction points, it’s easier to smoothout the entire contour, because the angles ingood j.p. will be small enough. So if P k is a good j.p., we don’t apply any changes to thePCC (see Figure 3 on the left). Otherwise, inorder to increase the regularity of the curve inthat particular point, we use biarcs. Indeed, aswe have already explained, biarcs are G con-tinuous, so they can regularize quite easily a C continuous curve.Suppose P k is a j.p. of two consecutive arcs A k − , A k and let δ be a positive fixed toler-ance. Then we can find two points Q k − , Q k respectively on the arcs A k − , A k satisfying (cid:107) Q k − − P k (cid:107) = (cid:107) Q k − P k (cid:107) = δ , i.e. Q k − , Q k are on the circumference with center P k and radius δ . The idea is to employ thesetwo points as initial and final points of a biarc,starting from Q k − and ending in Q k (see Fig-ure 3 on the right). Moreover, we can easilycompute the tangent vectors in Q k − and Q k because we know the centers and the radius of A k − , A k (a tangent vector of a circle is well-known). The biarc we build with this methodregularize the PCC locally, in a neighbourhoodof the j.p. P k .It’s highly recommended to choose δ not toohigh nor too small. In the first case, we canobtain a biarc that doesn’t satisfy the tolerancechosen at the beginning of the procedure —remember that we are trying to build a PCCwhose distance from the given points is lessthan a positive number we decide —. In thesecond case we smooth out the curve with sucha small radius that the approximating biarcwould be very similar to our j.p., i.e. useless.A second way to accomplish this task is to usesimple arcs: given two arcs (possibly degener-ate in segments) with a common j.p., it’s al-ways possible to build a simple arc approxi-mating the j.p. within the tolerance and tan-gent to the other ones. If R , R are the radiiof the two consecutive arcs, than we can in-crease or decrease them to obtain a new circle intersection whose distance from the two arcsis the same. This point will be the center of thetangent arc (see Figure 4). C continuity andreducing the number of arcs In the previous sections we discussed abouthow to find the longest arcs, how to fill thegaps among them and how to smooth the fi-nal curve. Now we want to describe a differ-ent method for approximating the data set, inorder to build a C continuous circular curvewithout any gap. This method tries to finda curve composed by a low number of arcs(not necessarily the minimum). This problemcomes from a technical motivation: we wantthe CNC machine to perform a reduced num-ber of movements. Clearly, minimizing thenumber of arcs means maximizing the num-ber of approximated points for some arc. Inthis method we will use again the matrix M .The idea is quite simple but efficient. To fixideas, suppose that M (1) = 5 . Then, we canbuild a least square arc starting from P , end-ing in P approximating P , P , P . Supposealso that M (3) = 3 , M (4) = 7 , M (5) = 4 . Let A be the arc starting from P , ending in P and let A be the arc starting from P , endingin P . Then A , A will be C continuous in P and will approximate eight points in total.Now let A be the arc starting from P , endingin P and let A be the arc starting from P ,ending in P . Then A , A will be C con-tinuous in P and will approximate ten pointsin total. So A , A will be better than A , A ,because they can approximate a higher num-ber of points. Starting from the example wecan develop the general case. Let S = M ( i ) , i = 1 , . . . , N − , where N is the number ofpoints to be approximated. Then we can finda least square arc starting from P i , ending in P i + S − , and a least square arc starting from12 i + S − , ending in P i + S + M ( i + S − − . Thesetwo arcs will approximate S + M ( i + S − − points in total. Now let j be an index suchthat i + 2 ≤ j ≤ i + M ( i ) − . We can buildwith the same technique two consecutive con-tinuous arcs, the first one from P i to P j andthe second one from P j to P j + M ( j ) − , approx-imating j − i +1+ M ( j ) − j − i + M ( j ) pointsin total. So, if j − i + M ( j ) > S + M ( i + S − − ,these two arcs will be better than the first ones.Also here we will work with a new matrix G :at the beginning we set G = M ; during thenext steps, we will set G ( j ) = 0 if G ( j ) is nota starting point for some arcs (as we have al-ready done in the previous sections).1. Set i = 1 , G = M ;2. While i ≤ N − do
3. Set S = M ( i ) , k = i , j = i + 2 ;4. Set Z = S + M ( i + S − − ;5. While j ≤ i + M ( i ) − do If j − i + M ( j ) > Zk = j ; Z = j − i + M ( j ) ;7. j = j + 1 ;8. For j = i + 1 , . . . , k − G ( j ) = 0 ;9. For j = k + 1 , . . . , k + W ( k ) − G ( j ) = 0 ;10. If i (cid:54) = kG ( i ) = k − i + 1 ;11. i = k + M ( k ) − ;The matrix G determined in the algorithm isused then in the real approximation, and we can apply again the algorithm studied previ-ously. As one can immediately notice, thismethod is better than the longest-arc-searchalgorithm if we take in consideration the reg-ularity. Indeed, the final curve will be C con-tinuous. The worst case occurs when the finalcircular spline doesn’t end precisely in the lastpoint P N but in P N − . In this situation wehave to build a segment from P N − to P N inorder to achieve C continuity.
10 Corner detection
Up to now we have supposed that all the sec-tions we want to approximate don’t containany corner points. In this paragraph we aregoing to explain how to detect corner pointsin the data set. This is not an easy problem,even when the points to be approximated arenot affected by errors. We have to distinguishtwo different cases: the first one occurs whenthe data set actually contains corner points, i.e.there are some points that are also corners; thesecond one occurs when the corners are not inthe data set, i.e. all the points are not cor-ners. It’s easy to understand that the secondcase is more difficult and challenging than thefirst one. This case is caused by detection er-rors: several laser detectors can’t give the ex-act position of corners. In presence of a corner,they create a sequence of points really close toit, but each of those points can’t be consideredas a good approximation of the corner. Thefirst case can be solved with one of the algo-rithms explained in [1], [2], [6], [11], [12].So, let us focus on the second case. The triv-ial idea is the following: suppose our data setis composed by N distinct points P i ; for eachpoint we define two unit vectors T ,i = P i − − P i (cid:107) P i − − P i (cid:107) , T ,i = P i +1 − P i (cid:107) P i +1 − P i (cid:107) , α i between the two vectors, using the formula cos α i = T ,i · T ,i . Then, we say that P i isa corner if cos α i ≥ − ε , where ε is a positivevalue in [0 , . With few examples one caneasily see that this condition is not sufficientto define corners, especially when the data setfalls in the second case (when points are notcorners). So we need another condition. Thesecond requirement comes from a geometricalinterpretation of curvature. We know that thecurvature of a smooth curve at each point isthe reciprocal of the radius of its osculating cir-cle. If the curvature at the point x is a smallnumber, then the curve will be similar to a linein x .Instead, if the curvature at x is high, then thecurve will be similar to a small circle in x , i.e.the curve will quickly change direction . Wewill use the same principle to determine cor-ner points.Assume now that all the points of the dataset have been captured by laser detector withconstant step. Then, the distance from twoconsecutive points will be approximately thesame. Let us choose two values ε ∈ [0 , , R max > and consider three consecutivepoints P i − , P i , P i +1 satisfying cos α i ≥ − ε .Let R be the radius of the circle passingthrough those three points. We say that P i isan anchor point if R < R max .Notice this: being an anchor point doesn’t im-ply being a corner point. If P i is an anchorpoint, then it’s quite probable that the actualcorner point is in a neighbourhood of P i . Inorder to localize the corner we need two an-chor points. Then, we compute α i +1 , α i − with the same technique described above. If cos α i +1 ≤ cos α i − , the two anchor points willbe P i and P i +1 . Otherwise, they will be P i and P i − . Here one can understand the assump-tion that the distance from any two consecu-tive points is approximately the same: if thepoints are detected with adaptive steps tech- niques, it’s hard to find a right value of R max .After finding all the anchor points (two ata time), we use them to find corner points.Clearly, the corners need to be placed betweenthese two points. The idea is to use anothertime the least square arcs. Suppose that P i and P i +1 are two consecutive anchor points.Let us consider the sets Λ i,k = { P i − j : j = 0 , . . . , k } , Λ i +1 ,h = { P i + j : j = 1 , . . . , h } , where k, h ∈ N and k ≥ , h ≥ . We wantto increase k, h in order to build two arcs A i,k , A i +1 ,h approximating respectively Λ i,k , Λ i +1 ,h within a given positive tolerance δ . Notice this:the arcs we are building are not constrainedarcs, i.e. they don’t pass through any twogiven points. This time the best solution comesfrom Taubin’s algorithm. For k = 2 , h = 3 ,the solution is trivial because the circle pass-ing through three distinct points always exists(degenerate in a line, eventually). So we try toincrease k and h (up to a maximum limit M )and we build Taubin’s circles. If the distancebetween Λ i,k and A i,k is less than δ , k need tobe increased.Otherwise, we decrease k and, in this case, webuild the three point circle. We do the samecontrol for A i +1 ,h . At the end of the process,we will find two right values for k and h andwe will build the two corresponding arcs A i,k , A i +1 ,h . These two arcs must intersect eachother in two points. We select the nearest pointto P i . This point will be the approximation ofthe corner point.Pay attention to the fact that if the circles don’tintersect each other, then there is somethingwrong in the data set or the algorithm detectedtwo false anchor points. Another considera-tion: the corner points will be located approxi-mately between the two anchor points becausewe used circles to find them: Taubin’s methodallows to create the best circle approximating14igure 5: Example of undefined corner (on the left) and approximation of it (red point on theright). Here a tolerance of δ = 1 mm was used. The green circles represent the anchor pointsdetected by our algorithm using R max = 20 mm and ε = 0 . . The blue lines are least squarecontrained arcs.a sequence of points, so it’s reasonable to as-sume that A i , A i +1 are locally good approxi-mations of the object contour.1. Find two consecutive anchor points P i , P i +1 ;2. Set k = 3 , h = 4 , f lag k = 1 , f lag h = 1 ,choose an integer M ≥ and a tolerance δ > ;3. While k, h ≤ M do
4. Find Taubin’s circles A i,k , A i +1 ,h approximating Λ i,k , Λ i +1 ,h ;5. If f lag k = 1 If dist ( A i,k , Λ i,k ) < δk = k + 1 ; else f lag k = 0 ;7. If f lag h = 1 If dist ( A i +1 ,h , Λ i +1 ,h ) < δh = h + 1 ; else f lag h = 0 ;9. If f lag k = 0 and f lag h = 0 return k − , h − ; break ;10. Return k − , h − .
11 Testing the methods
In this last section we want to give some exam-ples of application of the two methods. First,we want to find a good approximation for thecontour of the object represented in Figure 6on the next page. We give a first rough approx-imation with 134 points captured by a laserdetector (see Figure 6 again).As one can see, this is a bad situation: thepoints are not uniform, i.e. they were deter-mined with an adaptive method. We selectthis example because we want to show that thetwo methods also work in difficult cases. First,we try to improve the data set: since there arewide gaps between some points, we add six15igure 6: Metal 2D object to be approximated (on the left) and points captured by laser detector(on the right). Figure 7: Improvement of the data set.16igure 8: Approximation of the data set using longest arcs approach.points (following a procedure that uses arc in-tersections) and we obtain the data set repre-sented on Figure 7. Among the new points,the most important one is the green one, be-cause the gap in that section of the contouris too wide to hope for a correct approxima-tion. Secondly, we apply the detection algo-rithms and we find all five corner points. Weuse R max = 20 mm , ε = 0 . (for detection ofanchor points).Table 1: Results with the first method.Tolerance Arcs Segments Bad Points1.5 mm 14 1 21 mm 17 3 60.5 mm 22 3 6Notice this: adding the green points is crucial.Indeed, in this case the research of all the cor-ner points fails, because the initial data set is Table 2: Results with the second method.Tolerance Arcs Segments Bad Points1.5 mm 18 0 31 mm 21 1 20.5 mm 28 0 4ill-posed in a neighbourhood of one of the cor-ners, i.e. there are not enough points to definea good approximation. Using the two meth-ods, we find the results in Table Here the tolerance is expressed in millimetres. Re-member that 1 mm ≈ . , the final re-sult is not reliable, because in the best casewe would obtain a curve composed by arcs ap-proximating a low number of points (three orfour), so we would also have a high number ofbad points.Secondly, we want to approximate the contourof Figure 9. This time we don’t need to buildnew points, because the sections between twoconsecutive corners are quite short and linear.However, we need to find a good approxima-tion for corner points: as we can see in Figure10 there is the same problem we faced in theother example, so the corner points are not in-cluded in the data set. Using R max = 20 mm and ε = 0 . (for detection of anchor points)we find 18 corner points. After applying thesecond method, we have the results in Table 3.Table 3: Results with the second method onthe wooden object with ε = 0 . .Tolerance Arcs Segments Bad Points1.5 mm 43 0 91 mm 46 0 100.5 mm 51 0 10Table 4: Results with the second method onthe wooden object with ε = 0 . .Tolerance Arcs Segments Bad Points1.5 mm 38 0 41 mm 40 0 70.5 mm 43 0 7We can notice that the number of bad points in Table 3 is quite high. This is justified by thefact that the algorithm detected only 18 cor-ners, so at least three or four bad points shouldbe corner points. Indeed, using ε = 0 . for de-tection of anchor points we find 21 corners andthe results in Table 4.Finally, we want to test the method in a sim-ple case. The contour to be approximated isrepresented in Figure 11.Figure 11: Simple contour.In this situation, the sections between two con-secutive corners are simple lines, but there isstill the problem of detecting the true corners,as we can notice in Figure 12. Using again R max = 20 mm and ε = 0 . we find five cor-ners and the results of Table 5.Table 5: Results with the second method onthe simple contour.Tolerance Arcs Segments Bad Points1.5 mm 5 0 01 mm 5 0 00.5 mm 5 0 019igure 12: Also in this case the corners are not included in the data set. Conclusions
Two methods of approximation of a planardata set were presented. Both methods arebased on arcs and segments. The first onesearches the longest arcs in the data set andeventually fills the gaps among them usingbiarcs or segments. The second one createsa lower number of arcs but the final curve ismore regular and C continuous. Moreover,the number of segments is minimized. Finally,a method of corner detection was presented.The corners of the object contour are deter-mined using circle intersections and Taubin’salgorithm. Aknowledgments
The author is very greatful to Professor C.Dagnino and Dr. D. Inaudi for their supportand their suggestions.
References [1] D. C
HETVERIKOV , A simple and efficientalgorithm for detection of high curvaturepoints in planar curves , Springer, Lecture notes in computer science, vol. 2756,746-753, 2003[2] H. F
REEMAN , L. S. D
AVIS , A corner-finding algorithm for chain-coded curves ,IEEE Transactions on computers, 26:297-303, March 1977[3] D. S. M
EEK , D. J. W
ALTON , The familyof biarcs that matches planar, two point G Hermite data , Elsevier, Journal ofComputational and Applied Mathematics212, 31-45, 2008[4] C. J. O NG , Y. S. W ONG , H. T. L OH , X.G. H ONG , An optimization approach forbiarc curve-fitting of B-spline curves , El-sevier, Computer-Aided Design, Vol. 28,951-959, 1996[5] H. P
ARK , Error-bounded biarc approxima-tion of planar curves , Elsevier, Computer-Aided Design, Vol. 36, 1241-1251, 2004[6] G. V. P
EDROSA , C. A. Z. B
ARCELOS , Anisotropic diffusion for effective shapecorner point detection , Elsevier, PatterRecognition Letters 31, 1658-1664, 2010[7] L. A. P
IEGL , W. T
ILLER , Biarc ap-proximation of NURBS curves , Elsevier,20omputer-Aided Design, Vol. 34, 807-814, 2002[8] L. A. P
IEGL , W. T
ILLER , Data approxi-mation using biarcs , Springer-Verlag, Vol.18, 59-65, 2002[9] L. A. P
IEGL , W. T
ILLER , Least-SquaresB-Spline Curve Approximation with Ar-bitrary End Derivatives , Springer-Verlag,Vol. 16, 109-116, 2000[10] L. A. P
IEGL , W. T
ILLER , The NURBS book ,Springer-Verlag, 1995[11] A. R
OSENFELD , E. J
OHNSTON , An-gle detection on digital curves , IEEETransactions on computers, 22:875-878,September 1973[12] A. R
OSENFELD , J. S. W
ESZKA , An im-proved method of angle detection on digi-tal curves , IEEE Transactions on comput-ers, 24:940-941, September 1975[13] J. R. R
OSSIGNAC , A. A. G. R
EQUICHA , Piecewise-circular curves for geometricmodeling , IBM J. Res. Develop., Vol. 31,296-313, 1987[14] X. S
ONG , M. A
IGNER , F. C
HEN , B. J ÜT - TLER , Circular spline fitting using anevolution process , Elsevier, Journal ofComputational and Applied Mathematics231, 423-433, 2008[15] G. T
AUBIN , Estimation of planar curves,surfaces, and nonplanar space curves de-fined by implicit equations with applica-tions to edge and range image segmenta-tion , IEEE Transactions on Patter Analy-sis and Machine Intelligence, Vol. 13, No.11, November 1991[16] Y. J. T
SENG , Y. D. C