Bézier Curve Fitting With A Parallel Differential Evolution Algorithm

11th IFAC Workshop on Intelligent Manufacturing Systems The International Federation of Automatic Control May 22-24, 2013. São Paulo, Brazil

ThCT2.2

B´ ezier Curve Fitting With A Parallel Differential Evolution Algorithm Allan Y. Hasegawa ∗,1 Roberto S. U. Rosso Jr. ∗,2 Marcos Sales Guerra Tsuzuki ∗∗,3 ∗

Universidade do Estado de Santa Catarina, Joinville, Brazil. Computer Science Department ∗∗ Escola Polit´ecnica da Universidade de S˜ ao Paulo, S˜ ao Paulo, Brazil. Mechatronics and Mechanical Systems Engineering Department Computational Geometry Laboratory Abstract: The curve fitting problem is present in several areas, including signal processing, statistics and geometric modeling. Computationally expensive algorithms were proposed to solve this problem, some using artificial intelligence. A previous work implemented a sequential Differential Evolution (DE) algorithm to solve the problem of B´ezier Curve Fitting. This work implements and analyzes a solution to approximate B´ezier curves using a parallel Differential Evolution algorithm. Tests were performed to verify the objective function in several generations of the DE algorithm, as well as an analysis of this solution’s speedup and parallel efficiency. Keywords: curves; parallel; approximation; heuristics; algorithm. 1. INTRODUCTION The curve fitting problem with discrete points as input is an essential task in many engineering problems. Many curve fitting methods are known: spline interpolation, B´ezier piecewise curve interpolation and others (Bartels et al., 1987). Spline interpolation defines a second order continuous curve. The interpolated curve can be represented as a B´ezier piecewise curve in which the number of control points is three times the number of interpolating discrete points. Such fitting methods are used for reverse engineering where the geometry must be preserved. Other applications where human interaction is necessary, a smoothed curve with fewer control points is preferred. An unique B´ezier curve can be better manipulated by a user (Kitson, 1989; Piegl and Tiller, 1997). The curve fitting problem can be approached as an optimization problem in which the summation of the distances between each given discrete point and the fitting curve is minimum. Usually, meta-heuristics are used to solve such optimization problem: genetic algorithms (Afshar et al., 2011), simulated annealing (Tavares et al., 2011; Tsuzuki et al., 2006; Martins et al., 2001) and others. Pandunata and Shamsuddin (2010) proposed a new method to solve the curve fitting problem using the evolutionary algorithm called Differential Evolution (DE), a heuristic algorithm created for minimizing nonlinear and non differentiable continuous space functions. Pandunata and Shamsuddin (2010) used only the sequential variation of the DE algorithm. Meta-heuristics evaluate the objective function several times, and the curve-point distance evaluation is a very 1 2 3

e-mail: [email protected]. e-mail: [email protected]. e-mail: [email protected].

978-3-902823-33-5/13/$20.00 © 2013 IFAC

important procedure. Different methods have been used to solve such problem. Hu and Wallner (2005) proposed a Newton-Raphson based approach where the distance is precisely determined. Tavares et al. (2011) evaluated several points from the fitting curve and determined the nearest point pairs by comparing the two point sets (one from the fitting curve discretization and the input discrete point set). This work uses a faster solution, for each given discrete point, the the B´ezier curve parameter corresponding to the nearest point is approximately determined. This work describes an implementation and analysis, to solve the approximation curve fitting problem for B´ezier curves using the DE algorithm. The method proposed by Pandunata and Shamsuddin (2010) is expanded using a parallel solution coded in C++. Section 2 describes the B´ezier curve. Section 3 presents the curve fitting problem for a B´ezier curve as an optimization problem. Section 4 introduces the DE algorithm. Section 5 presents the method proposed and implemented in this work. Section 6 shows and discusses this work’s results. Finally, Section 7 presents the conclusions. This work’s source code is available at https://github. com/AranHase/pdebc, and released as an open-source project with the Apache License, Version 2.0 (2004). ´ 2. BEZIER CURVES The B´ezier curve approximates a set of control points. These control points define the curve shape, however, the curve may not pass through them. The only guarantee of interpolation happens in the first and last control points, because the curve starts and ends in them. Another important characteristic of the B´ezier curve is the global control, where a change in any control point will modify

