Multilayer embedded bat algorithm for B-spline curve reconstruction

Abstract

This paper presents a new method called multilayer embedded bat algorithm (ME-BAT) to solve the general curve reconstruction problem with free-form parametric B-splines. Opposed to previous approaches in the literature, this method computes the optimal values of all free variables (data parameters, breakpoints, and poles), a very difficult task because they are strongly intertwined in a highly nonlinear way. The method is based on the idea of applying the bat algorithm at different layers: a main bat algorithm at an upper layer to compute the breakpoints and a second bat algorithm at a lower layer to compute the data parameters. This second bat algorithm is embedded into the first one and executed for each breakpoint vector of the population and at each iteration step of the main algorithm. Then, the poles are calculated by least-squares minimization through SVD. The method has been applied to three real-world engineering examples. The experimental results show that the method performs very well, being able to recover the underlying shape of data with high accuracy. A comparison with eleven alternative methods (including six classical methods in the field and all the metaheuristic methods applied so far to this problem) shows that this method outperforms the previous approaches in the field for all instances in the benchmark.

Keywords

CAD/CAM reverse engineering curve reconstruction B-spline curves bio-inspired computation bat algorithm

1. Introduction

A striking feature of current manufacturing processes is the increasing trend of mass customization. Nowadays, products are manufactured in lesser amounts but with greater product diversity. A critical issue is how to use flexible computer-aided engineering systems to optimize mass customization without the corresponding increase in production costs [19]. In this context, design is becoming a major factor in the product development lifecycle. Individual customization implies that the geometric model of a product has to be frequently modified during the design process. As a result, computer-aided design and manufacturing (CAD/CAM) systems are becoming indispensable tools for industrial production and manufacturing.

A common strategy for mass customization in many industrial fields is to build prototypes to guide the early conceptual design. Typically, the process begins with a volume of raw material (such as clay, foam, wood, or plastics) from which a physical prototype is created. This physical model is digitally stored through data sampling by using technologies such as 3D laser scanning. The output is a cloud of data points. Then, curve and surface reconstruction techniques are applied to recover the underlying shape of data. This process, called reverse engineering, has many advantages: the digital model is much easier to understand and analyze than the physical model. It is also cheaper to customize, since any modification is simulated by computer and can readily be reversed if needed. Finally, digital models can be efficiently stored, archived and transferred among manufacturers, providers and customers, using the current high-speed networks and enhanced methods for feature-based data exchange and handling and big data cloud manufacturing [3, 34, 45].

Owing to these reasons, reverse engineering is the prevalent technology in many industrial fields such as CAD/CAM and computer-numerically controlled (CNC) machining and milling for automotive, aerospace and ship building industries. Another popular example is rapid prototyping (the generation of scale models of physical objects from CAD data), mostly due to the popularity and affordable prices of 3D printers, additive layer manufacturing technologies, and handheld scanning devices. Other examples include pressure machines in shoemaking industry, coordinate measuring machines for metrology and quality assessment, and computer tomography and magnetic resonance for medicine and biomedical engineering (e.g., prosthesis and customized medical implants).

In this context, this paper presents a method called multilayer embedded bat algorithm (ME-BAT) to solve the general curve reconstruction problem with parametric B-splines. The method relies on the bat algorithm, a powerful metaheuristic technique for continuous optimization problems. Opposed to all previous approaches, our method computes the optimal values of all free variables at full extent. The experimental results show that it performs very well, being able to recover the underlying shape of data with high accuracy. Comparison with eleven alternative methods shows that this method outperforms previous approaches in the field for all instances in the benchmark.

The structure of this paper is as follows: Section 2 describes the B-spline curve reconstruction problem and the state of the art in the field. Section 3 describes our proposed method in detail. The experimental results for three illustrative real-world examples are reported in Section 4. A comparison of this method with other alternative approaches is discussed in Section 5. The paper closes in Section 6 with the main conclusions and some plans for future work.

2. Related work

The problem of curve reconstruction from data points can be stated as follows:

Given a (generally large) collection of (usually noisy) data points, obtain the curve providing the best fitting to the input data points, according to an error metrics.

When data is acquired from accurate sources and smooth shapes, it can be performed through popular interpolation techniques. However, real-world data are usually affected by measurement noise, so interpolation methods fail since they force the fitting functions to pass through the noisy outliers. To overcome this drawback, approximation techniques are applied instead. In this case, curve reconstruction can be mathematically formulated as an optimization problem, generally a least-squares fitting problem [6, 26, 27].

Solving this optimization problem requires an adequate choice of the fitting functions. Arguably, the best ones are the free-form parametric curves, which are very flexible and well-suited for interactive design. Among them, the B-spline curves are the most powerful fitting curves, and are now included as basic primitives in all major CAD/CAM programs of the market. Based on these reasons, in this paper we address the curve reconstruction problem with B-splines curves, which are briefly described in next paragraphs. The interested reader is referred to [27] for further details.

Let $\mathbf{\Upsilon}=\{\upsilon_{0}=\alpha,\upsilon_{1},\dots,\upsilon_{\nu},% \upsilon_{\nu+1}=\beta\}$ be a strictly increasing sequence of non-negative real numbers called breakpoints on the compact interval $[\alpha,\beta]\subset\mathbb{R}$ . Note that in this paper vectors are denoted in bold. Without loss of generality, $[\alpha,\beta]$ can be assumed to be $[0,1]$ . For each sequence $\mathbf{\Upsilon}$ , the $k$ -th B-spline basis function $\psi_{k}^{\rho}(\omega)$ of order $\rho$ (or equivalently degree $\rho-1$ ) is defined recurrently as:

$\psi_{k}^{1}(\omega,\mathbf{\Upsilon})=\begin{cases}1&{\rm if}\,\upsilon_{k}% \leqslant\omega<\upsilon_{k+1}\\ 0&{\rm otherwise}\\ \end{cases}$ (1)

for $k=0,\dots,\nu$ , and:

$\displaystyle\psi_{k}^{\rho}(\omega,\mathbf{\Upsilon})={{\omega-\upsilon_{k}}% \over{\upsilon_{k+\rho-1}-\upsilon_{k}}}\psi_{k}^{\rho-1}(\omega,\mathbf{% \Upsilon})$ $\displaystyle\quad∼{}+{{\upsilon_{k+\rho}-\omega}\over{\upsilon_{k+\rho}-% \upsilon_{k}}}\psi_{k+1}^{\rho-1}(\omega,\mathbf{\Upsilon})$ (2)

for $\rho>1$ and $k=0,\dots,\nu-\rho+1$ . If necessary, the convention ${{{0}\over{0}}}=0$ in Eq. (2) is applied. Note that $\psi_{k}^{\rho}(\omega,\mathbf{\Upsilon})$ is a nonlinear function of $\mathbf{\Upsilon}$ .

A B-spline curve of order $\rho$ and breakpoints $\mathbf{\Upsilon}$ , denoted as $\mathbf{\Phi}(\omega,\mathbf{\Upsilon})$ , is given by:

$\mathbf{\Phi}(\omega,\mathbf{\Upsilon})={\displaystyle\sum\limits_{i=0}^{% \sigma}\mathbf{\Xi}_{i}\,\psi_{i}^{\rho}(\omega,\mathbf{\Upsilon})}$ (3)

where ${\mathbf{\Xi}}=\{\mathbf{\Xi}_{i}\}_{i=0,\dots,\sigma}$ are vector coefficients of the curve called poles and $\left\{\psi_{i}^{\rho}(\omega,\mathbf{\Upsilon})\right\}_{i}$ are the basis functions defined above. We impose the boundary conditions $\upsilon_{-\rho+1}=\upsilon_{-\rho+2}=\dots=\upsilon_{-2}=\upsilon_{-1}=% \upsilon_{0}=\alpha$ and $\upsilon_{\nu+1}=\upsilon_{\nu+2}=\dots=\upsilon_{\nu+\rho}=\beta$ , with the effect to force the curve to interpolate the first and last poles, i.e., $\mathbf{\Phi}(\alpha,\mathbf{\Upsilon})=\mathbf{\Xi}_{0}$ and $\mathbf{\Phi}(\beta,\mathbf{\Upsilon})=\mathbf{\Xi}_{\sigma}$ , a very desirable property in industrial applications of computer design and manufacturing and many other fields. Note however that the method described in this paper is actually independent of these constraints.

