- Email: [email protected]
Piecewise Bézier Curve Fitting by Multiobjective Simulated Annealing
12th IFAC Workshop on Intelligent Manufacturing Systems 12th IFAC Workshop on Intelligent Manufacturing Systems December 5-7, 2016. Austin, TX, USA 12th IFAC IFAC Workshop on Intelligent Intelligent Manufacturing Systems Systems 12th on Manufacturing December Workshop 5-7, 2016. Austin, TX, USA December December 5-7, 5-7, 2016. 2016. Austin, Austin, TX, TX, USA USA Available online at www.sciencedirect.com
ScienceDirect IFAC-PapersOnLine 49-31 (2016) 49–54
Piecewise B´ e zier Curve Fitting by Piecewise B´ e zier Curve Fitting by Piecewise B´ e zier Curve Fitting by Multiobjective Simulated Annealing Multiobjective Simulated Annealing Multiobjective Simulated Annealing ∗ 1
Edson Kenji Ueda ∗ 1 Marcos de Sales Guerra Tsuzuki ∗∗ 22 Edson Ueda ∗ 3 de Sales Guerra Tsuzuki ∗ 1 1 Marcos Rog´ eKenji rio Yugo Takimoto e Kubagawa Sato ∗∗ 44∗∗ 22 ∗ Edson Kenji Ueda Marcos de Sales Guerra ∗ 3 Andr´ Edson Kenji Ueda Marcos de Guerra Tsuzuki Tsuzuki Rog´ e rio Yugo Takimoto e Kubagawa Sato ∗Andr´ 5 Sales ∗ 3 Thiago de Castro Martins Eigi Miyagi ∗ 3 ∗Andr´ Rog´ e Takimoto e Sato∗∗ ∗∗66 44 5 Paulo Rog´ erio rio Yugo Yugo Takimoto Andr´ e Kubagawa Kubagawa Thiago de Castro Martins Paulo Eigi Miyagi ∗∗ 7 Sato ∗ 5 ∗ S´ılvio Ubertino RossoEigi Jr ∗∗ ∗ 5 Paulo ∗ 6 6 Thiago de Martins Miyagi 7 ThiagoRoberto de Castro Castro Martins Paulo Eigi Miyagi Roberto S´ ılvio Ubertino Rosso Jr ∗∗ 7 7 ∗∗ Roberto S´ S´ılvio ılvio Ubertino Ubertino Rosso Rosso Jr Jr Roberto ∗ ecnica da Universidade de S˜ ao Paulo, S˜ ao Paulo, Brazil. ∗ Escola Polit´ Escola Polit´ eecnica da Universidade de S˜ a o Paulo, S˜ a oDepartment. Paulo, Brazil. ∗ Mechatronics and Mechanical Systems Engineering ∗ Escola Polit´ cnica da Universidade de S˜ a o Paulo, S˜ a Paulo, Escola Polit´ ecnica da Universidade de S˜ aEngineering o Paulo, S˜ ao oDepartment. Paulo, Brazil. Brazil. Mechatronics and Mechanical Systems ∗∗ EstadualSystems de Santa Catarina, Brazil. Mechatronics and Mechanical Engineering Department. ∗∗ Universidade Mechatronics and Mechanical Systems Engineering Department. Universidade Estadual de Santa Catarina, Brazil. ∗∗ ∗∗ Universidade Estadual de Santa Catarina, Brazil. Universidade Estadual de Santa Catarina, Brazil. Abstract: The determination of an approximation curve from a given sequence of points is Abstract: The determination of an curve given of is an important task in CAD. This proposes an algorithm to aaadetermine a piecewise B´ezier Abstract: The determination of work an approximation approximation curve from from given sequence sequence of points points is Abstract: The determination of an approximation curve from given sequence of points is an important task in CAD. This work proposes an algorithm to determine a piecewise B´ e zier curve that approximates a sequence of proposes points. Itan is algorithm used a multiobjective simulated annealing an important task in CAD. This work to determine aa piecewise B´ eezier an important task in CAD. This work proposes an algorithm to determine piecewise B´ zier curve that approximates sequence ofbetween points. It used aa multiobjective annealing aiming at minimizing the a theis of pointssimulated and the curve, the curve that approximates a sequence points. is used simulated annealing curve that approximates adiscrepancy sequence of ofbetween points. It It isgiven used sequence a multiobjective multiobjective simulated annealing aiming at minimizing the discrepancy the given sequence of points and the curve, the curve length and the absolute difference of the curve length and length of the given sequence of aiming at minimizing the discrepancy between the given sequence of points and the curve, the aiming at minimizing the discrepancy between the given sequence of points and thesequence curve, the curve and difference of curve and length the of points. The discrepancy between the given sequence of points the of curve is determined curve length length and the the absolute absolute difference of the the curve length length andand length of the given given sequence by of curve length and the absolute difference of the curve length and length of the given sequence of points. The discrepancy between the given sequence of points and the curve is determined by the sumThe of the distance between eachgiven pointsequence from theof sequence andthe thecurve approximation curve, points. discrepancy the points and is determined by points. The discrepancy between the given sequence of points and the curve is determined by the sum of the distance between each point from the sequence and the approximation curve, and the distance from a point to the determined by an enhanced method in which the the sum of between each point from the and curve, the sum of the the distance distance between eachcurve pointis from the sequence sequence and the the approximation approximation curve, and the distance from a point to the curve is determined by an enhanced method in which the curve is discretized. and the distance from a point to the curve is determined by an enhanced method in which and the distance from a point to the curve is determined by an enhanced method in which the the curve is discretized. curve is is discretized. discretized. curve © 2016, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: B´ezier Curve, Piecewise Curve, Curve Fitting, Multiobjective Simulated Annealing, Keywords: B´ B´eezier zier Curve, Piecewise Piecewise Curve, Curve Curve Fitting, Multiobjective Multiobjective Simulated Annealing, Annealing, CAD. Keywords: Keywords: B´ezier Curve, Curve, Piecewise Curve, Curve, Curve Fitting, Fitting, Multiobjective Simulated Simulated Annealing, CAD. CAD. CAD. 1. INTRODUCTION ´ 2. BEZIER CURVE 1. INTRODUCTION INTRODUCTION ´´ 2. B EZIER CURVE 1. 2. B EZIER CURVE 1. INTRODUCTION ´ 2. BEZIER CURVE One task very important in engineering is to determine a A B´ezier curve is given by One task very important in engineering is to determine a by curve fromvery a given sequence of points, specially in reverseaa A n One task important in is A B´ B´eezier zier curve curve is is given given by One important in engineering engineering is to to determine determine given by curvetask fromvery a where given sequence ofused points, specially inobjects. reverse A B´ezier curve is n engineering, it can be to reconstruct curve from a given sequence of points, specially in reverse n P(u) = pi Bi,n (u), u ∈ [0, 1] (1) curve from a given sequence of points, specially in reverse n engineering, where it can be used to reconstruct objects. P(u) = p B u ∈ [0, 1] (1) i Bi,n (u), There are methods incan CAD that to uses curves toobjects. create engineering, where it be used reconstruct i=0 p P(u) = p (u), u ∈ [0, 1] (1) engineering, where it can be used to reconstruct objects. i i,n P(u) = B (u), u ∈ [0, 1] (1) i i,n There are methods in CAD that uses curves to create i=0 surfaces and in Tiller, and curves by determining There methods CAD that to i=0 There are are(Piegl methods in CAD 1996), that uses uses curves to create create i=0 surfaces (Piegl and Tiller, 1996), and by determining these curves theand object can1996), be reconstructed. Another where pi are the control points, n + 1 the total numbers surfaces (Piegl Tiller, and surfaces (Piegl and Tiller, and by by determining determining the control n +Bernstein 11 the total numbers where p i are these curves curves the object can1996), be reconstructed. reconstructed. Another of control points Bi,npoints, (u) is the polynomial the control points, n numbers p application is the path planning by determining points where these the object can be Another are the and control points, n+ +Bernstein 1 the the total total numbers where pii are these curvesis the object can be reconstructed. Another of control points and B (u) is the polynomial i,n application the path planning by determining points basis function that is defined as of control points and B (u) is the Bernstein polynomial (places) where vehicle had to pass through. i,n application is the path planning by determining points of control points and B (u) is the Bernstein polynomial i,n application is the path planning by determining points basis function that is defined as (places) where where the the vehicle vehicle had had to to pass pass through. through. basis function that as (places) function that is idefined defined n−i as n is (places) the vehicle had tofrom passathrough. i = 0, ..., n. (2) The taskwhere to determine a curve sequence of point basis B i,n (u) = n ui (1 − u)n−i , n i B (u) = u (1 − u) ,, ii = 0, ..., n. (2) The task to determine a curve from a sequence of point i n−i n i,n i n−i can be done by two different approaches, one by interpolaB (u) = u (1 − u) = 0, ..., n. (2) The task to determine a curve from a sequence of point i,n Bi,n (u) = ii u (1 − u) , i = 0, ..., n. (2) The task to by determine a curve from a sequence of point can be done two different approaches, one by interpola i tion (Maekawa et al., 2007; Gofuku et al., 2009; Jak´ o bczak, can be done by two different approaches, one by interpolan can be done by et two different approaches, one byJak´ interpolawhere ni is the binomial function, given by tion (Maekawa al., 2007; Gofuku et al., 2009; o bczak, 2015) and another approximation and where n is the tion (Maekawa et 2007; Gofuku 2009; oobczak, function, by given binomial tion et al., al., by 2007; Gofuku et et al., al., (Pandunata 2009; Jak´ Jak´ bczak, binomial by 2015)(Maekawa and another another by approximation (Pandunata and where where niii is is the the binomial function, given by 0given n! function, n Shamsuddin, 2010; Hasegawa et al., 2013; Adi et al., 2010). 2015) and by approximation (Pandunata and 2015) and another by approximation (Pandunata and ≡ 1. (3) = n! 0 n Shamsuddin, 2010; Hasegawa et al., 2013; Adi et al., 2010). The first determines a curve or function that pass by all the n! n 0 Shamsuddin, 2010; Hasegawa et al., 2013; Adi et al., 2010). i!(n − i)! i ≡ 1. (3) = n! 0 ≡ 1. n = Shamsuddin, 2010; Hasegawa etfunction al., 2013; Adipass et al., 2010). (3) The first determines a curve or that by all the i!(n − i)! 0 i ≡ 1. (3) = points, while the second determines a curve that is closer The first determines a curve or function that pass by all the i i!(n − i)! 0 The firstwhile determines a curve or function that pass by all the i!(n − i)! 0 i points, the second determines a curve that is closer to a sequence of points. This work presents a new method 2.1 Piecewise B´ e zier Curve Segments points, while second determines aa curve is closer points, while the the secondThis determines curveathat that is closer 2.1 Piecewise B´ezier Curve Segments to curve sequence of points. points. work presents new method of approximation that work uses piecewise eziermethod curves 2.1 to aa sequence of This presents aaB´ new Piecewise to a sequence of points. This work presents new method Piecewise B´ B´eezier zier Curve Curve Segments Segments of curve curve approximation that uses uses piecewisetoB´ B´ ezier zier curves curves 2.1 and a multiobjective simulated annealing approximate The characteristic of each control point of a B´ezier curve of approximation that piecewise e of curve approximation that uses piecewise B´ e zier curves multiobjective simulated annealing to approximate approximate The characteristic ofcan each control point of aa piecewise B´eezier zier curve aand sequence of points.simulated This workannealing is structured as follow. The modify all the curveof becontrol avoidedpoint by using cuand aa to each of B´ and a multiobjective multiobjective simulated annealing to approximate The characteristic characteristic ofcan each control point of a piecewise B´ezier curve curve a sequence of points. This work is structured as follow. modify all the curve be avoided by using cuIn section 2 is done a brief review of B´ e zier curve, the bic B´ e zier curve segments. This sequence of curve segment aa sequence of points. This work is structured as follow. modify all the curve can be avoided by using piecewise cusequence of points. This work is structured as follow. modify all the curve can be avoided by using piecewise cuIn section 2 is done a brief review of B´ e zier curve, the bic B´ e zier curve segments. This sequence of curve segment 1 curve fitting22 problem and piecewise cubic eziercurve, curve.the In bic is done with a weak C 1 continuity with of the surrounding In section is done aa brief review of B´ eeB´ zier B´ eezier curve segments. This sequence curve segment In section is done brief review of B´ zier curve, the bic B´ zier curve segments. This sequence of curve segment curve fitting fitting problem all and piecewise cubic to B´eebe zierminimized curve. In In is done where with aathe weak C 11 continuity with the surrounding section 3 is presented thepiecewise cost functions derivatives have the same but curve problem and cubic B´ zier curve. is with weak C with the surrounding curve fitting problem all and piecewise cubic to B´ebe zierminimized curve. In curves, is done done where with athe weak C continuity continuity with thedirection, surrounding section 3 is presented the cost functions curves, derivatives have the same direction, but as well as the problem parametrization. In section 4 the not necessary the same intensity. To create these sequence section 33 is presented all the cost functions to be minimized curves, where derivatives have the same direction, but section is presented all the cost functions to be minimized curves, where the the same derivatives haveTo the samethese direction, but as well well and as the the problem parametrization. In section 44 the the not necessary intensity. create sequence results finally a conclusion in section 5. of curves, the last control point of a curve is the first as as problem parametrization. In section not necessary necessary the the same same intensity. intensity. To To create create these these sequence sequence as well as the problem parametrization. In section 4 the not results and finally a conclusion in section 5. of curves, the last control point of a curve is the first control point the control following curve, to determinate results of curves, the last point of curve is results and and finally finally aa conclusion conclusion in in section section 5. 5. of curves, theof last control point of aaand curve is the the first first 1 email: [email protected]. control point of the following curve, and to determinate the second control point of a curve segment it is used the control point of the following curve, and to determinate 1 control point of the following curve, and to determinate 2 email: [email protected]. the second control point of aa curve segment it is used the 1 email: [email protected]. email: [email protected]. 1 equation the second control point of curve segment it is used 2 email: [email protected]. 3 the second control point of a curve segment it is used the the email: [email protected]. 2 email: [email protected]. equation email: [email protected]. 2 3 email: [email protected]. equation 4 [email protected]. p3i+4 = p3i+3 − (p3i+2 − p3i+3 ) ∗ β3i+4 , (4) equation 3 email: email: [email protected]. email: [email protected]. 3 4 p p (p p ∗∗ β (4) email: [email protected]. 3i+4 = 3i+3 − 3i+2 − 3i+3 ) 3i+4 ,, 5 email: [email protected]. 4 email: [email protected]. p = p − (p − p ) β (4) 3i+4 3i+3 3i+2 3i+3 3i+4 email: [email protected]. [email protected]. 4 p3i+4 = p3i+3 − (p3i+2 − p3i+3 ) ∗ β3i+4 , (4) 5 6 email: email: [email protected]. [email protected]. where i = 0, . . . , j − 1 and j − 1 is the number of segments 5 email: email: [email protected]. 5 6 email: [email protected]. where i = 0, . . . , j − 1 and j − 1 is the number of segments 7 email: [email protected]. 6 email: and βi iiis= email: [email protected]. [email protected]. continuity An11 example is shown in Fig. 1, where .. .. .. ,, jj − and is of 6 where =aa 0, 0, − 11factor. and jj − − is the the number number of segments segments 7 email: and β email: [email protected]. [email protected]. continuity factor. An example is shown in Fig. 1, i is 7 and β [email protected]. is a continuity factor. An example is shown 7 email: i and βi is a continuity factor. An example is shown in email: [email protected]. in Fig. Fig. 1, 1,
Copyright@ 57 Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © 2016 2016, IFAC IFAC (International Federation of Automatic Control) Copyright@ 2016 IFAC 57 Copyright@ 2016 IFAC 57 Peer review under responsibility of International Federation of Automatic Copyright@ 2016 IFAC 57 Control. 10.1016/j.ifacol.2016.12.160
2017 IFAC IMS 50 December 5-7, 2016. Austin, TX, USA
p1
Edson Kenji Ueda et al. / IFAC-PapersOnLine 49-31 (2016) 49–54
p2
p1 d1 , u1 = 0.24
p3 p0
P(u1 ) d2 , u2 = 0.47 p3 , P(u4 ), d4 , u4 = 1
p4
p0 , P(u0 ), d0 , u0 = 0
P(u2 ) P(u3 )
p5
p9
d3 , u3 = 0.76 Control Point (p) B´ ezier Curve (P)
p8
Discretized Points of the B´ ezier Curve (u)
p6
Data Points (d)
p7
Curve Error
Fig. 1. Sequence of three cubic B´ezier curve. Points p3 and p6 are transition points between B´ezier curves.
Fig. 2. Distance determined by an approximated method. The methods proposed by Hasegawa et al. (2013) and Pandunata and Shamsuddin (2010) could wrongly determine the points as shown in Fig. 2. Maekawa et al. (2007) has some advantages, although its processing cost is higher. On the other hand, the method proposed by Tavares et al. (2011) is the one with more advantages due to the balance between precision and processing time. In this method, the higher the discretization the better the quality and the processing time increases.
that is a composition of three segments of cubic B´ezier curve. 2.2 B´ezier curve fitting The objective of curve fitting problem is to define a curve that approximates a given sequence of points by determining its control points (p0 , ..., pn ). This determination is done by minimizing the discrepancy between the curve and the given sequence of points, which can be calculated as the distance between the curve and each point of the given sequence of points. The function that describes the function to be minimized is f (p0 , ..., pn ) =
m−1 k=1
p2
|dk − P(uk )|2
3. COST FUNCTIONS Equation (5) is the first cost function adopted. This expression evaluate just the discrepancy between curve and given points. Other functions could be evaluated. In this work it will be also evaluated the curve length and the absolute difference between the length of the curve and the given sequence of points. Once (5) defines a family of curves, the other two functions will evaluate the difference among the curves of the family.
(5)
where 0 < k < m , n+1 is the number of control points and d is the sequence of points with m+1 elements. |dk −P(uk )| is the distance between one point dk and the B´ezier curve. uk is the parameter value that defines which point of the curve is the closest to the point dk , which also can be interpreted as the projection of the point dk over the curve. Points d0 and dm are not added in the sum, once the B´ezier curve always interpolates the first and last control points.
The method to determine the discrepancy between point and curve is an enhancement of the method proposed by Tavares et al. (2011). In the determination of the point of a discretized curve that is closer to a given point dk , there are considered two triangles as shown is Fig. 3. Squares points are points of the discretized curve and circle point is a point of the given sequence of points.
There are several methods to determine the parameter uk . Hasegawa et al. (2013) and Pandunata and Shamsuddin (2010) determined this parameter in an approximated approach by using the chord length parametrization. This parametrization approximates the uk parameter as a fraction between the sum of the chords until the respective point dk and the sum of all the chords. Maekawa et al. (2007) used the Newton-Raphson method to determine the parameter uk which makes the segment that connect the points P(uk ) and dk orthogonal to the tangent at P(uk ) point. Tavares et al. (2011) used a method in which the curve is discretized with a sequence of points and all the points are measured to the given point dk and it is determined the point of the discretized curve that is closer to dk .
Firstly, a search is done to determine which point of the discretized curve is closer to the point of the sequence dk . This point is named P(uk ). As shown in Fig. 3, two triangles are defined. The first is defined by the points dk , P(uk ) and P(uk−1 ); the second is defined by dk , P(uk ) and P(uk+1 ). Then it is evaluated which triangle has the smallest area. Next, it is made a verification whether the triangle with smallest area is acute or obtuse, in specific the angles between the segments dk P(uk ) and P(uk )P(uk−1 ) for the first triangle, and dk P(uk ) and P(uk )P(uk+1 ) for the second. If the triangle with the smallest area is acute, it is calculated the height of the triangle related to the segment that is defined by two points of the discretized curve, and 58
2017 IFAC IMS December 5-7, 2016. Austin, TX, USA
Edson Kenji Ueda et al. / IFAC-PapersOnLine 49-31 (2016) 49–54
51
dk
Sequence of Points
P(uk−1 )
P(uk )
P10 P11
...
1st Curve
P1v-1 P1v
...
... j-2 Curves
...
Bézier curve
(b)
Fig. 5. (a) First limit to position of new curve transition; (b) Second limit to position of new curve transition. the dimension of a vector that is the parameter β of the expression 4. This algorithm can also process discrete parameters which is the integer numbers that represents an index associated with points of the sequence that there is a transition between curves.
dm-4 dm-3 dm-2 dm-1 dm
Pj0 Pj1
Control Point
Bézier curve
(a)
P(uk+1 )
Fig. 3. After finding the point P(uk ) of the discretized curve that is closer to the point dk , two triangles are defined. d0 d1 d2 d3 d4
Sequence of Points
Control Point
Using Fig. 1 as an example, that shows a composition of three cubic B´ezier curve, the continuous parameters are:
Pj1v-1 Pjv
(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)
Last Curve
Fig. 4. Local search for the closest point. this height is the distance between curve and point. If the triangle with the smallest area is obtuse, than it is used the other triangle to determine the distance between curve and point. And, if both triangles are obtuse, neither of them are used and it is adopted the distance between dk and P(uk ) as being the distance between discretized curve and point.
Coordinate x of point p1 (p1 x) Coordinate y of point p1 (p1 y) Coordinate x of point p2 (p2 x) Coordinate y of point p2 (p2 y) Continuity factor of point p4 (β4 ) Coordinate x of point p5 (p5 x) Coordinate y of point p5 (p5 y) Continuity factor of point p7 (β7 ) Coordinate x of point p8 (p8 x) Coordinate y of point p8 (p8 y)
and the discrete parameters are two index that are associated with some point of the given sequence, in which there is a transition of curves.
Another improvement is done in the first step, where it is searched the point of the discretized curve that is closer to a point of the sequence. Once the piecewise B´ezier curve is used, each segment approximates only a section of the sequence. Then the search is locally done by searching each point only in one segment as shown in the representation of the Fig. 4. In this example, points d5 and dm−4 are points in which there is a transition of B´ezier curves. So, for points d1 , d2 , d3 and d4 , the search for the closest point is done only in the first curve defined by points P11 , . . . , P1v−1 . And for points dm−4 , dm−3 , dm−2 and dm−1 , the search is done only in the last curve defined by points Pj1 , . . . , Pjv−1 .
3.2 Next Candidate Generation To generate a new candidate by modifying a continuous parameter it is used the method developed by Corana et al. (1987) that is given by xk+1 = xk + v · ∆ri · ei
(6)
where ∆ri is the step size in direction ei and v is an uniform random number in [-1, 1]. To modify a discrete parameter, it is defined a simple heuristic. Initially, the curves division are determined according to the division of the number of points of the sequence and the number of curve segments. So, the selection of which point is modified is not random. To modify an index associated with a transition point, it is generated a random value between [0, 1], and if this value is lower than 0.5, the index is decreased by one unity, otherwise it is increased by one unity. This operation modifies which point of the given sequence of points is going to be used as a transition point between two B´ezier curves. This heuristic is limited by two restrictions:
The determination of both lengths can be easily calculated. The length of the given sequence of points is calculated by the sum of all the segments that is determined by consecutive points. And, once the curves were discretized to determine the discrepancy, the length of the curve is determined by the sum of all the segments of the discretized curve. 3.1 Problem Parametrization It is used AMOSA (Archived MultiObjective Simulated Annealing) (Bandyopadhyay et al., 2008) a multiobjective simulated annealing to solve the problem of curve fitting. This algorithm can process very well the continuous variables, such as the control points of the curves and
• the change of position can not be more than four unities from the initial position. This change could not happen in the increment and decrease of the 59
2017 IFAC IMS 52 December 5-7, 2016. Austin, TX, USA
Edson Kenji Ueda et al. / IFAC-PapersOnLine 49-31 (2016) 49–54
index. This rule is adopted to avoid the creation of candidates that changes too much the initial and final position of one B´ezier curve segment. This limit is shown in Fig. 5(a), in which the blue square is the initial position of joint of two segments, and this position can be modified to any green square. • each parameter must have a minimum distance to the next and previous parameter, this limit is done to avoid the creation of too short curve segments that approximates only few points of the sequence. It is adopted a minimum distance of 4 unities, once cubic curves is used to approximate the sequence. And to approximates few points is a waste of computational time. This limit is shown in Fig. 5(b), in which the blue square is the initial position of the joint of two segments and this position can be modified to any green square. Notice that respecting the first limit, the change could be greater, however due to the fact that the end of the second segment is too close, the position of change is shorter. 1 Curve 2 Curves 3 Curves
55 F2
(a)
50
(b)
Fig. 7. Fitting curve for the parabola example with 1, 2 and 3 curve segments. Each segment is in different colors and its control points are diamond shaped in the respective color and the given sequence of points is the black circles. (a) The curve of lower discrepancy of F1 x F2; (b) The curve of lower discrepancy of F1 x F3.
45 0
2,000 4,000 6,000 F1 (a)
120
F3
0.6 0.4
F2
1 Curve 2 Curves 3 Curves
5 Curves 6 Curves 7 Curves
100 80
0.2 0
0 0
10
20
1,000
2,000
F1
30
(a)
F1 (b)
1 F3
Fig. 6. Pareto front of the parabola example. (a) Curve discrepancy (F1) x fitting curve length (F2); (b) Curve discrepancy (F1) x absolute difference between the length of the fitting curve and the length of the sequence of points (F3).
5 Curves 6 Curves 7 Curves
0.5 0 0
4. RESULTS
20
40
60
80
F1 The developed method was tested in three examples. In each example the curve with desired format is created by discretizing it with 50 points and a noise is added in each coordinate. The chosen formats to be tested are a parabola, a curve with inflection and a curve with auto-intersection. To determine the curve discrepancy each curve segment was discretized with 100 points. The objective functions to be minimized are the curve discrepancy
(b)
Fig. 8. Pareto front of the example of the curve with auto-intersection. (a) Curve discrepancy (F1) x fitting curve length (F2); (b) Curve discrepancy (F1) x absolute difference between the length of the fitting curve and the length of the sequence of points (F3).
60
2017 IFAC IMS December 5-7, 2016. Austin, TX, USA
Edson Kenji Ueda et al. / IFAC-PapersOnLine 49-31 (2016) 49–54
53
both images, it is possible to observe that the developed method successfully approximates the given sequence of points with 1, 2 and 3 curves.
(F1), the fitting curve length (F2) and the absolute difference between the length of the fitting curve and the length of the sequence of points (F3). The AMOSA method was applied to F1 with F2; and to F1 with F3. Each example is tested with three different number of curve segments.
The second example is a curve with auto-intersection and it is used 5, 6 and 7 curve segments to approximate the given sequence of points. Fig. 8 shows the Pareto fronts. In Fig. 8(a), it is possible to observe once again that the increase of the discrepancy decreases the fitting curve length. However, in this example, the higher the number of curves the higher the curve length. Fig. 8(b) shows that there is almost no relation between the discrepancy and the absolute difference of lengths. There is no perceptive change in the discrepancy with a change in the difference of lengths. Fig. 9 shows the curves with lower discrepancy by minimizing F1 and F2 in Fig. 9(a); and minimizing F1 and F3 in Fig. 9(b). 2 Curves 3 Curves 4 Curves
F2
120
(a)
100 80 60 40
0
0.5
1
1.5 ·10
F1 (a)
F3
0.6
2 4
2 Curves 3 Curves 4 Curves
0.4 0.2 0 5
10
15
20
F1 (b) (b)
Fig. 10. Pareto front of the example of the curve with inflection. (a) Curve discrepancy (F1) x fitting curve length (F2); (b) Curve discrepancy (F1) x absolute difference between the length of the fitting curve and the length of the sequence of points (F3).
Fig. 9. Fitting curve for the auto-intersection example with 5, 6 and 7 curves. Each segment is in different colors and its control points are diamond shaped in the respective color and the given sequence of points is the black circles. (a) The curve of lower discrepancy of F1 x F2; (b) The curve of lower discrepancy of F1 x F3.
The third example is a curve with inflection and it is used 2, 3 and 4 curve segments to approximate the given sequence of points. Fig. 10 shows the Pareto fronts. Fig. 10(a) shows the same behavior of the Pareto front of the example of the curve with auto-intersection. Fig. 10(b) shows again that there is almost no relation between the discrepancy and the absolute difference of lengths. However, in this example, the higher the number of curves the lower the minimum value of discrepancy. Fig. 11 shows the curves with lower discrepancy by minimizing F1 and F2 in Fig. 11(a); and minimizing F1 and F3 in Fig. 11(b). In Fig. 11(a) all the curves are successfully created, but all curves in Fig. 11(b) a slight hook is presented at the end of the last piece of curve. These hooks are created because
The first example is a parabola and it is used 1,2 and 3 curve segments to approximate the given sequence of points. Fig. 6 shows the Pareto fronts and it is possible to observe from Fig. 6(a) that the lower the discrepancy the higher the length of the fitting curve. The number of curves has low influence over the discrepancy due to the simplicity of this example. Fig. 6(b) shows that by minimizing the discrepancy you increase the difference of lengths, however after a value there is no change of the difference of lengths when you increase the discrepancy. Fig. 7 shows the curves with lower discrepancy by minimizing F1 and F2 in Fig. 7(a); and minimizing F1 and F3 in Fig. 7(b). In 61
2017 IFAC IMS 54 December 5-7, 2016. Austin, TX, USA
Edson Kenji Ueda et al. / IFAC-PapersOnLine 49-31 (2016) 49–54
5. CONCLUSION
the length of the sequence of points is slightly longer than the desired curve.
An algorithm to solve the curve fitting problem was developed, that used piecewise cubic B´ezier curves with a weak C 1 continuity and a multiobjective simulated annealing. The cost functions analyzed were the fitting curve discrepancy, the fitting curve length and the absolute difference of the fitting curve length and given points length. These developed methods were tested and applied in three different examples and almost all the examples were successfully approximated. As future works, the developed method used several numbers of curves, however it is necessary to study a method to determine the optimum number of curves. ACKNOWLEDGEMENTS MSG Tsuzuki and TC Martins were partially supported by CNPq (grants 310.663/2013-0 and 306.415/2012-7). RY Takimoto and EK Ueda were partially supported by CNPq (grants 150508/2015-8 and 140518/2015-0). AK Sato was supported by CAPES/PNPD. REFERENCES (a)
Adi, D., Shamsuddin, S., and Hashim, S.Z.M. (2010). Nurbs curve approximation using particle swarm optimization. In Computer Graphics, Imaging and Visualization (CGIV), 2010 Seventh International Conference on, 73–79. Bandyopadhyay, S., Saha, S., Maulik, U., and Deb, K. (2008). A simulated annealing-based multiobjective optimization algorithm: AMOSA. Evolutionary Computation, IEEE Transactions on, 12(3), 269–283. Corana, A., Marchesi, M., Martini, C., and Ridella, S. (1987). Minimizing multimodal functions of continuous variables with the simulated annealing algorithm. ACM Transactions on Mathematical Software, 13, 262–280. Gofuku, S., Tamura, S., and Maekawa, T. (2009). Pointtangent/point-normal b-spline curve interpolation by geometric algorithms. Computer-Aided Design, 41, 412– 422. Hasegawa, A.Y., Jr., R.S.U.R., and Tsuzuki, M.S.G. (2013). Differential evolution optimization for bezier curve fitting. In 11th IFAC Workshop on Intelligent Manufacturing Systems, 233–238. Jak´obczak, D.J. (2015). Probabilistic 2D point interpolation and extrapolation via data modeling. Informatica (Slovenia), 39(1), 53–61. Maekawa, T., Matsumoto, Y., and Namik, K. (2007). Interpolation by geometric algorithm. Computer-Aided Design, 39, 313–323. Pandunata, P. and Shamsuddin, S. (2010). Differential evolution optimization for bezier curve fitting. In Computer Graphics, Imaging and Visualization (CGIV), 2010 Seventh International Conference on, 68–72. Piegl, L. and Tiller, W. (1996). Algorithm for approximate nurbs skinning. Computer-Aided Design, 28(9), 699 – 706. Tavares, R., Martins, T., and Tsuzuki, M. (2011). Simulated annealing with adaptive neighborhood: A case study in off-line robot path planning. Expert Systems with Applications, 38(4), 2951 – 2965.
(b)
Fig. 11. Fitting curve for the inflection example with 2, 3 and 4 curve segments. Each segment is in different colors and its control points are diamond shaped in the respective color and the given sequence of points is the black circles. (a) The curve of lower discrepancy of F1 x F2; (b) The curve of lower discrepancy of F1 x F3.
62
ScienceDirect IFAC-PapersOnLine 49-31 (2016) 49–54
Piecewise B´ e zier Curve Fitting by Piecewise B´ e zier Curve Fitting by Piecewise B´ e zier Curve Fitting by Multiobjective Simulated Annealing Multiobjective Simulated Annealing Multiobjective Simulated Annealing ∗ 1
Edson Kenji Ueda ∗ 1 Marcos de Sales Guerra Tsuzuki ∗∗ 22 Edson Ueda ∗ 3 de Sales Guerra Tsuzuki ∗ 1 1 Marcos Rog´ eKenji rio Yugo Takimoto e Kubagawa Sato ∗∗ 44∗∗ 22 ∗ Edson Kenji Ueda Marcos de Sales Guerra ∗ 3 Andr´ Edson Kenji Ueda Marcos de Guerra Tsuzuki Tsuzuki Rog´ e rio Yugo Takimoto e Kubagawa Sato ∗Andr´ 5 Sales ∗ 3 Thiago de Castro Martins Eigi Miyagi ∗ 3 ∗Andr´ Rog´ e Takimoto e Sato∗∗ ∗∗66 44 5 Paulo Rog´ erio rio Yugo Yugo Takimoto Andr´ e Kubagawa Kubagawa Thiago de Castro Martins Paulo Eigi Miyagi ∗∗ 7 Sato ∗ 5 ∗ S´ılvio Ubertino RossoEigi Jr ∗∗ ∗ 5 Paulo ∗ 6 6 Thiago de Martins Miyagi 7 ThiagoRoberto de Castro Castro Martins Paulo Eigi Miyagi Roberto S´ ılvio Ubertino Rosso Jr ∗∗ 7 7 ∗∗ Roberto S´ S´ılvio ılvio Ubertino Ubertino Rosso Rosso Jr Jr Roberto ∗ ecnica da Universidade de S˜ ao Paulo, S˜ ao Paulo, Brazil. ∗ Escola Polit´ Escola Polit´ eecnica da Universidade de S˜ a o Paulo, S˜ a oDepartment. Paulo, Brazil. ∗ Mechatronics and Mechanical Systems Engineering ∗ Escola Polit´ cnica da Universidade de S˜ a o Paulo, S˜ a Paulo, Escola Polit´ ecnica da Universidade de S˜ aEngineering o Paulo, S˜ ao oDepartment. Paulo, Brazil. Brazil. Mechatronics and Mechanical Systems ∗∗ EstadualSystems de Santa Catarina, Brazil. Mechatronics and Mechanical Engineering Department. ∗∗ Universidade Mechatronics and Mechanical Systems Engineering Department. Universidade Estadual de Santa Catarina, Brazil. ∗∗ ∗∗ Universidade Estadual de Santa Catarina, Brazil. Universidade Estadual de Santa Catarina, Brazil. Abstract: The determination of an approximation curve from a given sequence of points is Abstract: The determination of an curve given of is an important task in CAD. This proposes an algorithm to aaadetermine a piecewise B´ezier Abstract: The determination of work an approximation approximation curve from from given sequence sequence of points points is Abstract: The determination of an approximation curve from given sequence of points is an important task in CAD. This work proposes an algorithm to determine a piecewise B´ e zier curve that approximates a sequence of proposes points. Itan is algorithm used a multiobjective simulated annealing an important task in CAD. This work to determine aa piecewise B´ eezier an important task in CAD. This work proposes an algorithm to determine piecewise B´ zier curve that approximates sequence ofbetween points. It used aa multiobjective annealing aiming at minimizing the a theis of pointssimulated and the curve, the curve that approximates a sequence points. is used simulated annealing curve that approximates adiscrepancy sequence of ofbetween points. It It isgiven used sequence a multiobjective multiobjective simulated annealing aiming at minimizing the discrepancy the given sequence of points and the curve, the curve length and the absolute difference of the curve length and length of the given sequence of aiming at minimizing the discrepancy between the given sequence of points and the curve, the aiming at minimizing the discrepancy between the given sequence of points and thesequence curve, the curve and difference of curve and length the of points. The discrepancy between the given sequence of points the of curve is determined curve length length and the the absolute absolute difference of the the curve length length andand length of the given given sequence by of curve length and the absolute difference of the curve length and length of the given sequence of points. The discrepancy between the given sequence of points and the curve is determined by the sumThe of the distance between eachgiven pointsequence from theof sequence andthe thecurve approximation curve, points. discrepancy the points and is determined by points. The discrepancy between the given sequence of points and the curve is determined by the sum of the distance between each point from the sequence and the approximation curve, and the distance from a point to the determined by an enhanced method in which the the sum of between each point from the and curve, the sum of the the distance distance between eachcurve pointis from the sequence sequence and the the approximation approximation curve, and the distance from a point to the curve is determined by an enhanced method in which the curve is discretized. and the distance from a point to the curve is determined by an enhanced method in which and the distance from a point to the curve is determined by an enhanced method in which the the curve is discretized. curve is is discretized. discretized. curve © 2016, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: B´ezier Curve, Piecewise Curve, Curve Fitting, Multiobjective Simulated Annealing, Keywords: B´ B´eezier zier Curve, Piecewise Piecewise Curve, Curve Curve Fitting, Multiobjective Multiobjective Simulated Annealing, Annealing, CAD. Keywords: Keywords: B´ezier Curve, Curve, Piecewise Curve, Curve, Curve Fitting, Fitting, Multiobjective Simulated Simulated Annealing, CAD. CAD. CAD. 1. INTRODUCTION ´ 2. BEZIER CURVE 1. INTRODUCTION INTRODUCTION ´´ 2. B EZIER CURVE 1. 2. B EZIER CURVE 1. INTRODUCTION ´ 2. BEZIER CURVE One task very important in engineering is to determine a A B´ezier curve is given by One task very important in engineering is to determine a by curve fromvery a given sequence of points, specially in reverseaa A n One task important in is A B´ B´eezier zier curve curve is is given given by One important in engineering engineering is to to determine determine given by curvetask fromvery a where given sequence ofused points, specially inobjects. reverse A B´ezier curve is n engineering, it can be to reconstruct curve from a given sequence of points, specially in reverse n P(u) = pi Bi,n (u), u ∈ [0, 1] (1) curve from a given sequence of points, specially in reverse n engineering, where it can be used to reconstruct objects. P(u) = p B u ∈ [0, 1] (1) i Bi,n (u), There are methods incan CAD that to uses curves toobjects. create engineering, where it be used reconstruct i=0 p P(u) = p (u), u ∈ [0, 1] (1) engineering, where it can be used to reconstruct objects. i i,n P(u) = B (u), u ∈ [0, 1] (1) i i,n There are methods in CAD that uses curves to create i=0 surfaces and in Tiller, and curves by determining There methods CAD that to i=0 There are are(Piegl methods in CAD 1996), that uses uses curves to create create i=0 surfaces (Piegl and Tiller, 1996), and by determining these curves theand object can1996), be reconstructed. Another where pi are the control points, n + 1 the total numbers surfaces (Piegl Tiller, and surfaces (Piegl and Tiller, and by by determining determining the control n +Bernstein 11 the total numbers where p i are these curves curves the object can1996), be reconstructed. reconstructed. Another of control points Bi,npoints, (u) is the polynomial the control points, n numbers p application is the path planning by determining points where these the object can be Another are the and control points, n+ +Bernstein 1 the the total total numbers where pii are these curvesis the object can be reconstructed. Another of control points and B (u) is the polynomial i,n application the path planning by determining points basis function that is defined as of control points and B (u) is the Bernstein polynomial (places) where vehicle had to pass through. i,n application is the path planning by determining points of control points and B (u) is the Bernstein polynomial i,n application is the path planning by determining points basis function that is defined as (places) where where the the vehicle vehicle had had to to pass pass through. through. basis function that as (places) function that is idefined defined n−i as n is (places) the vehicle had tofrom passathrough. i = 0, ..., n. (2) The taskwhere to determine a curve sequence of point basis B i,n (u) = n ui (1 − u)n−i , n i B (u) = u (1 − u) ,, ii = 0, ..., n. (2) The task to determine a curve from a sequence of point i n−i n i,n i n−i can be done by two different approaches, one by interpolaB (u) = u (1 − u) = 0, ..., n. (2) The task to determine a curve from a sequence of point i,n Bi,n (u) = ii u (1 − u) , i = 0, ..., n. (2) The task to by determine a curve from a sequence of point can be done two different approaches, one by interpola i tion (Maekawa et al., 2007; Gofuku et al., 2009; Jak´ o bczak, can be done by two different approaches, one by interpolan can be done by et two different approaches, one byJak´ interpolawhere ni is the binomial function, given by tion (Maekawa al., 2007; Gofuku et al., 2009; o bczak, 2015) and another approximation and where n is the tion (Maekawa et 2007; Gofuku 2009; oobczak, function, by given binomial tion et al., al., by 2007; Gofuku et et al., al., (Pandunata 2009; Jak´ Jak´ bczak, binomial by 2015)(Maekawa and another another by approximation (Pandunata and where where niii is is the the binomial function, given by 0given n! function, n Shamsuddin, 2010; Hasegawa et al., 2013; Adi et al., 2010). 2015) and by approximation (Pandunata and 2015) and another by approximation (Pandunata and ≡ 1. (3) = n! 0 n Shamsuddin, 2010; Hasegawa et al., 2013; Adi et al., 2010). The first determines a curve or function that pass by all the n! n 0 Shamsuddin, 2010; Hasegawa et al., 2013; Adi et al., 2010). i!(n − i)! i ≡ 1. (3) = n! 0 ≡ 1. n = Shamsuddin, 2010; Hasegawa etfunction al., 2013; Adipass et al., 2010). (3) The first determines a curve or that by all the i!(n − i)! 0 i ≡ 1. (3) = points, while the second determines a curve that is closer The first determines a curve or function that pass by all the i i!(n − i)! 0 The firstwhile determines a curve or function that pass by all the i!(n − i)! 0 i points, the second determines a curve that is closer to a sequence of points. This work presents a new method 2.1 Piecewise B´ e zier Curve Segments points, while second determines aa curve is closer points, while the the secondThis determines curveathat that is closer 2.1 Piecewise B´ezier Curve Segments to curve sequence of points. points. work presents new method of approximation that work uses piecewise eziermethod curves 2.1 to aa sequence of This presents aaB´ new Piecewise to a sequence of points. This work presents new method Piecewise B´ B´eezier zier Curve Curve Segments Segments of curve curve approximation that uses uses piecewisetoB´ B´ ezier zier curves curves 2.1 and a multiobjective simulated annealing approximate The characteristic of each control point of a B´ezier curve of approximation that piecewise e of curve approximation that uses piecewise B´ e zier curves multiobjective simulated annealing to approximate approximate The characteristic ofcan each control point of aa piecewise B´eezier zier curve aand sequence of points.simulated This workannealing is structured as follow. The modify all the curveof becontrol avoidedpoint by using cuand aa to each of B´ and a multiobjective multiobjective simulated annealing to approximate The characteristic characteristic ofcan each control point of a piecewise B´ezier curve curve a sequence of points. This work is structured as follow. modify all the curve be avoided by using cuIn section 2 is done a brief review of B´ e zier curve, the bic B´ e zier curve segments. This sequence of curve segment aa sequence of points. This work is structured as follow. modify all the curve can be avoided by using piecewise cusequence of points. This work is structured as follow. modify all the curve can be avoided by using piecewise cuIn section 2 is done a brief review of B´ e zier curve, the bic B´ e zier curve segments. This sequence of curve segment 1 curve fitting22 problem and piecewise cubic eziercurve, curve.the In bic is done with a weak C 1 continuity with of the surrounding In section is done aa brief review of B´ eeB´ zier B´ eezier curve segments. This sequence curve segment In section is done brief review of B´ zier curve, the bic B´ zier curve segments. This sequence of curve segment curve fitting fitting problem all and piecewise cubic to B´eebe zierminimized curve. In In is done where with aathe weak C 11 continuity with the surrounding section 3 is presented thepiecewise cost functions derivatives have the same but curve problem and cubic B´ zier curve. is with weak C with the surrounding curve fitting problem all and piecewise cubic to B´ebe zierminimized curve. In curves, is done done where with athe weak C continuity continuity with thedirection, surrounding section 3 is presented the cost functions curves, derivatives have the same direction, but as well as the problem parametrization. In section 4 the not necessary the same intensity. To create these sequence section 33 is presented all the cost functions to be minimized curves, where derivatives have the same direction, but section is presented all the cost functions to be minimized curves, where the the same derivatives haveTo the samethese direction, but as well well and as the the problem parametrization. In section 44 the the not necessary intensity. create sequence results finally a conclusion in section 5. of curves, the last control point of a curve is the first as as problem parametrization. In section not necessary necessary the the same same intensity. intensity. To To create create these these sequence sequence as well as the problem parametrization. In section 4 the not results and finally a conclusion in section 5. of curves, the last control point of a curve is the first control point the control following curve, to determinate results of curves, the last point of curve is results and and finally finally aa conclusion conclusion in in section section 5. 5. of curves, theof last control point of aaand curve is the the first first 1 email: [email protected]. control point of the following curve, and to determinate the second control point of a curve segment it is used the control point of the following curve, and to determinate 1 control point of the following curve, and to determinate 2 email: [email protected]. the second control point of aa curve segment it is used the 1 email: [email protected]. email: [email protected]. 1 equation the second control point of curve segment it is used 2 email: [email protected]. 3 the second control point of a curve segment it is used the the email: [email protected]. 2 email: [email protected]. equation email: [email protected]. 2 3 email: [email protected]. equation 4 [email protected]. p3i+4 = p3i+3 − (p3i+2 − p3i+3 ) ∗ β3i+4 , (4) equation 3 email: email: [email protected]. email: [email protected]. 3 4 p p (p p ∗∗ β (4) email: [email protected]. 3i+4 = 3i+3 − 3i+2 − 3i+3 ) 3i+4 ,, 5 email: [email protected]. 4 email: [email protected]. p = p − (p − p ) β (4) 3i+4 3i+3 3i+2 3i+3 3i+4 email: [email protected]. [email protected]. 4 p3i+4 = p3i+3 − (p3i+2 − p3i+3 ) ∗ β3i+4 , (4) 5 6 email: email: [email protected]. [email protected]. where i = 0, . . . , j − 1 and j − 1 is the number of segments 5 email: email: [email protected]. 5 6 email: [email protected]. where i = 0, . . . , j − 1 and j − 1 is the number of segments 7 email: [email protected]. 6 email: and βi iiis= email: [email protected]. [email protected]. continuity An11 example is shown in Fig. 1, where .. .. .. ,, jj − and is of 6 where =aa 0, 0, − 11factor. and jj − − is the the number number of segments segments 7 email: and β email: [email protected]. [email protected]. continuity factor. An example is shown in Fig. 1, i is 7 and β [email protected]. is a continuity factor. An example is shown 7 email: i and βi is a continuity factor. An example is shown in email: [email protected]. in Fig. Fig. 1, 1,
Copyright@ 57 Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © 2016 2016, IFAC IFAC (International Federation of Automatic Control) Copyright@ 2016 IFAC 57 Copyright@ 2016 IFAC 57 Peer review under responsibility of International Federation of Automatic Copyright@ 2016 IFAC 57 Control. 10.1016/j.ifacol.2016.12.160
2017 IFAC IMS 50 December 5-7, 2016. Austin, TX, USA
p1
Edson Kenji Ueda et al. / IFAC-PapersOnLine 49-31 (2016) 49–54
p2
p1 d1 , u1 = 0.24
p3 p0
P(u1 ) d2 , u2 = 0.47 p3 , P(u4 ), d4 , u4 = 1
p4
p0 , P(u0 ), d0 , u0 = 0
P(u2 ) P(u3 )
p5
p9
d3 , u3 = 0.76 Control Point (p) B´ ezier Curve (P)
p8
Discretized Points of the B´ ezier Curve (u)
p6
Data Points (d)
p7
Curve Error
Fig. 1. Sequence of three cubic B´ezier curve. Points p3 and p6 are transition points between B´ezier curves.
Fig. 2. Distance determined by an approximated method. The methods proposed by Hasegawa et al. (2013) and Pandunata and Shamsuddin (2010) could wrongly determine the points as shown in Fig. 2. Maekawa et al. (2007) has some advantages, although its processing cost is higher. On the other hand, the method proposed by Tavares et al. (2011) is the one with more advantages due to the balance between precision and processing time. In this method, the higher the discretization the better the quality and the processing time increases.
that is a composition of three segments of cubic B´ezier curve. 2.2 B´ezier curve fitting The objective of curve fitting problem is to define a curve that approximates a given sequence of points by determining its control points (p0 , ..., pn ). This determination is done by minimizing the discrepancy between the curve and the given sequence of points, which can be calculated as the distance between the curve and each point of the given sequence of points. The function that describes the function to be minimized is f (p0 , ..., pn ) =
m−1 k=1
p2
|dk − P(uk )|2
3. COST FUNCTIONS Equation (5) is the first cost function adopted. This expression evaluate just the discrepancy between curve and given points. Other functions could be evaluated. In this work it will be also evaluated the curve length and the absolute difference between the length of the curve and the given sequence of points. Once (5) defines a family of curves, the other two functions will evaluate the difference among the curves of the family.
(5)
where 0 < k < m , n+1 is the number of control points and d is the sequence of points with m+1 elements. |dk −P(uk )| is the distance between one point dk and the B´ezier curve. uk is the parameter value that defines which point of the curve is the closest to the point dk , which also can be interpreted as the projection of the point dk over the curve. Points d0 and dm are not added in the sum, once the B´ezier curve always interpolates the first and last control points.
The method to determine the discrepancy between point and curve is an enhancement of the method proposed by Tavares et al. (2011). In the determination of the point of a discretized curve that is closer to a given point dk , there are considered two triangles as shown is Fig. 3. Squares points are points of the discretized curve and circle point is a point of the given sequence of points.
There are several methods to determine the parameter uk . Hasegawa et al. (2013) and Pandunata and Shamsuddin (2010) determined this parameter in an approximated approach by using the chord length parametrization. This parametrization approximates the uk parameter as a fraction between the sum of the chords until the respective point dk and the sum of all the chords. Maekawa et al. (2007) used the Newton-Raphson method to determine the parameter uk which makes the segment that connect the points P(uk ) and dk orthogonal to the tangent at P(uk ) point. Tavares et al. (2011) used a method in which the curve is discretized with a sequence of points and all the points are measured to the given point dk and it is determined the point of the discretized curve that is closer to dk .
Firstly, a search is done to determine which point of the discretized curve is closer to the point of the sequence dk . This point is named P(uk ). As shown in Fig. 3, two triangles are defined. The first is defined by the points dk , P(uk ) and P(uk−1 ); the second is defined by dk , P(uk ) and P(uk+1 ). Then it is evaluated which triangle has the smallest area. Next, it is made a verification whether the triangle with smallest area is acute or obtuse, in specific the angles between the segments dk P(uk ) and P(uk )P(uk−1 ) for the first triangle, and dk P(uk ) and P(uk )P(uk+1 ) for the second. If the triangle with the smallest area is acute, it is calculated the height of the triangle related to the segment that is defined by two points of the discretized curve, and 58
2017 IFAC IMS December 5-7, 2016. Austin, TX, USA
Edson Kenji Ueda et al. / IFAC-PapersOnLine 49-31 (2016) 49–54
51
dk
Sequence of Points
P(uk−1 )
P(uk )
P10 P11
...
1st Curve
P1v-1 P1v
...
... j-2 Curves
...
Bézier curve
(b)
Fig. 5. (a) First limit to position of new curve transition; (b) Second limit to position of new curve transition. the dimension of a vector that is the parameter β of the expression 4. This algorithm can also process discrete parameters which is the integer numbers that represents an index associated with points of the sequence that there is a transition between curves.
dm-4 dm-3 dm-2 dm-1 dm
Pj0 Pj1
Control Point
Bézier curve
(a)
P(uk+1 )
Fig. 3. After finding the point P(uk ) of the discretized curve that is closer to the point dk , two triangles are defined. d0 d1 d2 d3 d4
Sequence of Points
Control Point
Using Fig. 1 as an example, that shows a composition of three cubic B´ezier curve, the continuous parameters are:
Pj1v-1 Pjv
(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)
Last Curve
Fig. 4. Local search for the closest point. this height is the distance between curve and point. If the triangle with the smallest area is obtuse, than it is used the other triangle to determine the distance between curve and point. And, if both triangles are obtuse, neither of them are used and it is adopted the distance between dk and P(uk ) as being the distance between discretized curve and point.
Coordinate x of point p1 (p1 x) Coordinate y of point p1 (p1 y) Coordinate x of point p2 (p2 x) Coordinate y of point p2 (p2 y) Continuity factor of point p4 (β4 ) Coordinate x of point p5 (p5 x) Coordinate y of point p5 (p5 y) Continuity factor of point p7 (β7 ) Coordinate x of point p8 (p8 x) Coordinate y of point p8 (p8 y)
and the discrete parameters are two index that are associated with some point of the given sequence, in which there is a transition of curves.
Another improvement is done in the first step, where it is searched the point of the discretized curve that is closer to a point of the sequence. Once the piecewise B´ezier curve is used, each segment approximates only a section of the sequence. Then the search is locally done by searching each point only in one segment as shown in the representation of the Fig. 4. In this example, points d5 and dm−4 are points in which there is a transition of B´ezier curves. So, for points d1 , d2 , d3 and d4 , the search for the closest point is done only in the first curve defined by points P11 , . . . , P1v−1 . And for points dm−4 , dm−3 , dm−2 and dm−1 , the search is done only in the last curve defined by points Pj1 , . . . , Pjv−1 .
3.2 Next Candidate Generation To generate a new candidate by modifying a continuous parameter it is used the method developed by Corana et al. (1987) that is given by xk+1 = xk + v · ∆ri · ei
(6)
where ∆ri is the step size in direction ei and v is an uniform random number in [-1, 1]. To modify a discrete parameter, it is defined a simple heuristic. Initially, the curves division are determined according to the division of the number of points of the sequence and the number of curve segments. So, the selection of which point is modified is not random. To modify an index associated with a transition point, it is generated a random value between [0, 1], and if this value is lower than 0.5, the index is decreased by one unity, otherwise it is increased by one unity. This operation modifies which point of the given sequence of points is going to be used as a transition point between two B´ezier curves. This heuristic is limited by two restrictions:
The determination of both lengths can be easily calculated. The length of the given sequence of points is calculated by the sum of all the segments that is determined by consecutive points. And, once the curves were discretized to determine the discrepancy, the length of the curve is determined by the sum of all the segments of the discretized curve. 3.1 Problem Parametrization It is used AMOSA (Archived MultiObjective Simulated Annealing) (Bandyopadhyay et al., 2008) a multiobjective simulated annealing to solve the problem of curve fitting. This algorithm can process very well the continuous variables, such as the control points of the curves and
• the change of position can not be more than four unities from the initial position. This change could not happen in the increment and decrease of the 59
2017 IFAC IMS 52 December 5-7, 2016. Austin, TX, USA
Edson Kenji Ueda et al. / IFAC-PapersOnLine 49-31 (2016) 49–54
index. This rule is adopted to avoid the creation of candidates that changes too much the initial and final position of one B´ezier curve segment. This limit is shown in Fig. 5(a), in which the blue square is the initial position of joint of two segments, and this position can be modified to any green square. • each parameter must have a minimum distance to the next and previous parameter, this limit is done to avoid the creation of too short curve segments that approximates only few points of the sequence. It is adopted a minimum distance of 4 unities, once cubic curves is used to approximate the sequence. And to approximates few points is a waste of computational time. This limit is shown in Fig. 5(b), in which the blue square is the initial position of the joint of two segments and this position can be modified to any green square. Notice that respecting the first limit, the change could be greater, however due to the fact that the end of the second segment is too close, the position of change is shorter. 1 Curve 2 Curves 3 Curves
55 F2
(a)
50
(b)
Fig. 7. Fitting curve for the parabola example with 1, 2 and 3 curve segments. Each segment is in different colors and its control points are diamond shaped in the respective color and the given sequence of points is the black circles. (a) The curve of lower discrepancy of F1 x F2; (b) The curve of lower discrepancy of F1 x F3.
45 0
2,000 4,000 6,000 F1 (a)
120
F3
0.6 0.4
F2
1 Curve 2 Curves 3 Curves
5 Curves 6 Curves 7 Curves
100 80
0.2 0
0 0
10
20
1,000
2,000
F1
30
(a)
F1 (b)
1 F3
Fig. 6. Pareto front of the parabola example. (a) Curve discrepancy (F1) x fitting curve length (F2); (b) Curve discrepancy (F1) x absolute difference between the length of the fitting curve and the length of the sequence of points (F3).
5 Curves 6 Curves 7 Curves
0.5 0 0
4. RESULTS
20
40
60
80
F1 The developed method was tested in three examples. In each example the curve with desired format is created by discretizing it with 50 points and a noise is added in each coordinate. The chosen formats to be tested are a parabola, a curve with inflection and a curve with auto-intersection. To determine the curve discrepancy each curve segment was discretized with 100 points. The objective functions to be minimized are the curve discrepancy
(b)
Fig. 8. Pareto front of the example of the curve with auto-intersection. (a) Curve discrepancy (F1) x fitting curve length (F2); (b) Curve discrepancy (F1) x absolute difference between the length of the fitting curve and the length of the sequence of points (F3).
60
2017 IFAC IMS December 5-7, 2016. Austin, TX, USA
Edson Kenji Ueda et al. / IFAC-PapersOnLine 49-31 (2016) 49–54
53
both images, it is possible to observe that the developed method successfully approximates the given sequence of points with 1, 2 and 3 curves.
(F1), the fitting curve length (F2) and the absolute difference between the length of the fitting curve and the length of the sequence of points (F3). The AMOSA method was applied to F1 with F2; and to F1 with F3. Each example is tested with three different number of curve segments.
The second example is a curve with auto-intersection and it is used 5, 6 and 7 curve segments to approximate the given sequence of points. Fig. 8 shows the Pareto fronts. In Fig. 8(a), it is possible to observe once again that the increase of the discrepancy decreases the fitting curve length. However, in this example, the higher the number of curves the higher the curve length. Fig. 8(b) shows that there is almost no relation between the discrepancy and the absolute difference of lengths. There is no perceptive change in the discrepancy with a change in the difference of lengths. Fig. 9 shows the curves with lower discrepancy by minimizing F1 and F2 in Fig. 9(a); and minimizing F1 and F3 in Fig. 9(b). 2 Curves 3 Curves 4 Curves
F2
120
(a)
100 80 60 40
0
0.5
1
1.5 ·10
F1 (a)
F3
0.6
2 4
2 Curves 3 Curves 4 Curves
0.4 0.2 0 5
10
15
20
F1 (b) (b)
Fig. 10. Pareto front of the example of the curve with inflection. (a) Curve discrepancy (F1) x fitting curve length (F2); (b) Curve discrepancy (F1) x absolute difference between the length of the fitting curve and the length of the sequence of points (F3).
Fig. 9. Fitting curve for the auto-intersection example with 5, 6 and 7 curves. Each segment is in different colors and its control points are diamond shaped in the respective color and the given sequence of points is the black circles. (a) The curve of lower discrepancy of F1 x F2; (b) The curve of lower discrepancy of F1 x F3.
The third example is a curve with inflection and it is used 2, 3 and 4 curve segments to approximate the given sequence of points. Fig. 10 shows the Pareto fronts. Fig. 10(a) shows the same behavior of the Pareto front of the example of the curve with auto-intersection. Fig. 10(b) shows again that there is almost no relation between the discrepancy and the absolute difference of lengths. However, in this example, the higher the number of curves the lower the minimum value of discrepancy. Fig. 11 shows the curves with lower discrepancy by minimizing F1 and F2 in Fig. 11(a); and minimizing F1 and F3 in Fig. 11(b). In Fig. 11(a) all the curves are successfully created, but all curves in Fig. 11(b) a slight hook is presented at the end of the last piece of curve. These hooks are created because
The first example is a parabola and it is used 1,2 and 3 curve segments to approximate the given sequence of points. Fig. 6 shows the Pareto fronts and it is possible to observe from Fig. 6(a) that the lower the discrepancy the higher the length of the fitting curve. The number of curves has low influence over the discrepancy due to the simplicity of this example. Fig. 6(b) shows that by minimizing the discrepancy you increase the difference of lengths, however after a value there is no change of the difference of lengths when you increase the discrepancy. Fig. 7 shows the curves with lower discrepancy by minimizing F1 and F2 in Fig. 7(a); and minimizing F1 and F3 in Fig. 7(b). In 61
2017 IFAC IMS 54 December 5-7, 2016. Austin, TX, USA
Edson Kenji Ueda et al. / IFAC-PapersOnLine 49-31 (2016) 49–54
5. CONCLUSION
the length of the sequence of points is slightly longer than the desired curve.
An algorithm to solve the curve fitting problem was developed, that used piecewise cubic B´ezier curves with a weak C 1 continuity and a multiobjective simulated annealing. The cost functions analyzed were the fitting curve discrepancy, the fitting curve length and the absolute difference of the fitting curve length and given points length. These developed methods were tested and applied in three different examples and almost all the examples were successfully approximated. As future works, the developed method used several numbers of curves, however it is necessary to study a method to determine the optimum number of curves. ACKNOWLEDGEMENTS MSG Tsuzuki and TC Martins were partially supported by CNPq (grants 310.663/2013-0 and 306.415/2012-7). RY Takimoto and EK Ueda were partially supported by CNPq (grants 150508/2015-8 and 140518/2015-0). AK Sato was supported by CAPES/PNPD. REFERENCES (a)
Adi, D., Shamsuddin, S., and Hashim, S.Z.M. (2010). Nurbs curve approximation using particle swarm optimization. In Computer Graphics, Imaging and Visualization (CGIV), 2010 Seventh International Conference on, 73–79. Bandyopadhyay, S., Saha, S., Maulik, U., and Deb, K. (2008). A simulated annealing-based multiobjective optimization algorithm: AMOSA. Evolutionary Computation, IEEE Transactions on, 12(3), 269–283. Corana, A., Marchesi, M., Martini, C., and Ridella, S. (1987). Minimizing multimodal functions of continuous variables with the simulated annealing algorithm. ACM Transactions on Mathematical Software, 13, 262–280. Gofuku, S., Tamura, S., and Maekawa, T. (2009). Pointtangent/point-normal b-spline curve interpolation by geometric algorithms. Computer-Aided Design, 41, 412– 422. Hasegawa, A.Y., Jr., R.S.U.R., and Tsuzuki, M.S.G. (2013). Differential evolution optimization for bezier curve fitting. In 11th IFAC Workshop on Intelligent Manufacturing Systems, 233–238. Jak´obczak, D.J. (2015). Probabilistic 2D point interpolation and extrapolation via data modeling. Informatica (Slovenia), 39(1), 53–61. Maekawa, T., Matsumoto, Y., and Namik, K. (2007). Interpolation by geometric algorithm. Computer-Aided Design, 39, 313–323. Pandunata, P. and Shamsuddin, S. (2010). Differential evolution optimization for bezier curve fitting. In Computer Graphics, Imaging and Visualization (CGIV), 2010 Seventh International Conference on, 68–72. Piegl, L. and Tiller, W. (1996). Algorithm for approximate nurbs skinning. Computer-Aided Design, 28(9), 699 – 706. Tavares, R., Martins, T., and Tsuzuki, M. (2011). Simulated annealing with adaptive neighborhood: A case study in off-line robot path planning. Expert Systems with Applications, 38(4), 2951 – 2965.
(b)
Fig. 11. Fitting curve for the inflection example with 2, 3 and 4 curve segments. Each segment is in different colors and its control points are diamond shaped in the respective color and the given sequence of points is the black circles. (a) The curve of lower discrepancy of F1 x F2; (b) The curve of lower discrepancy of F1 x F3.
62