233

the entire curve. The parametrization of the curve can be made with the equation n X p(u) = pi Bi,n (u), u ∈ [0, 1] (1) i=0

where the variable pi represents the control points, (n − 1) is the curve degree, see Mortenson (1997). basis functions can be determined by   n i Bi,n (u) = u (1 − u)n−i , i = 0, . . . n i
where ni is the binomial coefficient function:   n n! = i i!(n − i)!

and The (2)

(3)

´ 3. BEZIER CURVE FITTING The B´ezier curve fitting problem tries to reduce the distance between the curve and a set of geometric data points. It is important to know the distance between the curve and the data points. Then, assuming a number n of control points, a set d of data points with m elements, the equation d0 = p(0), dm = p(1), m−1 X f (p1 , · · · , pm−1 ) = kdk − p(u)k2 (4)

Fig. 1. Calculating the distance between a B´ezier Curve, with 4 control points, and 5 data points.

k=1

calculates the curve error, or how far away the curve is from the data points (Pandunata and Shamsuddin, 2010), where p is the B´ezier Curve defined by Equation (1). The values d0 and dm are predefined, because it is known that the B´ezier curve will always start and end on the first and last control points, respectively. The curve error is determined by summing the distance of the data points to the B´ezier curve. Tavares et al. (2011) evaluated the B´ezier Curve at several discrete points and determined the nearest discrete point to each of the inputted points. This approach has drawbacks since the objective function is influenced by the number of points used to discretize the B´ezier curve. A better approach is to evaluate and estimate the parameter associated with the point on the B´ezier Curve that is nearest to a given point. Different methods to create those parametric values will make the algorithm approximate different parts of the curve to the data points, thus changing the curve’s final form. The three most common methods are: equally spaced, chord length and centripetal method. The chord length method is the most widely used (Piegl and Tiller, 1997), and the one used in this work. The chord length method uses Equation (5) u0 = 0 um = 1 |dk − dk−1 | uk = uk−1 + P , k ∈ (1, . . . m) (5) m |dj − dj−1 | j=1

to calculate the distance between each data point, then normalize these values between [0, 1] (Pandunata and Shamsuddin, 2010). For example, to calculate the distance from the B´ezier Curve in Figure 1 to a set of data points (the d variable in 978-3-902823-33-5/13/$20.00 © 2013 IFAC

Initialization - Mutation - Crossover - Selection

Fig. 2. DE algorithm steps (Pandunata and Shamsuddin, 2010). Figure 1), the parameterization values are first generated using Equation (5). In Figure 1, these values are shown with the variable u. Each data point d is related to a single parameterization value u between [0, 1]. The next step is to calculate the points of the curve related by these parameterization values. Therefore, Equation (1) is used with the B´ezier Curve control points, variable c, and a parameterization value u, to find the variable p. Finally, the error of the curve is the sum of the difference between each data point d to its related B´ezier Curve point p. These differences are drawn with dashed lines in Figure 1. The B´ezier curve fitting problem can be reduced to a problem of minimizing the result of Equation (4), while having the position of the B´ezier curve control points as variables. 4. DIFFERENTIAL EVOLUTION Differential Evolution, or DE, is a heuristic algorithm created for minimizing nonlinear and non differentiable continuous space functions (Storn and Price, 1997). This algorithm can be divided into four steps, as in Figure 2. The following subsections describe all four steps of this algorithm, and then finish with the explanation of the parallel model. 4.1 Initialization The first step, initialization, starts a set of N vectors

234