Let us now consider an input set of measured data points $\{\mathbf{\Theta}_{k}\}_{k=1,\dots,\chi}$ , obtained by 3D laser scanning, coordinate measuring machines, industrial computer tomography, or other digitizing methods, as it happens in many industrial and engineering problems. The goal of B-spline curve reconstruction consists of obtaining the B-spline curve $\mathbf{\Phi}(\omega,\mathbf{\Upsilon})$ of order $\rho$ approximating better the data points $\{\mathbf{\Theta}_{k}\}_{k}$ . Since the fitting curve is parametric, a suitable parameterization of data points is needed, so that each data point $\mathbf{\Theta}_{k}$ is associated with a data parameter value, denoted as $\xi_{k}\in[\alpha,\beta]$ . Due to the boundary conditions introduced above, we can take $\mathbf{\Phi}(\xi_{1},\mathbf{\Upsilon})=\mathbf{\Theta}_{1}$ and $\mathbf{\Phi}(\xi_{\chi},\mathbf{\Upsilon})=\mathbf{\Theta}_{\chi}$ and perform approximation on the remaining parameters:

$\mathbf{\Theta}_{k}\approx\mathbf{\Phi}(\xi_{k},\mathbf{\Upsilon})={% \displaystyle\sum\limits_{i=0}^{\sigma}\mathbf{\Xi}_{i}\,\psi_{i}^{\rho}(\xi_{% k},\mathbf{\Upsilon})}$ (4)

for $k=2,\dots,\chi-1$ . Equation (4) can be written in matrix notation as:

${\mathbf{\Theta}}={\bf\Psi}.{\bf\Xi}$ (5)

where ${\mathbf{\Theta}}=({\mathbf{\Theta}}_{2},\dots,{\mathbf{\Theta}}_{\chi-1})^{T}$ , ${\bf\Xi}=({\bf\Xi}_{0},\dots,{\bf\Xi}_{\sigma})^{T}$ , ${\bf\Psi}=(\{\psi_{i}^{\rho}(\xi_{l},\mathbf{\Upsilon})\}_{i=0,\dots,\sigma\,;% \,l=2,\dots,\chi-1})$ is the matrix of sampled B-spline basis functions, and $(.)^{T}$ represents the transpose of a vector or matrix. The dimension of the search space $D\subset\mathbb{R}^{n}$ in Eq. (5) is $n(\sigma+1)+\nu+\chi-2$ , which could be very large (from thousands to millions) for very complicated shapes.

System Eq. (5) is generally overdetermined. Consequently, the matrix of basis functions ${\bf\Psi}$ is not invertible and no direct solution can be obtained. In this case, consider the least squares approximation of Eq. (4), defined as the minimization problem given by:

$\Lambda=\underset{\overset{\{\mathbf{\Xi}_{i}\}_{i}}{\overset{\{\upsilon_{j}\}% _{j}}{\{\xi_{k}\}_{k}}}}{\text{minimize}}\left(\sum\limits_{k=2}^{\chi-1}\Bigg% {\|}{\mathbf{\Theta}}_{k}{-}{\sum\limits_{i=0}^{\sigma}\mathbf{\Xi}_{i}\psi_{i% }^{\rho}(\xi_{k},\mathbf{\Upsilon})}\Bigg{\|}_{2}\right)\!\!$ (6)

where ${\|.\|_{2}}$ represents the Euclidean norm. Since the input consists only of the data points, solving Eq. (6) requires to compute all free variables of the problem. This is a very difficult and challenging task, because they are related to each other in a highly nonlinear and complicated way. This explains why previous methods in the literature make a “reasonable” choice of values for some of those variables, which are then used to compute the remaining ones. This procedure is strongly dependent on the choice of values for some variables, a task that requires a lot of expertise. As a result, this process remains manual and no automatic procedure can be applied. Moreover, the problem is nonconvex (as the curve is a nonconvex function of the breakpoints) and multimodal (as there could be several local and/or global minima of the objective function). Therefore, Eq. (6) represents a strongly nonconvex nonlinear high-dimensional multimodal continuous optimization problem.

2.1 Previous work

Table 1
Metaheuristic methods applied to B-spline curve reconstruction.

Method & reference	Year	Disc./cont.	Parametric	$\{\xi_{k}\}_{k}$	$\{\upsilon_{j}\}_{j}$	Additional comments
Genetic algorithms [41]	1999	Discrete	$\times$	$\times$	✓	Only explicit curves supported
Genetic algorithms [42]	2003	Continuous	$\times$	$\times$	✓	Only explicit curves supported
Artificial immune systems [30]	2009	Discrete	✓	$\times$	✓	Large discretization errors
Particle swarm optimization [29]	2010	Continuous	✓	✓	$\times$	Averaging method for breakpoints
Particle swarm optimization [7]	2011	Continuous	$\times$	$\times$	✓	Only explicit curves supported
Gaussian mixture model [46]	2011	Continuous	✓	$\times$	✓	Only closed curves supported
Firefly algorithm [8]	2013	Continuous	✓	✓	$\times$	De Boor method for breakpoints
PESA [31]	2013	Discrete	✓	$\times$	✓	Large discretization errors
Artificial immune systems [10]	2015	Continuous	$\times$	$\times$	✓	Only explicit curves supported
Memetic firefly algorithm [9]	2016	Continuous	$\times$	$\times$	✓	Only explicit curves supported
ME-BAT algorithm (this work)	2017	Continuous	✓	✓	✓	Fully general method

Curve reconstruction with free-form parametric curves has been an active field of research for a very long time. First approaches in the field were based on classical mathematical optimization methods [6]. Early approaches were developed in the 60 s and 70 s for Bézier curves, where the bottleneck is the parameterization step. Due to its simplicity, the uniform parameterization (where all data parameters are equally spaced) is a popular choice. However, in many practical settings (e.g., CNC machining and milling for manufacturing, metrology and quality assessment) it is more convenient to consider arc-length parameterization. This method was later extended to the centripetal method, a variation of the arc-length parameterization that yields better results for shapes with sharp turns [27].

It was soon noticed however that Bézier curves are affected by the so-called global control: moving a pole automatically modifies the entire curve, making interactive design very tedious and difficult to handle. Due to these reasons, Bézier curves have been replaced in industrial applications by the more powerful B-spline curves. Many B-spline curve reconstruction approaches focus on the problem of data parameterization and skip the (more difficult) problem of computing the breakpoints. In general, these approaches proceed by setting the number of breakpoints a priori and then computing their location according to some formula [16, 24, 25, 32]. The simplest way to do it is to consider equally spaced values (uniform breakpoint vector). However, this approach may lead to singular systems of equations and does not reflect the real distribution of data. A more refined procedure consists of the averaging method and its variations [27], which allocate breakpoints to ensure that every interval span contains at least one parametric value.

It has been shown that B-spline curve reconstruction improves dramatically if the breakpoints are treated as free variables of the problem [2, 17]. Typically this is carried out by changing their number through either breakpoint insertion or removal. These methods require terms or parameters (tolerance errors, smoothing factors, cost functions, initial breakpoint locations) whose values are chosen subjectively [6, 11, 17]. Therefore, they fail to automatically generate a good breakpoint vector. In addition, these methods tend to be prone to errors and time-consuming, and yield an unnecessarily large number of breakpoints. Other approaches use error bounds [24], curvature-based squared distance minimization [32], or dominant points [25]. In general, they perform well but require some strong constraints (such as high differentiability, closed curves, noiseless data) that are difficult to meet for many real-world applications. Other methods use curvature information extracted from input data [14, 18, 21], so they are restricted to smooth data points and strongly affected by the noise in data.

In the last few years, different artificial intelligence techniques have been applied to this problem. Some examples are neural networks [12], Kohonen’s SOM (Self-Organizing Maps) nets [13], functional networks [15], and support vector machines [16]. All these methods are restricted to particular (very simple) cases and do not address the general case.

A very promising line of research is given by the bio-inspired metaheuristic techniques, which have been intensively applied to solve difficult optimization problems [1, 4, 5, 20, 22, 33, 39, 49]. Table 1 shows all the metaheuristic methods for B-spline curve reconstruction reported in the literature (arranged in rows in chronological order). For each method, the table shows (in columns) the method used and its bibliographic reference, the year it was introduced, and its discrete/continuous nature. Then, three true/false questions follow to enrich the discussion: parametric B-spline curves supported; data parameters $\{\xi_{k}\}_{k}$ computed; and breakpoints $\{\upsilon_{j}\}_{j}$ computed. Answer true is marked with a check (✓); otherwise with symbol $\times$ . Last column provides additional comments. Roughly speaking, these methods can be classified into two groups: discrete approaches [30, 31, 41] and continuous approaches [7, 8, 9, 10, 29, 42, 46]. Methods in the former group convert the original continuous problem into a discrete combinatorial optimization problem solved by either artificial immune systems [30] or genetic algorithms [31, 41]. As expected, this conversion process introduces large discretization errors, making them both inaccurate and unreliable for real-world problems. The continuous methods avoid the discretization errors but they generally restricted to the simpler case of explicit functions (that do not require data parameterization) and cannot deal with the parametric case [7, 9, 10, 42]. The methods in [8, 29, 46] consider the case of parametric B-spline curves but only compute some of the variables, either the data parameters [8, 29] or the breakpoints [46], so they cannot solve the general case. To summarize, no single method described in the literature so far is able to solve the curve reconstruction problem with parametric B-spline curves in all its generality. As shown in last row of Table 1, our method is aimed at filling this gap.

3. The proposed method

3.1 The bat algorithm

The bat algorithm is a powerful bio-inspired metaheuristic algorithm for optimization based on the echolocation behavior of bats and originally proposed in 2010 by Yang [35, 36]. Despite the short time since its appearance, it has already been applied to several engineering and industrial problems. See also [37] for a comprehensive, updated review of the bat algorithm and all its variants and applications.

This paper considers the standard bat algorithm (as described in the original paper in [35]), whose pseudo-code is shown in Algorithm 1. Basically, the algorithm considers an initial population of $\mathcal{P}$ individuals (bats) representing potential solutions of the optimization problem. Each bat has a location ${\bf x}_{i}$ and velocity ${\bf v}_{i}$ , initialized with random values. Also a pulse frequency, pulse rate, and loudness, which are computed for each individual bat (lines 2–6). Then, the swarm evolves over generations (line 7) until the maximum number of generations, $\mathcal{G}_{\max}$ , is reached (line 22). For each generation $g$ and each bat (line 8), new frequency, location and velocity are computed (lines 9 and 10) as:

$\displaystyle f_{i}^{g}=f_{\min}^{g}+\beta(f_{\max}^{g}-f_{\min}^{g})$ (7) $\displaystyle{\bf v}_{i}^{g}={\bf v}_{i}^{g-1}+[{\bf x}_{i}^{g-1}-{\bf x^{*}}]% \,f_{i}^{g}$ (8) $\displaystyle{\bf x}_{i}^{g}={\bf x}_{i}^{g-1}+{\bf v}_{i}^{g}$ (9)

where $\beta\in[0,1]$ follows the random uniform distribution, and ${\bf x^{*}}$ represents the current global best location (solution). The superscript $(.)^{g}$ is used to denote the current generation $g$ . Then, a local solution around ${\bf x^{*}}$ is selected (lines 11–14). The search is intensified by a local random walk (line 15): once a solution is selected among the current best solutions, it is perturbed locally through a random walk of the form: ${\bf x}_{new}={\bf x}_{old}+\epsilon\mathcal{A}^{g}$ where $\epsilon$ is a uniform random number on the interval $[-1,1]$ and $\mathcal{A}^{g}$ is the average loudness of all the bats at generation $g$ .

Algorithm 1

Bat Algorithm

(Initial Parameters)Population size: $\mathcal{P}$ Maximum number of generations: $\mathcal{G}_{\max}$ Loudness and pulse rate: $\mathcal{A}$ , $r$ Maximum frequency: $f_{\max}$ Objective function: $\phi(\bf{x})$ , with ${\bf x}=(x_{1},\dots,x_{d})^{T}$ Random vectors: $\Theta=(\theta_{1},\dots,\theta_{\mathcal{P}})$ $\Psi=(\psi_{1},\dots,\psi_{\mathcal{P}})\text{ with }\theta_{k},\psi_{k}\in U(% 0,1)$

$g\leftarrow 0$ // $g$ : generation index

for $i=1$ to $\mathcal{P}$ do

Initialize the location and velocity ${\bf x}_{i}$ and ${\bf v}_{i}$

Define pulse frequency $f_{i}$ at ${\bf x}_{i}$

Initialize pulse rates $r_{i}$ and loudness $\mathcal{A}_{i}$

end for

while $g\leqslant\mathcal{G}_{\max}$ do

for $i=1$ to $\mathcal{P}$ do

Generate new solutions by adjusting frequency,

10:

and updating locations and velocities

11:

if $\theta_{i}>r_{i}$ then

12:

$\mathbb{s}^{\textit{best}}\leftarrow\mathbb{s}^{g}$ //select current best global solution

13:

$\mathbb{ls}^{\textit{best}}\leftarrow\mathbb{ls}^{g}$ //generate a local solution around $\mathbb{s}^{\textit{best}}$

14:

end if

15:

Generate a new solution by local random walk

16:

if $\psi_{i}<\mathcal{A}_{i}$ and $\phi(\bf{x}_{i})<\phi(\bf{x^{*}})$ then

17:

Accept new solutions

18:

Increase $r_{i}$ and decrease $\mathcal{A}_{i}$

19:

end if

20:

end for

21:

$g\leftarrow g+1$

22:

end while

23:

Rank the bats and find current best $\bf{x^{*}}$

24:

return $\bf{x^{*}}$

If the new solution achieved is better than the previous best one, it is probabilistically accepted depending on the value of the loudness. In that case, the algorithm increases the pulse rate and decreases the loudness (lines 16–19). The evolution rules for loudness and pulse rate are as: $\mathcal{A}_{i}^{g+1}=\alpha\mathcal{A}_{i}^{g}$ and $r_{i}^{g+1}=r_{i}^{0}[1-\exp(-\gamma g)]$ respectively, where $\alpha$ and $\gamma$ are constants. This process is repeated for the given number of generations.

3.2 Overview of the method

The method introduced in this paper applies the bat algorithm to solve the difficult problem given by Eq. (6). This requires to optimize all variables of the problem, instead of some of them as previous methods do. This cannot be done by a simple application of the bat algorithm; a more sophisticated approach is needed.

The method proposed addresses this problem by applying the bat algorithm at two different levels (called layers). The first (upper) layer computes the optimal vector of breakpoints. To this aim, the method is initialized with a random population of breakpoint vectors. For each, a second population of parameter vectors is randomly generated and then optimized through bat algorithm to obtain the best parameter vector associated with such a breakpoint vector. After this initialization stage, bat algorithm is applied to the population of breakpoint vectors for a given number of iterations. At each iteration step, each individual (the breakpoint vector) is modified according to Eqs (7)–(9), then, the change is probabilistically accepted according to the improvement of the fitness function and the value of the loudness. However, the fitness function depends not only on the breakpoint vector but also on the parameter vector. Therefore, any update of the breakpoint vector means that a new optimal parameter vector is to be computed. This task is carried out at a second (lower) layer, consisting of another bat algorithm embedded into the first (main) one. In other words, for each iteration step of the bat algorithm for the breakpoint vector population, another bat algorithm for a new parameter vector population is executed. The best parameter vector is selected to evaluate the fitness and the main bat algorithm continues until convergence. To reflect this nested structure of layers where the second bat algorithm is embedded into the first one at each iteration step, this procedure is called multilayer embedded bat algorithm (ME-BAT for short).

The output of this ME-BAT algorithm is the optimal pair of parameter and breakpoint vectors. They are then used to compute the optimal values for the poles. The resulting system of equations is overdetermined, so they are computed by using the least-squares method through either LU decomposition or SVD. This final procedure yields the B-spline curve that fits the data points better in the least-squares sense.

3.3 Description of the method

In this section, the method is formalized and each step of the method is discussed in detail.

3.3.1 Input data

The input of the method consists of:

the order of the B-spline curve, $\rho$ ,

the input data, $\{\mathbf{\Theta}_{k}\}_{k=1,\dots,\chi}$ , and