xi,G , i = 0, 1, 2, . . . N − 1 (6) for use as a population for each generation G. Usually, this population is generated using random values (Storn and Price, 1997). 4.2 Mutation The basic principle of the DE algorithm is the method for generating vectors for trial. This process creates new vectors via an operation between two elements of the population, then combines characteristics from a third member. The equation v = xr3,G + F · (xr1,G − xr2,G ) (7) r1, r2, r3 ∈ [0, N − 1], r1, r2, r3 ∈ N0 , r1 6= r2, r1 6= r3, r2 6= r3, F >0 presents this process in a formal way. The variables r1, r2 and r3 are chosen randomly. N is the number of vectors and F is a constant called the mutation factor (Storn and Price, 1997; Pandunata and Shamsuddin, 2010). 4.3 Crossover The crossover, also known as recombination, idea is inserting mutations into the population, but still retaining some original vectors. This process is defined by the equation  v if r ≤ CR u = (x ) if r > CR , CR ∈ [0, 1] (8) i,G where r is a random value between [0, 1]. CR is the crossover probability and v is the variable from Equation (7) (Storn and Price, 1997; Pandunata and Shamsuddin, 2010). 4.4 Selection The last step verifies which vectors will survive for the next generation. For a generation G, the vectors u and xi,G , from Equation (8), are placed for trials. As the DE is essentially a minimization algorithm, the selection of which vector will be part of the generation G + 1 is done using the objective function, where the vector that has the lowest result will be chosen. This process is defined by the equation  u if f (u) ≤ f (xi,G ) xi,G+1 = (9) xi,G if f (u) > f (xi,G ) where f is the objective function of the minimization problem (Storn and Price, 1997; Pandunata and Shamsuddin, 2010). 4.5 Parallel Differential Evolution The DE algorithm can be easily parallelized because each member of the population is analyzed individually. Typically you have two models of parallelism. The first creates a queue where the population is stored, and each individual is processed by one processor, thus resulting in a fine grained parallelism. However, this method can be problematic when there is a limited number of processors, 978-3-902823-33-5/13/$20.00 © 2013 IFAC

Fig. 3. Illustration of the migration process. or when the objective function requires global information about the population (Tasoulis et al., 2004). The second model divides the population into subgroups and then processes each subgroup in parallel. The communication between subgroups happens during a new stage called migration, which allows the best candidate of each subgroup to migrate to a neighboring subgroup. The migration is the process that creates a copy of a vector in one population, then sends it to another population to replace a random vector in the destination. The populations are organized as a circular queue to maintain a migration order. Each group offers its best candidate for migration, however, any subgroup can refuse to receive members from others populations. The chances of that happening are determined by the variable φ ∈ [0, 1], where, at each iteration, a random value between [0, 1] is compared against φ, and if it is inferior, the migration will be accepted (Tasoulis et al., 2004). This model is demonstrated in Figure 3. In this example, the population is divided into three subgroups. At generation G, each thread does a DE step, individually and in parallel, for its local population, thus creating the population of generation G+1. After all threads are at generation G + 1, the migration step starts. The best candidate from each subgroup migrates to the next subgroup. This process is repeated until all generations have been created. 5. THE PROPOSED METHOD The B´ezier curve fitting problem can be solved with the DE algorithm following the same steps as in Figure 2. By using Equation (4), the curve fitting problem becomes a minimizing problem, where the objective is to reduce the distance between the curve and the data points. The B´ezier curve’s control points become the variables of the system. Each control point of the B´ezier curve is processed individually and sequentially. The first and last control points are placed at the position of the first and last data point, respectively, and both are ignored until the end of the algorithm. The remaining control points are then approximated using the DE algorithm. The steps of initialization, mutation and crossover are processed as described in Section 4. The initial population created by the initialization step, Section 4.1, is generated using random values. In the selection step, in Equation (9), the function f is replaced by Equation (4). This step is the

235

foreach DE Generation do foreach Curve Control Point do foreach Vector in Population do for k < number of Data Points do // Equation (5) m ← number of Data Points d ← array of Data Points |dk −dk−1 | uk ← uk−1 + P m |dj −dj−1 |

j=1

// Equation (1) p ← array of Curve Control Points n ← number Control Points −1 for i < number Control Points do n! t ← t + pi ∗ i!(n−i)! ∗ uik ∗ (1 − uk )n−i // Equation (4) curves error ← curves error + kdk − tk2

Algorithm 1. Naive implementation to calculate the curve error.

// Cache Equation (5) m ← number of Data Points d ← array of Data Points for k < number of Data Points do |dk −dk−1 | uk ← uk−1 + P m |dj −dj−1 |