the number of poles, $\sigma$ .

The issue of finding suitable values for the parameters of the method will be discussed in Section 3.4.

3.3.2 Breakpoint placement – data parameterization

As indicated above, these two tasks are dependent of each other and must be conducted simultaneously. They are carried out by applying the new approach introduced in this paper, ME-BAT. It consists of two stages: initialization and main procedure.

Stage 1: Initialization. Firstly, the method is initialized with a random population of breakpoint vectors:

$\bar{\mathbf{\Upsilon}}^{(0)}=(\mathbf{\Upsilon}_{1}^{(0)},\dots,\mathbf{% \Upsilon}_{\eta}^{(0)})$ (10)

where the right superscript ${(.)}^{(0)}$ is used onwards to indicate the number of iterations of the bat algorithm for the optimization of the breakpoint vectors, and $\eta$ is the breakpoint population size. The bar above any vector $\bar{(.)}$ is used to indicate a population of vectors, as in the left-hand side of Eq. (10). For each breakpoint vector $\mathbf{\Upsilon}_{j}^{(0)}$ , $(j=1,\dots,\eta)$ , a population of $\lambda$ parameter vectors denoted as:

$\bar{\mathbf{\Delta}}_{j}^{(0)}=(\mathbf{\Delta}_{j,1}^{(0)},\dots,\mathbf{% \Delta}_{j,\lambda}^{(0)})$

is randomly generated. Then, an initial optimization on this population of parameter vectors is performed through a bat algorithm. This procedure is denoted as:

$\mathscr{B}_{j}(\{\bar{\mathbf{\Delta}}_{j}^{(0)},\delta\};\Lambda_{\mathbf{% \Upsilon}_{j}^{(0)}})={}^{\odot}\mathbf{\Delta}_{j}^{(1)}$ (11)

where $\mathscr{B}(\left\{{\it pop},{\it iter}\right\};{\it ff})$ is the operator that applies a bat algorithm on a population pop for iter iterations, ff is used to indicate the fitness function, and the subscript $j$ indicates the association with the vector $\mathbf{\Upsilon}_{j}^{(0)}$ , which is used to compute the fitness function, denoted here as $\Lambda_{\mathbf{\Upsilon}_{j}^{(0)}}$ , according to Eq. (6). The output is the optimal parameter vector, represented by ${}^{\odot}\mathbf{\Delta}_{j}^{(1)}$ , where the left superscript ${}^{\odot}(.)$ is used onwards to indicate the optimum. Bat algorithm in Eq. (11) can be further decomposed as:

$\begin{array}[]{l}\mathscr{B}_{j}\left(\left\{\bar{\mathbf{\Delta}}_{j}^{(0)},% \delta\right\};\Lambda_{\mathbf{\Upsilon}_{j}^{(0)}}\right)=\\ \hskip 28.452756pt\overset{\delta}{\underset{i=1}{\bigcirc}}\left\{\left[% \Lambda_{\mathbf{\Upsilon}_{j}^{(0)}}\left(\oplus\left({}^{[i]}\mathbf{\Delta}% _{j,l}^{(0)}\right)\right)\right]_{l=1}^{\lambda}\right\}\\ \end{array}$ (12)

where $\overset{\delta}{\underset{i=1}{\bigcirc}}$ is used to represent the composition of iterations $i=1,\dots,\delta$ and the operator $\oplus({\it arg})$ updates the individual arg according to Eqs (7)–(9), i.e., lines 9-10 of Algorithm 1. Equation (12) is useful to better understand how this bat algorithm works. For each iteration $i=1,\dots,\delta$ , the operator $\oplus$ is applied onto ${}^{[i]}\mathbf{\Delta}_{j,l}^{(0)}$ , where the left superscript ${}^{[i]}(.)$ is used to indicate the iterations of the bat algorithm applied to the population $\bar{\mathbf{\Delta}}_{j}^{(0)}$ . The output of this operator, ${}^{[i+1]}\mathbf{\Delta}_{j,l}^{(0)}$ , is used to compute the fitness function through the operator $\Lambda_{\mathbf{\Upsilon}_{j}^{(0)}}$ , which means that Eq. (6) is computed for the pair $({}^{[i+1]}\mathbf{\Delta}_{j,l}^{(0)},\mathbf{\Upsilon}_{j}^{(0)})$ and then probabilistically accepted according to lines 11–19 in Algorithm 1. Then, the population at the last iteration, ${}^{[\delta]}\bar{\mathbf{\Delta}}_{j}^{(0)}$ is evaluated and the best individual, ${}^{\odot}\mathbf{\Delta}_{j}^{(1)}$ , is selected. Finally, take $\mathbf{\Upsilon}_{j}^{(1)}=\mathbf{\Upsilon}_{j}^{(0)}$ .

The output of this first stage is given by the collection of pairs $({}^{\odot}\mathbf{\Delta}_{j}^{(1)},\mathbf{\Upsilon}_{j}^{(1)})$ , for $j=1,\dots,\eta$ .

Stage 2: Main procedure. At this stage, the ME-BAT approach is applied. It consists of applying the bat algorithm at different layers. The upper layer computes the optimal vector of breakpoints through the bat algorithm for a given number of iterations, $\kappa$ . The initial population of breakpoint vectors is given by $\bar{\mathbf{\Upsilon}}^{(1)}=\{\mathbf{\Upsilon}_{j}^{(1)}\}_{j}$ . At each iteration step $g=1,\dots,\kappa$ , each individual $\mathbf{\Upsilon}_{j}^{(1)}$ is modified according to Eqs (7)–(9), then the change is probabilistically accepted according to the improvement of the fitness function and the value of the loudness. Mathematically, this means:

$\mathscr{B}\left(\{\bar{\mathbf{\Upsilon}}^{(1)},\kappa\};(\Lambda_{{}^{\odot}% \mathbf{\Delta}_{j}^{(g+1)}})_{g}\right)={}^{\odot}\mathbf{\Upsilon}$ (13)

where the meaning of all symbols is already known. Further decomposition of this bat algorithm is:

$\begin{array}[]{l}\mathscr{B}\left(\{\bar{\mathbf{\Upsilon}}^{(1)},\kappa\};(% \Lambda_{{}^{\odot}\mathbf{\Delta}_{j}^{(g+1)}})_{g}\right)=\\ \hskip 28.452756pt\overset{\kappa}{\underset{g=1}{\bigcirc}}\left\{[\Lambda_{{% }^{\odot}\mathbf{\Delta}_{j}^{(g+1)}}(\oplus(\mathbf{\Upsilon}_{j}^{(g)}))]_{j% =1}^{\eta}\right\}\\ \end{array}$ (14)

where, once again, all symbols are known. However, the fitness operator $\Lambda_{{}^{\odot}\mathbf{\Delta}_{j}^{(g+1)}}$ requires further explanation. As previously remarked, the fitness function of this reconstruction problem depends not only on the breakpoint vector but also on the parameter vector. Therefore, any update of the breakpoint vector at step $g$ from $\mathbf{\Upsilon}_{j}^{(g)}$ to $\mathbf{\Upsilon}_{j}^{(g+1)}=\oplus(\mathbf{\Upsilon}_{j}^{(g)})$ implies that a new optimal parameter vector has to be computed, for $g=1,\dots,\kappa$ . For instance, think about the parameter vector ${}^{\odot}\mathbf{\Delta}_{j}^{(1)}$ . Although it is optimal for its associated breakpoint vector $\mathbf{\Upsilon}_{j}^{(1)}$ , it is no longer optimal for the updated individual $\mathbf{\Upsilon}_{j}^{(2)}$ . Therefore, a new optimal parameter vector ${}^{\odot}\mathbf{\Delta}_{j}^{(2)}$ is needed. This task is carried out at a lower layer, consisting of another bat algorithm embedded into the main one. Therefore, for each iteration step of the bat algorithm for the breakpoint vector population, a second bat algorithm for a new parameter vector population is executed as:

$\displaystyle{}^{\odot}\mathbf{\Delta}_{j}^{(g+1)}=\mathscr{B}_{j}\left(\{{}^{% [0]}\bar{\mathbf{\Delta}}_{j}^{(g)},\delta\};\Lambda_{{}\mathbf{\Upsilon}_{j}^% {(g+1)}}\right)$ $\displaystyle=\overset{\delta}{\underset{i=1}{\bigcirc}}\left\{[\Lambda_{% \mathbf{\Upsilon}_{j}^{(g+1)}}(\oplus({}^{[i]}\mathbf{\Delta}_{j,l}^{(g)}))]_{% l=1}^{\lambda}\right\}$ (15)

that is, at iteration step $g$ , the the new breakpoint vector $\mathbf{\Upsilon}_{j}^{(g+1)}$ is used to compute the fitness. As result, the best parameter vector ${}^{\odot}\mathbf{\Delta}_{j}^{(g+1)}$ associated with $\mathbf{\Upsilon}_{j}^{(g+1)}$ is selected and the main bat algorithm continues until convergence.

The final output of this process is the best pair of parameter vector and breakpoint vector: $\left({}^{\odot}\mathbf{\Delta},{}^{\odot}\mathbf{\Upsilon}\right)$ .

3.3.3 Pole computation

The output of the previous step is now used to compute the poles of the curve. Using matrix notation and inserting the optimal values for data parameters and breakpoints, Eq. (6) can be written as:

$\underset{\{\mathbf{\Xi}_{i}\}_{i}}{min}\|{\mathbf{\Theta}}-{\bf\Psi}.{\bf\Xi}% \|^{2}$ (16)

where the optimization is now performed only on the poles $\{\mathbf{\Xi}_{i}\}_{i}$ . The necessary condition for ${\mathbf{\Xi}}$ to be the solution of Eq. (16) is that:

${\bf\Psi}^{T}.{\bf\Psi}.{\bf\Xi}={\bf\Psi}^{T}.{\bf\Theta}$ (17)

which leads to the normal equation:

${\bf\Gamma}.{\bf\Xi}={\bf\Pi}$ (18)

where

${\bf\Gamma}={\left[\sum\limits_{k=2}^{\chi-1}\psi_{l}^{\rho}(\xi_{k},\mathbf{% \Upsilon})\psi_{i}^{\rho}(\xi_{k},\mathbf{\Upsilon})\right]}$

and

${\bf\Pi}={\left[\sum\limits_{k=2}^{\chi-1}{\bf\Theta}_{k}\psi_{i}^{\rho}(\xi_{% k},\mathbf{\Upsilon})\right]}$

for $i,l=0,\dots,n$ . Note that ${\bf\Psi}^{T}.{\bf\Psi}$ is a symmetric square matrix and positive semidefinite, so system Eq. (17) always has a solution. It can be solved numerically by Gaussian elimination or by using the singular value decomposition (SVD) (see [28] for details). The algebraic solution of Eq. (17) is given by: ${\bf\Xi}={\bf\Psi}^{+}.{\bf\Theta}$ , where ${\bf\Psi}^{+}$ denotes the Moore-Penrose pseudo-inverse of ${\bf\Psi}$ . If the matrix ${\bf\Psi}$ has rank $\sigma+1$ , ${\bf\Gamma}$ is non-singular (or equivalently, ${\bf\Psi}^{T}.{\bf\Psi}$ is invertible) and Eq. (18) defines ${\bf\Xi}$ uniquely. Also, an explicit formula for ${\bf\Psi}^{+}$ is available as ${\bf\Psi}^{+}=({\bf\Psi}^{T}.{\bf\Psi})^{-1}.{\bf\Psi}^{T}={\bf\Gamma}^{-1}.{% \bf\Psi}^{T}$ . Note that ${\bf\Psi}^{+}$ is then a left inverse of ${\bf\Psi}$ . Otherwise, ${\bf\Psi}^{+}$ is defined by choosing ${\bf\Xi}$ to minimize $\|{\bf\Xi}\|^{2}$ among the solutions of Eq. (17).

3.4 Parameter tuning

Parameter tuning is a critical issue for metaheuristic techniques, as their performance is strongly dependent on the parameter values. Such values are also problem-dependent, so they cannot be determined in advance. The approach in this paper is a mixed one: good values described in the literature are used as an initial seed to guide the process; then, they are refined by numerous computer simulations for different parameter values. As a result, the final choice is necessarily empirical. The different parameters used in this method are arranged in rows in Table 2. For each parameter, the table shows (in columns) its symbol, meaning, range of values, and the parameter value chosen in this paper.

Table 2
Parameters and values used in the ME-BAT algorithm

Symbol	Meaning	Range	Value
$\eta$	Breakpoint population size	20–100	80
$\lambda$	Parameter population size	20–100	50
$\kappa$	# of iterations (main bat)	100–2000	1000
$\delta$	# of iterations (embedded bat)	50–500	100
${\mathcal{A}}^{0}$	Initial loudness	(0, 2)	0.5
${\mathcal{A}}_{\min}$	Minimum loudness	[0, 1]	0
$r^{0}$	Initial pulse rate	[0, 1]	0.5
$f_{\max}$	Maximum frequency	[0, 10]	2
$\alpha$	Multiplicative factor	(0, 1)	0.6
$\gamma$	Exponential factor	[0, 1]	0.4

The most important parameters of the method are:

Population size (for breakpoint vectors, $\eta$ , and parameter vectors, $\lambda$ ): Increasing the population size decreases the number of iterations, but it also increases the number of function evaluations. Therefore, a trade-off between both situations is required. In this work, the populations vary from 20 to 100 bats for the breakpoint and the parameter vectors. We got that a value of 80 for the former and 50 for the later provides the best efficiency for the method. This disparity is because the optimal value of the embedded bat algorithm is less critical, as it changes frequently during execution. So, it is better to take fewer individuals to hasten the execution without loss of quality.

Maximum number of iterations (for breakpoint vectors, $\kappa$ , and parameter vectors, $\delta$ ): Again, two values are considered: one for the main bat algorithm computing the breakpoints and another for the embedded bat algorithm for the data parameters. The first one is more critical, so its range was taken from 100 to 2000 iterations. We found that the method converged in less than 500 iterations in all executions, so this parameter was set to 1000 iterations. For $\delta$ the range analyzed varies from 50 to 500. Since no variations were noticed for values larger than 100, it was finally set to 100.

Initial loudness, ${\mathcal{A}}^{0}$ , minimum loudness, ${\mathcal{A}}_{\min}$ , and multiplicative factor, $\alpha$ : From the simulations we found that results do not change significantly for values of ${\mathcal{A}}^{0}$ in the range (0, 2), and of ${\mathcal{A}}_{\min}$ in the interval [0, 1], meaning that these parameters are very robust against variations on such intervals. So they are set to 0.5 and 0 respectively, as suggested by several authors. We also found that $\alpha=0.6$ is a suitable value for this problem.

Initial pulse rate, $r^{0}$ , and exponential factor, $\gamma$ : of these two parameters, the initial pulse rate is the most relevant. It was set to 0.5, meaning that the selection has an equal probability of change in the long term. On the contrary, parameter $\gamma$ has little relevance in the simulation results, as it only affects the very early iterations. It was set to 0.4 as suggested by previous works.

Figure 1.

Butterfly example: (top) initial noisy data points; (middle) reconstructed curve with the control polygon and the poles; (bottom) best reconstructed B-spline curve.

Figure 2.

Convergence diagram for the butterfly example: RMSE vs. iterations of main bat algorithm, $g$ .

4. Experimental results

The method introduced in previous section has been applied to three real-world examples from different engineering fields: a butterfly model, a hand drill, and a fan blade. In all cases, the input data is a collection of scanned noisy data points. They are fitted through a single parametric B-spline curve. The order of the curve, $\rho$ , and the number of poles, $\sigma$ , are freely chosen by the user. Choice of $\rho$ is important: low-order polynomials give little flexibility in controlling the shape of the curve while high-order polynomials can introduce unwanted wiggles and require more computation; therefore, too low and too high values for $\rho$ must be avoided. A classical choice is to consider fourth-order B-splines. This is also the value taken in this paper. Note, however, that the method does not depend on the value of $\rho$ . The best value of $\sigma$ is determined empirically by taking different values for the computer simulations and selecting that with the best fitting error. The performance of the method is checked by computing the RMSE (root-mean square error), given by:

$\textit{RMSE}={\sqrt{{{\Lambda}\over{\chi}}}}$ (19)