j=1

B ← Cache the Basis Functions foreach DE Generation do foreach cp in Curve Control Point do // Calculate Equation (1) partially for k < number of Data Points do foreach cpa in Curve Control Point do if cpa 6= cp then ptk ← ptk + cpa ∗ Bk foreach v in Population do // Equation (4) for k < number of Data Points do p ← ptk + v ∗ Bk curves error ← curves error + kdk − pk2

Algorithm 2. Implementation used in this work to calculate the curve error.

one that has the highest consumption of computational power because the algorithm needs to recalculate the curve’s error, with Equation (4), every time the mutation step modifies an element of the population. To make matters worse, Equation (4) needs the variable u, generated by Equation (5), to expand the function p(u) with Equation (1). The pseudo code in Algorithm 1 demonstrates this complexity, showing how line 13, relating to Equation (4), is processed G ∗ P ∗ C 2 ∗ D times, where G is the number of generations, P is the population size, C is the number of control points of the B´ezier Curve and D is the number of Data Points. To reduce the complexity of the selection step, the variables uk from Equation (5) can be generated only once, at 978-3-902823-33-5/13/$20.00 © 2013 IFAC

spawn N threads // N is the number of processes foreach cp in CP do // CP is the set of the curve’s // Control Points initialize a population for cp subdivide cp’s population into N subgroups assign the right population’s subgroup to each thread while stop condition not yet met do foreach cp in CP do pre-compute equations foreach thread do // this scope is done in parallel perform a DE step on cp’s population M ← best candidate for migration // Migration step foreach m in M do r ← random number ∈ [0, 1] if r ≤ φ then migrate m to the next process population Algorithm 3. High level description of this proposed method.

the start of the algorithm, and stored in memory because there are no changes to the data points during the curve fitting process. Equation (1) also does not need a complete recalculation every time the selection step is invoked. Taking advantage of the characteristic that the data points are constant, the basis functions Bi,n (u), from Equation (1), can be precomputed. Another characteristic that can be used to reduce complexity is that each control point is executed in sequence. While a control point is being evaluated, the others do not change, as well as the basis functions and data points, allowing the pre-computation of the Equation (1) partially. Algorithm 2 demonstrates this optimization in a pseudo code, showing how line 17, relating to Equation (4), is processed G ∗ C ∗ P ∗ D times, where G is the number of generations, C is the number of control points in the B´ezier Curve, P is the population size and D is the number of Data Points. The parallelism model adopted is that of division of the population into subgroups and migration, as described in Section 4.5. The way in which this division occurs determines the load balancing algorithm. However, as the processing time for different vectors of the population tends not to have significant variations, in this work the vectors are evenly distributed for each process. The message exchange between processes is minimal, and happens only in the stage of migration. Algorithm 3 presents a high level description of this proposed method. 6. TESTS AND RESULTS This work was implemented in the programming language C++. The software environment used for testing was: Linux Mint 14 64bits running on the console with a bare minimum of background services. The hardware configuration used was: an Intel i7-2630m CPU with HyperThreading and 6Gbytes of DDR3-1333Mhz RAM. The application was compiled with the optimization flags “-

236

Fig. 4. Fitting of a curve with 25 data points and 20 control points after 500 generations. The circles connected by red lines are the data points, while the blue path is the B´ezier Curve. Fig. 6. Processing time for curve fitting when varying the number of control points.

Fig. 5. Behavior of the objective function in different generations.

Fig. 7. Speedup when fitting a curve with 14 control points to 14 data points.