For all examples in this paper, 20 independent executions have been carried out to remove any spurious effect derived from stochasticity.

4.1 Example 1: Butterfly

First example corresponds to the shape of a butterfly model firstly introduced in [40] as an example about real-time NURBS interpolation for high-speed motion control in the context of CNC machine tools for manufacturing. Since then, it has been used as a suitable benchmark instance for various industrial applications [38, 43, 44], including the use of a commercial CNC system to produce a machined contour of this shape [43]. Originally, the curve was described as a cubic NURBS curve with 51 poles and non-unit weights at some corner points to strengthen the pull of the curve towards the control points. A NC code reporting the parameter values of the curve can be found in [40]. This example has been primarily chosen not only because it is a real-world example widely used as a benchmark in previous works but also because it is mathematically described as a NURBS curve. It is therefore a challenging shape for the method introduced in this paper, which does not consider weights.

Figure 3.

Hand drill example: (top) initial noisy data points; (middle) reconstructed curve with the control polygon and the poles; (bottom) best reconstructed B-spline curve.

Figure 4.

Convergence diagram for the hand drill example: RMSE vs. iterations of main bat algorithm, $g$ .

Figure 1(top) shows the collection of 472 scanned data points. Middle and bottom pictures in this figure show the best fitting B-spline curve, obtained for the value $\sigma=118$ . The figure in the middle also shows the poles and the control polygon of the curve, which are removed in the bottom figure for better visualization of the reconstructed shape. As the reader can see, the method is able to recover the underlying shape of the data points. This good visual quality is confirmed by the numerical results, with a RMSE value for the best fitting curve of 4.03041 $\times$ 10 ${}^{-3}$ . Note that, although the proposed method rely exclusively on polynomial functions, the original shape (given by a rational curve) is still reconstructed with high accuracy. Note also that the number of poles is substantially larger than the original curve in [40]. This means that this high accuracy is achieved at the expense of adding more parameters to the model. Figure 2 shows the convergence diagram for this example, depicting the RMSE vs. the number of iterations for the main bat algorithm.

Figure 5.

Whole three-blade fan of Example 3.

4.2 Example 2: Hand drill

Cross-section curves are very important in several applied fields. For example, in the medical area it is possible to reconstruct the outer surface of the organ under study from a set of parallel slices corresponding to different levels. Parametric curves are also used to represent active contours in ultrasound imaging segmentation [23]. Similarly, in CAD/CAM many objects are designed by defining a number of cross-sections; the shape of the surface in-between is computed by some interpolation scheme. The second example corresponds to the profile curve of the central cross-section of a commercial hand drill of a popular manufacturer. Figure 3(top) shows the collection of 466 scanned data points. This example is a good candidate to check the performance of our method against the noise, as data points are affected by measurement noise of varying intensity for different regions. The best fitting curve, obtained for $\sigma=$ 48 is displayed in the middle and bottom parts of this figure, whose meaning is similar to the previous example. Note that the B-spline curve yields smooth shapes and tends to remove the wandering parts of data. This fact is particularly noticeable in the upper left part of the shape in Fig. 3(top). The best RMSE for this example is: 6.17533 $\times$ 10 ${}^{-3}$ . This example shows the robustness of the proposed method against measurement noise. It also shows that the method is able to adapt to complicated shapes with several concavity changes. The convergence of the method is shown in Fig. 4.

Figure 6.

Lower fan blade example: (top) initial noisy data points; (middle) reconstructed curve with the control polygon and the poles; (bottom) best reconstructed B-spline curve.

Table 3

RMSE of the three examples for different noise intensity

Noise intensity rate ${}^{1}$	Butterfly	Hand drill	Fan blade
Initial data	4.03041e-3	6.17533e-3	1.68241e-3
0.5%	4.03283e-3	6.19883e-3	1.68833e-3
1%	4.13238e-3	6.35751e-3	1.70176e-3
2.5%	4.90582e-3	7.31422e-3	1.98423e-3
5%	6.60128e-3	9.79778e-3	2.59529e-3
10%	9.26833e-3	1.33282e-2	3.42791e-3

${}^{1}$ The initial data in first row are already noisy. Values in this column means this extra noise is added to the noisy initial data.

Figure 7.

Convergence diagram for the fan blade example: RMSE vs. iterations of main bat algorithm, $g$ .

4.3 Example 3: Fan blade

Last example corresponds to a commercial three-blade fan shown in Fig. 5. The 271 sampled data points correspond to the lower blade of the fan, displayed in Fig. 6(top). The other two blades of the fan are similar and can readily be obtained by simple rotation, so they are not considered here. The blade is reconstructed with a B-spline curve with 31 poles, displayed in middle and bottom figures, with a very good RMSE of 1.68241 $\times$ 10 ${}^{-3}$ for the best fitting curve. The convergence diagram is shown in Fig. 7.

4.4 Robustness against noise

To check the robustness of our method against noise, the (already noisy) initial data for the three examples have been perturbed by an additive Gaussian white noise of intensity rate $I_{r}$ (in percent) with respect to the initial signal. Then, our method is applied to these new data for five different values of $I_{r}$ , ranging from low intensity ( $I_{r}=$ 0.5%) to high intensity ( $I_{r}=$ 10%). Table 3 shows the RMSE fitting errors for the three examples in this paper. In all cases the RMSE increases by less than one order of magnitude as the noise increases. We conclude that the method is very robust for intensity noise up to 10% of the main signal.

4.5 Implementation issues

All computations in this paper have been performed on a 2.6 GHz. Intel Core i7 processor with 8 GB of RAM. The source code has been implemented by the authors in the popular scientific program Matlab, version 2013b. Equation (17) was solved through specialized Matlab commands for Gaussian elimination with partial pivoting and singular value decomposition (SVD) for squared and non-squared systems, respectively. In this paper, SVD has been primarily used since it provides the best numerical answer when the exact solution is not possible. Matlab also provides routines for ill-conditioned matrices, a situation that can happen in practice, for instance, when one or several singular values in SVD decomposition are null or very near to zero. Advisable answer to this problem is to set reciprocals of such singular values to zero. Matlab command svd handles this problem automatically. Another problem is that, although matrix ${\bf\Gamma}$ is nonnegative, symmetric and banded, its computation explicitly to solve Eq. (18) is not computationally efficient and often a source of numerical rounding errors. A better alternative is to use the QR decomposition of ${\bf\Psi}$ . In case ${\bf\Gamma}$ is well-conditioned and positive definite (that is, it has full rank), Eq. (18) can be solved directly by using the Cholesky decomposition. See [28] for details on these numerical routines and their implementation.

Table 4
RMSE with different methods (best results in bold)

Method for $\{\upsilon_{j}\}_{j}\,,\{\xi_{k}\}_{k}$	Butterfly	Error rate (%) ${}^{1}$	Hand drill	Error rate (%)	Fan blade	Error rate (%)
Uniform/uniform	0.01248773	309.83	0.00927372	150.17	0.00375599	223.25
Arc-length/uniform	0.01154352	286.41	0.00966291	156.47	0.00238136	141.54
Centripetal/uniform	0.01162628	288.46	0.00944385	152.92	0.00292093	173.61
Uniform/averaging	0.01191354	295.59	0.00931234	150.79	0.00384619	128.61
Arc-length/averaging	0.01163537	288.68	0.00906036	146.71	0.00280943	166.98
Centripetal/averaging	0.01157693	287.24	0.00912911	147.83	0.00372863	221.62
Centripetal/ClonalG [30]	0.01178512	292.40	0.01154432	186.94	0.00421173	250.33
PSO/averaging [29]	0.00927325	230.08	0.00875366	141.75	0.00310192	184.37
Arc-length/GMM [46]	0.00877482	217.71	0.00848763	137.44	$\times$	$\times$
Firefly algorithm [8]	0.00961356	238.52	0.00907538	146.96	0.00243651	144.82
PESA [31]	0.01063242	263.80	0.00941121	152.40	0.00385437	229.09
ME-BAT (this method)	0.00403041	–	0.00617533	–	0.00168241	–

${}^{1}$ The error rate (in %) is computed with respect to the RMSE of the best method (last row of the table).

Table 5

CPU time (in minutes) for the three examples in this paper