O3 -march=native -fomit-frame-pointer -pipe” with the Clang 3.3 compiler, see clang (2013), and the llvm libc++ Revision 177372, see libc++ (2013). The tests shown in this section were performed using the value 0.8 for the chance of mutation, CR, and 0.5 for the mutation factor, F . A population of size 128 was used for each B´ezier Curve control point. The stop condition is the completion of 500 generations. The value of these constants reflect the choices made in the Pandunata and Shamsuddin’s (2010) work. The data points used are shown in Figure 4. The first test checks the behavior of the objective function when approximating a curve with 20 control points to 25 data points, as in Figure 4. The results of the objective function, Equation (4), for 500 generations can be seen in Figure 5, where it can be seen that the curve makes most of its changes by the 64th generation, and then slowly improves until the last generation. The next test analyzes the processing time for approximating a curve for 500 generations. Figure 6 shows the results for a linear increase in the number of control points. These times were collected by running the parallel algorithm with 4 processes. By reducing the complexity of Equation (1), the growth of the processing time becomes linear when increasing the number of control points. The last analysis tests the parallel version of this work. Figure 7 shows the speedup of the algorithm when fitting a curve with 14 control points to 14 data points. Up to four 978-3-902823-33-5/13/$20.00 © 2013 IFAC

Fig. 8. Efficiency when fitting a curve with 14 control points to 14 data points. processes there is an increase in speedup, however, from the eighth process, this improvement begins to decline. Figure 8 shows the efficiency of the parallel algorithm. 7. CONCLUSIONS This work implemented an algorithm to solve the B´ezier curve fitting problem by using the parallel DE algorithm. The obtained B´ezier curve can be interactively manipulated. On a PC with four physical and eight logical cores, the speedup result proved that this solution accelerates the algorithm when using up to eight processes when compared with the sequential solution.

237

In future works, the authors intend to improve the efficiency of the parallel algorithm. Also, the length of the curve will be considered in future works, in which the fitting curve should have minimum distance to the given points and smaller lenght. ACKNOWLEDGEMENTS The authors RSU Rosso Jr. and AY Hasegawa wish to thank the support from the Santa Catarina State University. This work was supported in part by the National Technological Agency. MSG Tsuzuki was partially supported by the CNPq (grant 309570/2010-7). REFERENCES Afshar, N., Soryani, M., and Rahmani, A. (2011). Curve fitting using coevolutionary genetic algorithms. In SEMCCO (2), 201–210. Apache (2004). Apache license, version 2.0. URL http:// www.apache.org/licenses/LICENSE-2.0.html. [access 01-05-2013]. Bartels, R., Beatty, J., and Barsky, B. (1987). An Introduction to Splines for Use in Computer Graphics And Geometric Modeling. Morgan Kaufmann. Clang (2013). clang: a c language family frontend for llvm. URL http://clang.llvm.org/. [access 01-11-2013]. Hu, S.M. and Wallner, J. (2005). A second order algorithm for orthogonal projection onto curves and surfaces. Computer Aided Geometric Design, 22(3), 251 – 260. Kitson, F. (1989). An algorithm for curve and surface fitting using b-splines. In International Conference on Acoustics, Speech, and Signal Processing, 1207 –1210. libc++ (2013). “libc+” c++ standard library. URL http://libcxx.llvm.org/. [access 01-11-2013]. Martins, T., Moraes, V., and Tsuzuki, M. (2001). Simulated annealing applied to path planning. In 2nd IFAC Workshop on Intelligent Assembly and Disassembly. Mortenson, M.E. (1997). Geometric modeling (2nd ed.). John Wiley & Sons, Inc., New York, NY, USA. 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. (1997). The Nurbs Book. Monographs in Visual Communication Series. Springer-Verlag GmbH. Storn, R. and Price, K. (1997). Differential evolution a simple and efficient heuristic for global optimization over continuous spaces. J. of Global Optimization, 11(4), 341–359. Tasoulis, D., Pavlidis, N., Plagianakos, V., and Vrahatis, M. (2004). Parallel differential evolution. In Evolutionary Computation, 2004. CEC2004. Congress on, volume 2, 2023 – 2029 Vol.2. Tavares, R.S., Martins, T.C., and Tsuzuki, M.S.G. (2011). Simulated annealing with adaptive neighborhood: A case study in off-line robot path planning. Expert Systems and Applications, 38(4), 2951–2965. Tsuzuki, M.S.G., Martins, T.C., and Takase, F.K. (2006). Robot path planning using simulated annealing. In 12th IFAC/IFIP/IFORS/IEEE/IMS Symposium Information Control Problems in Manufacturing, 175–180. 978-3-902823-33-5/13/$20.00 © 2013 IFAC

238