Method	Butterfly	Hand drill	Fan blade
Centripetal/ClonalG [30]	20–28	130–180	90–110
PSO/averaging [29]	13–16	60–75	50–60
Arc-length/GMM [46]	18–23	90–110	$\times$
Firefly algorithm [8]	14–17	70–90	55–75
PESA [31]	20–30	120–150	75–90
ME-BAT (this method)	22–25	120–130	70–80

4.6 Computational times

It is well-known that metaheuristic techniques are not well suited for real-time applications. Also, it is very difficult to determine the CPU time in advance, as it can vary greatly depending on the complexity of the model, size of input data, population size, number of iterations, hardware and programming language used for execution, and many other factors. Table 5 shows the CPU time (in minutes) for all the metaheuristic methods applied so far to the parametric B-spline curve reconstruction problem, including our method (reported in last row). Although our CPU times might seem very large, they are not surprising at all considering the large number of executions of the bat algorithm required by the multilayer embedded structure of our approach. The comparison with the CPU times of other alternative methods in Table 5 show that they are of similar order as ours, although some are actually better. However, such better CPU times are obtained at the expense of much larger fitting errors (see Section 5 for details). The emphasis of our method is on quality, not on computational speed. Note also that the embedded bat algorithms in Eqs (11) and (12) are executed for parameter vector populations associated with independent breakpoint vectors and similar happens for Eq. (3.3.2). This means that they can be computed simultaneously, so the method is well suited for partial parallelization.

5. Comparison with other methods

As shown in Table 1, most of the previous work about B-spline curve reconstruction reported in the literature is focused on the simplest case: explicit curves. For the parametric case, only partial solutions are provided: in general, either the data parameters or the breakpoints are not optimized but fixed. To the best of authors knowledge, no single metaheuristic technique addressed the general B-spline curve reconstruction problem so far. This fact is not accidental, but a clear indication of the difficulty of this task.

This lack of previous references prevents the comparison with any other general method. It is possible, however, to make a comparison with other methods commonly accepted in the field, even although they provide only partial solutions. This includes the three classical parameterization methods indicated in Section 2.1: uniform, arc-length, and centripetal, along with two methods for breakpoint placement: uniform and averaging. Since solving the general case requires the full computation of data parameters and breakpoints, these partial approaches have to be combined, yielding six feasible couples. They are listed in rows 1–6 in Table 4. Five other methods for parametric B-spline curves (listed in rows 7–11 of Table 4) can be found in the literature. The first one applies an artificial immune system called ClonalG to compute the breakpoints from a centripetal parameterization [30]. The second method computes the breakpoints by the averaging method, while particle swarm optimization (PSO) is applied to the data parameters [29]. The method in [46] relies on the arc-length parameterization and applies an estimation of distribution algorithm (EDA) with Gaussian mixture distributions to compute the breakpoints. This method can only be applied to closed curves. Fourth method applies the firefly algorithm to data parameterization while the breakpoints are obtained by using De Boor’s method [8]. The immune approach in [30] is replaced in [31] by a Pareto Envelope-Based Selection Algorithm (PESA), a type of genetic algorithm for multi-objective optimization. The ME-BAT method is considered in last row.

Table 4 reports the RMSE values for all these methods (in rows) and the three examples (in columns) of this paper. Best results are highlighted in bold for easier identification. We also show (in columns) the error rate (in percent) with respect to our method for better assessment. As the reader can see, the ME-BAT algorithm outperforms all other methods for all instances in the benchmark. Among the classical methods, the uniform parameterization is generally the worst, while the arc-length parameterization is slightly better than the centripetal, but they still perform very similarly. Although based on a powerful metaheuristic technique, the method in [30] does not improve the classical methods. This can be explained by the fact that the method is not continuous, but based on the conversion of the original continuous optimization problem into a combinatorial optimization problem. This conversion process introduces large discretization errors, even for very simple shapes. This discrete approach is also applied in [31], although the fitting errors improve in all cases. The methods in [8, 29, 46] perform pretty well: they improve the results of the classical methods and those in [30, 31] significantly. Among them, the method in [46] is the best for the two first examples. However, it cannot be applied to the third one (an open curve). Still, their results are worse than those of the ME-BAT method. As a conclusion, the method introduced here outperforms all previous methods in the literature on our benchmark.

6. Conclusions and future work

This paper presents a method called multilayer embedded bat algorithm (ME-BAT) to solve the general reconstruction problem with parametric B-spline curves. The method relies on a powerful metaheuristic technique called bat algorithm aimed at solving difficult continuous optimization problems. Opposed to the previous methods in the literature, the proposed method computes the optimal values of all free variables (data parameters, breakpoints, and poles) at full extent. This is very difficult task because these variables are strongly intertwined in a nonlinear and complicated way, so they cannot be obtained independently. The proposed approach is based on the idea of applying the bat algorithm at different layers, with a main bat algorithm at an upper layer to compute the breakpoints and a second bat algorithm at a lower layer to compute the data parameters. This second bat algorithm is embedded into the first one and executed for each individual breakpoint of the population and at each iteration step of the main bat algorithm. Then, the poles are calculated by least-squares minimization through SVD. The method has been applied to three real-world examples from engineering fields. The experimental results show that the method performs very well, being able to recover the underlying shape of data with high accuracy. A comparison with eleven alternative methods (including six classical methods in the field and all the metaheuristic methods applied to this problem so far) shows that the proposed method outperforms previous methods in the literature for the examples in this benchmark.

Main contributions of the paper are:

A new bat algorithm-based scheme (ME-BAT) for continuous optimization where the variables are related to each other in a highly nonlinear and complicated way. This methodology is very general and can readily be adapted for application to any other problem exhibiting similar features.

The application of this scheme to the general curve reconstruction problem with free-form parametric B-splines. Opposed to all previous approaches, the proposed method computes all free variables of the problem at full extent. This is the first reported method able to solve this problem in all its generality.

Main limitations of the method are the computation time and the determination of the optimal number of poles. The first problem can be alleviated by partial parallelization of the embedded bat algorithm, but it is still unclear how much improvement could be obtained. An exciting trend in this regard is the increasing power of GPUs with their highly parallel structure, making them more efficient than CPUs for algorithms where large blocks of code or data can be processed in parallel (see [47, 48] for two recent examples of the use of GPUs to optimize the performance of metaheuristic techniques). The determination of the optimal number of poles might be addressed by adding a careful chosen penalty term to current fitness function. These tasks are part of the future research in the field.

This method could be extended to NURBS at the expense of introducing another set of free variables (the weights) that must be computed as well. The method should be modified to account for these extra variables. Theoretically, the addition of a new layer might probably solve the problem but the computation time could become unacceptable for many practical problems. Clearly, further research is still needed to solve these and other challenging issues related to curve reconstruction. The application of this method to other engineering problems is also part of our future work.

Footnotes

Acknowledgments

This research is supported by the Computer Science National Program of the Spanish Ministry of Economy and Competitiveness, Project Ref. TIN2012-30768 and Toho University (Funabashi, Japan). Special thanks are due to the editors and the four anonymous reviewers for their encouraging and constructive comments and very helpful feedback that allowed us to improve our paper significantly.

References

Boulkaibet

Mthembu

De Lima Neto

Marwala

. Finite element model updating using fish school search and volitive particle swarm optimization. Integrated Computer-Aided Engineering. 2015; 22(4): 361-376.

Burchard

. Splines (with optimal knots) are better. Applicable Analysis. 1974; 3: 309-319.

Cai

. Encryption based partial sharing of CAD models. Integrated Computer-Aided Engineering. 2015; 22(3): 243-260.

Cheng

Zhang

Caraffini

Neri

. Multicriteria adaptive differential evolution for global numerical optimization. Integrated Computer-Aided Engineering. 2015; 22(2): 103-107.

Coletta

LFS

Hruschka

Acharya

Ghosh

. Using metaheuristics to optimize the combination of classifier and cluster ensembles. Integrated Computer-Aided Engineering. 2015; 22(3): 229-242.

Dierckx

. Curve and surface fitting with splines. Oxford University Press, Oxford, 1993.

Gálvez

Iglesias

. Efficient particle swarm optimization approach for data fitting with free knot B-splines. Computer-Aided Design. 2011; 43(12): 1683-1692.

Gálvez

Iglesias

. From nonlinear optimization to convex optimization through firefly algorithm and indirect approach with applications to CAD/CAM. The Scientific World Journal. 2013; 10. Article ID 283919.

Gálvez

Iglesias

. New memetic self-adaptive firefly algorithm for continuous optimization. International Journal of Bio-Inspired Computation. 2016; 8(5): 300-317.

10.

Gálvez

Iglesias

Avila

Otero

Arias

Manchado

. Elitist clonal selection algorithm for optimal choice of free knots in B-spline data fitting. Applied Soft Computing. 2015; 26: 90-106.

11.

Goldenthal

Bercovier

. Spline curve approximation and design by optimal control over the knots. Computing. 2004; 72: 53-64.

12.

Yan

. Neural network approach to the reconstruction of free-form surfaces for reverse engineering. Computer-Aided Design. 1995; 27(1): 59-64.

13.

Hoffmann

. Numerical control of Kohonen neural network for scattered data approximation. Numerical Algorithms. 2005; 39: 175-186.

14.

Hölzle

. Knot adjustment for piecewise polynomial approximation of curves. Computer Aided Design. 2003; 15(5): 295-296.

15.

Iglesias

Echevarría

Gálvez

. Functional networks for B-spline surface reconstruction. Future Generation Computer Systems. 2004; 20(8): 1337-1353.

16.

Jing

Sun

. Fitting B-spline curves by least squares support vector machines. Proc of the 2nd, Int Conf on Neural Networks & Brain. Beijing (China). IEEE Press, 2005; 905-909.

17.

Jupp

DLB

. Approximation to data by splines with free knots. SIAM Journal of Numerical Analysis. 1978; 15: 328-343.

18.

Zhao

Goh

. Adaptive knot placement in B-spline curve approximation. Computer-Aided Design. 2005; 37: 791-797.

19.

Liu

Gao

Shen

. An integrated feature-based dynamic control system for on-line machining, inspection and monitoring. Integrated Computer-Aided Engineering. 2015; 22(2): 187-200.

20.

Lostado

Martinez

Donald

BJM

Villanueva

. Combining soft computing techniques and the finite element method to design and optimize complex welded products. Integrated Computer-Aided Engineering. 2015; 22(4): 153-170.

21.

Kruth

. Parameterization of randomly measured points for least squares fitting of B-spline curves and surfaces. Computer Aided Design. 1995; 27(9): 663-675.

22.

Mencía

Sierra

Mencía

Varela

. Genetic algorithms for the scheduling problem with arbitrary precedence relations and skilled operators. Integrated Computer-Aided Engineering. 2016; 23(3): 269-285.

23.

Pan

Yuan

. Using shapes correlation for active contour segmentation of uterine fibroid ultrasound images in computer-aided therapy. Applied Mathematics – A Journal of Chinese Universities. 2016; 31(1): 37-52.

24.

Park

. An error-bounded approximate method for representing planar curves in B-splines. Computer Aided Geometric Design. 2004; 21: 479-497.

25.

Park

Lee

. B-spline curve fitting based on adaptive curve refinement using dominant points. Computer-Aided Design. 2007; 39: 439-451.

26.

Patrikalakis

Maekawa

. Shape interrogation for computer aided design and manufacturing. Springer Verlag, Heidelberg: 2002.

27.

Piegl

Tiller

. The NURBS Book. Springer Verlag, Berlin Heidelberg: 1997.

28.

Press

Teukolsky

Vetterling

Flannery

. Numerical Recipes (2nd edition). Cambridge University Press, Cambridge, 1992.

29.

Sun

Tao

Wei

Xia

. B-spline curve fitting based on adaptive particle swarm optimization algorithm. Applied Mechanics and Materials. 2010; 20-23: 1299-1304.

30.

Ulker

Arslan

. Automatic knot adjustment using an artificial immune system for B-spline curve approximation. Information Sciences. 2009; 179: 1483-1494.

31.

Ulker

. B-Spline curve approximation using Pareto envelope-based selection algorithm-PESA. International Journal of Computer and Communication Engineering. 2013; 2(1): 60-63.

32.

Wang

Pottmann

Liu

. Fitting B-spline curves to point clouds by curvature-based squared distance minimization. ACM Transactions on Graphics. 2006; 25(2): 214-238.

33.

Tseng

JCR

Tsai

. A hybrid linear text segmentation algorithm using hierarchical agglomerative clustering and discrete particle swarm optimization. Integrated Computer-Aided Engineering. 2014; 21(1): 35-46.

34.

Zhang

. Service-oriented feature-based data exchange for cloud-based design and manufacturing. IEEE Transactions on Services Computing. (in press) DOI 10.1109/TSC.2015.2501981.

35.

Yang

. A new metaheuristic bat-inspired algorithm. In: Gonzalez

, et al., Nature Inspired Cooperative Strategies for Optimization (NISCO 2010). Studies in Computational Intelligence, Springer Berlin, 284, Springer, 2010; 65-74.

36.

Yang

Gandomi

. Bat algorithm: a novel approach for global engineering optimization. Engineering Computations. 2012; 29(5): 464-483.

37.

Yang

. Bat algorithm: Literature review and applications. Int J Bio-Inspired Computation. 2013; 5(3): 141-149.

38.

Yang

Shen

Yuan

Gao

. Curve fitting and optimal interpolation for CNC machining under confined error using quadratic B-splines. Computer-Aided Design. 2015; 66: 62-72.

39.

Yang

Emmerich

Bäck

Kok

. Multi-objective inventory routing with uncertain demand using population-based metaheuristics. Integrated Computer-Aided Engineering. 2016; 23(3): 205-220.

40.

Yau

Lin

Tsai

. Real-time NURBS interpolation using FPGA for high speed motion control. Computer-Aided Design. 2006; 38: 1123-1133.

41.

Yoshimoto

Moriyama

Harada

. Automatic knot adjustment by a genetic algorithm for data fitting with a spline. Proc of Shape Modeling International ’99, IEEE Computer Society Press, 1999; 162-169.

42.

Yoshimoto

Harada

Yoshimoto

. Data fitting with a spline using a real-coded algorithm. Computer-Aided Design. 2003; 35: 751-760.

43.

Zhang

Yuan

Gao

. Efficient algorithm for time-optimal feedrate planning and smoothing with confined chord error and acceleration. Int. Journal of Advanced Manufacturing Technology. 2013; 66(9-12): 1685-1697.

44.

Zhang

Guo

Gao

. Cubic spline trajectory generation with axis jerk and tracking error constraints. International Journal of Precision Engineering and Manufacturing. 2013; 14(7): 1141-1146.

45.

Zhang

Han

. Quantitative optimization of interoperability during feature-based data exchange. Integrated Computer-Aided Engineering. 2015; 23(1): 31-50.

46.

Zhao

Zhang

Yang

. Adaptive knot adjustment using a GMM-based continuous optimization algorithm in B-spline curve approximation. Computer-Aided Design. 2011; 43: 598-604.

47.

Zhou

Qiu

. Optimization of parallel iterated local search algorithms on graphics processing unit. Journal of Supercomputing. 2016; 72(6): 2394-2416.

48.

Zhou

Qiu

. Dynamic strategy based parallel ant colony optimization on GPUs for TSPs. Science China Information Sciences. 2017; 60: 068102:1-068102:3.

49.

Zhu

Xiao

Wang

Zhang

. Global path planning of wheeled robots using multi-objective memetic algorithms. Integrated Computer-Aided Engineering. 2015; 22(4): 387-404.

Multilayer embedded bat algorithm for B-spline curve reconstruction

Abstract

Keywords

1. Introduction

2. Related work

Table 1 Metaheuristic methods applied to B-spline curve reconstruction.

3.1 The bat algorithm

3.3 Description of the method

3.3.1 Input data

3.3.2 Breakpoint placement – data parameterization

Table 2 Parameters and values used in the ME-BAT algorithm

4.4 Robustness against noise

4.5 Implementation issues

Table 4 RMSE with different methods (best results in bold)

5. Comparison with other methods

6. Conclusions and future work

Footnotes

Acknowledgments

References

Table 1
Metaheuristic methods applied to B-spline curve reconstruction.

Table 2
Parameters and values used in the ME-BAT algorithm

Table 4
RMSE with different methods (best results in bold)