#### Abstract

This paper introduces a new method to compute the approximating explicit B-spline curve to a given set of noisy data points. The proposed method computes all parameters of the B-spline fitting curve of a given order. This requires to solve a difficult continuous, multimodal, and multivariate nonlinear least-squares optimization problem. In our approach, this optimization problem is solved by applying the firefly algorithm, a powerful metaheuristic nature-inspired algorithm well suited for optimization. The method has been applied to three illustrative real-world engineering examples from different fields. Our experimental results show that the presented method performs very well, being able to fit the data points with a high degree of accuracy. Furthermore, our scheme outperforms some popular previous approaches in terms of different fitting error criteria.

#### 1. Introduction

The problem of recovering the shape of a curve/surface, also known as *curve/surface reconstruction*, has received much attention in the last few years [1–14]. A classical approach in this field is to construct the curve as a cross-section of the surface of an object. This is a typical problem in many research and application areas such as medical science and biomedical engineering, in which a dense cloud of data points of the surface of a volumetric object (an internal organ, for instance) is acquired by means of noninvasive techniques such as computer tomography, magnetic resonance imaging, and ultrasound imaging. The primary goal in these cases is to obtain a sequence of *cross-sections* of the object in order to construct the surface passing through them, a process called *surface skinning*.

Another different approach consists of reconstructing the curve directly from a given set of data points, as it is typically done in *reverse engineering* for computer-aided design and manufacturing (CAD/CAM), by using 3D laser scanning, tactile scanning, or other digitizing devices [15, 16]. Depending on the nature of these data points, two different approaches can be employed: *interpolation* and *approximation*. In the former, a parametric curve is constrained to pass through all input data points. This approach is typically employed for sets of data points that are sufficiently accurate and smooth. On the contrary, approximation does not require the fitting curve to pass through the data points, but just close to them, according to prescribed distance criteria. Such a distance is usually measure along the normal vector to the curve at that point. The approximation approach is particularly well suited when data are not exact but subjected to measurement errors. Because this is the typical case in many real-world industrial problems, in this paper we focus on the approximation scheme to a given set of noisy data points.

There two key components for a good approximation: a proper choice of the approximating function and a suitable parameter tuning. The usual models for curve approximation are the free-form curves, such as Bézier, B-spline, and NURBS [17–26]. In particular, B-splines are the preferred approximating functions because of their flexibility, their nice mathematical properties, and its ability to represent accurately a large variety of shapes [27–29]. The determination of the best parameters of such a B-spline fitting curve is troublesome, because the choice of knot vectors has proved to be a strongly nonlinear continuous optimization problem. It is also multivariate, as it typically involves a large number of unknown variables for a large number of data points, a case that happens very often in real-world examples. It is also overdetermined, because we expect to obtain the approximating curve with many fewer parameters that the number of data points. Finally, the problem is known to be multimodal; that is, the least-squares objective function can exhibit many local optima [30, 31], meaning that the problem might have several (global and/or local) good solutions.

##### 1.1. Aims and Structure of This Paper

To summarize the previous paragraphs, our primary goal is to compute accurately all parameters of the approximating B-spline curve to a set of input data points. Unfortunately, we are confronted with a very difficult overdetermined, continuous, multimodal, and multivariate nonlinear optimization problem. In this paper, we introduce a new method to solve this fitting problem by using explicit free-form polynomial B-spline curves. Our method applies a powerful metaheuristic nature-inspired algorithm, called firefly algorithm, introduced recently by professor Xin-She Yang (Cambridge University) to solve optimization problems. The trademark of the firefly algorithm is its search mechanism, inspired by the social behavior of the swarms of fireflies and the phenomenon of bioluminescent communication. The paper shows that this approach can be effectively applied to obtain a very accurate explicit B-spline fitting curve to a given set of noisy data points, provided that an adequate representation of the problem and a proper selection of the control parameters are carried out. To check the performance of our approach, it has been applied to some illustrative examples of real-world engineering problems from different fields. Our results show that the method performs very well, being able to yield a very good fitting curve with a high degree of accuracy.

The structure of this paper is as follows: in Section 2 previous work for computing the explicit B-spline fitting curve to a set of data points is briefly reported. Then, some basic concepts about B-spline curves and the optimization problem to be solved are given in Section 3. Section 4 describes the fundamentals of the firefly algorithm, the metaheuristic used in this paper. The proposed firefly method for data fitting with explicit B-spline curves is described in detail in Section 5. Then, some illustrative examples of its application to real-world engineering problems from different fields along with the analysis of the effect of the variation of firefly parameter on the method performance and some implementation details are reported in Section 6. A comparison of our approach with other alternative methods is analyzed in detail in Section 7. The paper closes with the main conclusions of this contribution and our plans for future work in the field.

#### 2. Previous Work

Many previous approaches have addressed the problem of obtaining the fitting curve to a given set of data points with B-spline curves (the reader is kindly referred to [32–34] for a comprehensive review of such techniques). Unfortunately, many works focus on the problem of data parameterization and they usually skip the problem of computing the knots. But since in this work we consider explicit B-spline curves as approximating functions, our problem is just the opposite; the focus is not on data parameterization but on the knot vector. Two main approaches have been described to tackle this issue: fixed knots and free knots. In general, fixed-knot approaches proceed by setting the number of knots *a priori* and then computing their location according to some prescribed formula [6, 9, 18, 19, 22]. The simplest way to do it is to consider equally spaced knots, the so-called *uniform knot vector*. Very often, this approach is not appropriate as it may lead to singular systems of equations and does not reflect the real distribution of data. A more refined procedure consists of the *de Boor method* and its variations [32], which obtain knot vectors so that every knot span contains at least one parametric value.

However, it has been shown since the 70s that the approximation of functions by splines improves dramatically if the knots are treated as free variables of the problem [30, 31, 35–38]. This is traditionally carried out by changing the number of knots through either knot insertion or knot removal. Usually, these methods require terms or parameters (such as tolerance errors, smoothing factors, cost functions, or initial guess of knot locations) whose determination is often based on subjective factors [31, 39–43]. Therefore, they fail to *automatically* generate a good knot vector. Some methods yield unnecessary redundant knots [44]. Finally, other approaches are generally restricted to smooth data points [20, 45–47] or closed curves [11].

In the last few years, the application of Artificial Intelligence techniques has allowed the researchers and practitioners to achieve remarkable results regarding this data fitting problem. Most of these methods rely on some kind of neural networks [1, 8] and its generalization, the functional networks [3–5, 7, 21, 23, 48]. Other approaches are based on the application of metaheuristic techniques, which have been intensively applied to solve difficult optimization problems that cannot be tackled through traditional optimization algorithms. Recent schemes in this area involve particle swarm optimization [10, 13, 24], genetic algorithms [12, 49–51], artificial immune systems [25, 52], estimation of distribution algorithms [11], and hybrid techniques [26]. As it will be discussed later on the approach introduced in this paper; also it belongs to this class of methods (see Sections 4 and 5 for more details).

#### 3. Basic Concepts and Definitions

##### 3.1. B-Spline Curves

Let be a nondecreasing sequence of real numbers on the interval called *knots*. is called the *knot vector*. Without loss of generality, we can assume that . The th* B-spline basis function ** of order * (or equivalently, degree ) can be defined by the Cox-de Boor recursive relations [38] as
for with . Note that th B-spline basis function of order 1 in (1) is a piecewise constant function with value 1 on the interval called the *support* of and zero elsewhere. This support can be either an interval or reduced to a point, as knots and must not necessarily be different. The number of times a knot appears in the knot vector is called the *multiplicity* of the knot and has an important effect on the shape and properties of the associated basis functions (see, for instance, [29, 32] for details). If necessary, the convention in (2) is applied.

A* polynomial B-spline curve of order * (or degree ) is a piecewise polynomial curve given by
where are coefficients called *control points* as they roughly determine the shape of the curve and are the basis functions defined above. For a proper definition of a B-spline curve in (3), the following relationship must hold: (see [32] for further details).

In many practical applications, the first and last knots of are repeated as many times as the order as follows: , (such a knot vector is usually called a *nonperiodic knot vector*). In general, a B-spline curve does not interpolate any of its control points; interpolation only occurs for nonperiodic knot vectors. In such a case, the B-spline curve does interpolate the first and last control points [32]. Since they are the most common in computer design and manufacturing and many other fields, in this work we will restrict ourselves to nonperiodic knot vectors. Note however that our method does not preclude any other kind of knot vectors to be used instead.

##### 3.2. The Optimization Problem

Similar to [49], let us assume that the data to be fitted are given as where is the underlying (unknown) function of the data and is the measurement error and we assume that . In this regard, a convenient model function, , for (4) is given by (3), where both the vector of coefficients and the vector of knots are considered as tuning parameters and denotes the transpose of a vector or matrix. In the least-squares minimization of problem (4), we determine the control points of the B-spline curve by minimizing the least-squares error, , defined as the sum of squares of the residuals as follows: Note that the objective function to be optimized, , depends on the B-spline coefficients and the interior knots. This problem is well known to be a continuous multimodal and multivariate nonlinear minimization problem [31]. In order to address this problem, in this paper we apply a firefly algorithm and obtain an accurate choice of both and .

#### 4. Firefly Algorithm

The firefly algorithm is a metaheuristic nature-inspired algorithm introduced in 2008 by Yang to solve optimization problems [53, 54]. The algorithm is based on the social flashing behavior of fireflies in nature. The key ingredients of the method are the variation of light intensity and formulation of attractiveness. In general, the attractiveness of an individual is assumed to be proportional to their brightness, which in turn is associated with the encoded objective function.

In the firefly algorithm, there are three particular idealized rules, which are based on some of the major flashing characteristics of real fireflies [53]. They are as follows.(1) All fireflies are unisex, so that one firefly will be attracted to other fireflies regardless of their sex.(2) The degree of attractiveness of a firefly is proportional to its brightness, which decreases as the distance from the other firefly increases due to the fact that the air absorbs light. For any two flashing fireflies, the less brighter one will move towards the brighter one. If there is not a brighter or more attractive firefly than a particular one, it will then move randomly.(3) The brightness or light intensity of a firefly is determined by the value of the objective function of a given problem. For instance, for maximization problems, the light intensity can simply be proportional to the value of the objective function.

The distance between any two fireflies and , at positions and , respectively, can be defined as a Cartesian or Euclidean distance as follows: where is the th component of the spatial coordinate of the th firefly and is the number of dimensions.

In the firefly algorithm, as attractiveness function of a firefly one should select any monotonically decreasing function of the distance to the chosen firefly, for example, the exponential function as follows: where is the distance defined in (6), is the initial attractiveness at , and is an absorption coefficient at the source which controls the decrease of the light intensity.

The movement of a firefly which is attracted by a more attractive (i.e., brighter) firefly is governed by the following evolution equation: where the first term on the right-hand side is the current position of the firefly, the second term is used for considering the attractiveness of the firefly to light intensity seen by adjacent fireflies, and the third term is used for the random movement of a firefly in case there are not any brighter ones. The coefficient is a randomization parameter determined by the problem of interest, while is a random number generator uniformly distributed in the space .

#### 5. The Proposed Method

In this section, the firefly algorithm described in previous paragraphs is applied to obtain the approximating explicit B-spline curve of a given set of noisy data points expressed according to (3)-(4). To this aim, in this work we consider four different fitness functions. The first objective function corresponds to the evaluation of the least-squares function, given by (5). Since this error function does not consider the number of data points, we also compute the RMSE (root-mean squared error) given by

Note, however, that the best error might be obtained in either (5) or (9) at the expense of an unnecessarily large number of variables. To overcome this limitation, we also compute two other fitness functions: *Akaike Information Criterion* (AIC) [55, 56] and *Bayesian Information Criterion* (BIC) [57]. Both are penalized information-theoretical criteria aimed at finding the best approximating model to the true data through an adequate trade-off between fidelity and simplicity. They are comprised of two terms: the first one accounts for the fidelity of the model function, while the second one is a penalty term to minimize the number of free parameters of the model, . They are given, respectively, by

A simple comparison between (10) shows that both criteria are quite similar, but BIC applies a larger penalty than AIC. Thus, as other factors are equal, it tends to select simpler models than AIC (i.e., models with fewer parameters). Both methods provide very similar results for smooth underlying functions. Differences might arise, however, for functions exhibiting sharp features, discontinuities, or cusps. In particular, AIC tends to yield unnecessary redundant knots and, therefore, BIC becomes more adequate for such models. A great advantage of AIC and BIC criteria is that they avoid the use of subjective parameters such as error bounds or smoothing factors. Furthermore, they provide a simple and straightforward procedure to determine the best model: the smaller their value, the better fitness. Because of that, we use them to select the best model for our examples and to compare our results to those of other alternative approaches.

Before searching a solution to the problem, two sets of control parameters must be set up. The first set corresponds to the curve parameters: the number of control points and the order of the B-spline curve. In general, they are freely chosen by the user according to the complexity of the shape of the underlying function of data and the requirements of the particular problem under study, with the only constraint that . The second set of parameters relates to the firefly algorithm itself. As usual when working with nature-inspired metaheuristic techniques, the choice of control parameters is very important, as it strongly influences the performance of the method. Although some guidelines are given in the literature to tackle this issue, such a selection is problem-dependent and, therefore, it remains empirical to a large extent. In this paper, our parameter tuning is based on a combination of both theoretical results found in the literature of the field and a large collection of empirical results. From the former, it is well known that the value of absorption coefficient, , is important in determining the speed of convergence and the behavior of the method. In theory, , but in practice values in the extremes should be avoided. When the attractiveness becomes constant, which is equivalent to say that the light intensity does not decrease with distance, so the fireflies are visible anywhere, while that corresponds to the case where the attractiveness function in (7) becomes a Dirac -function. This case is equivalent to the situation where the fireflies move in a very foggy environment, where they cannot see the other fireflies and hence they move randomly. In other words, this corresponds to the case of a completely random search method. Based on this analysis, we carried out thousands of simulations to determine a suitable value for and found that the value provides a quick convergence of the algorithm to the optimal solution. Regarding the initial attractiveness, , some theoretical results suggest that is a good choice for many optimization problems. We also take this value in this paper, with very good results, as it will be discussed in next section. For the potential coefficient, , any positive value can be used. But it is noticed that the light intensity varies according to the inverse square law, so we set up accordingly. Finally, a stochastic component is necessary in the firefly method in order to allow new solutions appear and avoid to getting stuck in a local minimum. However, larger values introduces large perturbations on the evolution of the firefly and, therefore, delay convergence to the global optima. In this work, we carried out several simulations to determine the optimal value for this parameter and finally we set up its value of for this randomization parameter. The reader is referred to Section 6.4 for a more detailed discussion about how the variation of this parameter affects the behavior of our method.

Once the parameters have been set up, the firefly algorithm is executed for a given number of fireflies, representing potential solutions of the optimization problem. Each firefly in our method is a real-valued vector of length , initialized with uniformly distributed random numbers on the parametric domain and then sorted in increasing order. In this work, an initial population, , of fireflies is considered. Its value is set up to in all examples of this paper. We also tried larger populations of fireflies (up to individuals) but found that our results do not change significantly. Since larger populations mean larger computation times with no remarkable improvement at all, we found this value to be appropriate in our simulations.

This initial population is then subjected to the firefly mechanism and least-squares minimization, yielding a collection of new solutions to the optimization problem (5). In our executions we compute independently all four fitness functions indicated above, meaning that our simulations are to be repeated for each fitness function. On the other hand, each simulation is repeated several times according to a parameter in order to remove the spurious effects of randomization typically arising in any stochastic process. In this paper, we consider a value of , meaning that each simulation is independently carried out 80 times to obtain the final results. Each execution is performed for a given number of iterations, . In general, the firefly algorithm does not need a large number of iterations to reach the global optima. This also happens in this case. In all our computer simulations, we found empirically that is a suitable value, and larger values for this parameter do not improve our simulation results. Once the algorithm is executed for the given number of iterations, the firefly with the best (i.e., minimum) fitness value is selected as the best solution to the optimization problem.

#### 6. Experimental Results

The method described in the previous section has been applied to several examples from different fields. To keep the paper in manageable size, in this section we describe only three of them: the response of an RC circuit, a volumetric moisture content example, and the side profile section of a car body. These examples have been carefully chosen for two reasons: firstly, to show the diversity of situations to which our algorithm can be applied, and secondly, because they correspond to interesting real-world engineering problems rather than to academic examples. As we will show below, our method provides a very accurate solution to these problems in terms of explicit B-spline curves. We think we have provided enough examples to convince our readers about the wide range of applicability of our approach.

##### 6.1. Response of an RC Circuit

The resistor-capacitor (RC) circuit is an electronic circuit comprised of resistors and capacitors driven by a voltage or current source. RC circuits are very popular, the most common being the high-pass and low-pass filters. The simplest version is the first-order RC circuit, where just one capacitor and one resistor are considered. In this example, we consider a signal of data points of a series first-order RC circuit. They are depicted as a collection of red cross symbols in Figure 1. We apply our method to this sequence of data points by using a third-order B-spline curve with . The best fitting curve we obtain is displayed as a blue solid line in Figure 1. This example has been primarily chosen because the underlying curve of data shows a qualitatively different behavior at both ends, with a highly oscillating region on the left and a steady-state behavior on the right. All these features make this function a good candidate to check the performance of our approach. Even in this case, our method yields a fitting of the curve to the data points. The corresponding fitting errors for the best execution and the average of 20 executions are reported in the last column of Table 3. They are in good agreement with the visual results in the figure about the accurate fitting of data points.

##### 6.2. Volumetric Moisture Content

This example uses a set of measurements of a volumetric moisture content, as described in [41]. This example has been traditionally used as a benchmark to test different curve approximation methods (see [41] for details). Figure 2 shows the collection of data points along with their best fourth-order B-spline fitting curve for obtained with our firefly-based method for this data set. The figure clearly shows the good approximation of data points. Note that the fitting curve does not generally pass through the data points (i.e., it is truly an approximating curve, not an interpolating one). Note also that the fitting curve captures the tendency of points very well even with only a few parameters, thus, encoding economically the primary information of data. The corresponding fitting errors for the best execution and the average of 20 executions are again reported in the last column of Table 3.

##### 6.3. Side Profile Section of a Car Body

This example consists of the data fitting of a set of noisy data points from the *In* ( axis) side profile section of a clay model of a notchback three-box sedan car body. The data points were obtained by a layout machine by the Spanish car maker provider Candemat years ago. Data points of the side profile section are sorted in standard order, that is, from forward ( axis) on the left to backward ( axis) on the right [58]. This example has been chosen because it includes areas of varying slope, ranging from the strong slope at both ends (upwards on the left, downwards on the right) to the soft slope in middle part, with a horizontal area between pillars and and in the cargo box area, and a soft upward area in the car hood. Consequently, it is also a very good candidate to check the performance of our approach.

We have applied our firefly-based fitting method by using a fourth-order B-spline curve with . The best fitting curve along with the data points is displayed in Figure 3. Once again, note the excellent fitting of the curve to the data points. The corresponding fitting errors for the best execution and the average of 20 executions are again reported in the last column of Table 3.

##### 6.4. Analysis of Variation of Parameter

In this section we present a detailed discussion about the effects of the variation of the parameter on the performance of our method. To this purpose, we carried out several simulation trials and computed the fitting error for the eight error criteria used in this work. In this section, we show the results obtained for the first and third examples of this paper.Results for the second example are not reported here because the low number of input data points makes it particularly resilient to variations of parameter . As a result, the fitting errors for different values of are very similar to each other and to those reported in the second row of Table 2.

Table 1 shows the variation of the eight fitting error criteria for the RC circuit example and parameter varying from to with step-size . Once again, we carried out 20 independent simulations for each value of and obtained the best and the average values of the 20 executions. Best results for each error criterion are highlighted in bold for the sake of clarity. We can see that the values and lead to the best results for the best and average fitting error values, respectively. Figure 4 shows graphically, from left to right and top to bottom, the best (in blue) and average (in red) values of the , RMSE, AIC, and BIC fitting errors, respectively. From the figure, we can conclude that both and are good choices, as they lead to the best results for this example. Note that values lower and larger than this optimal values yield worse fitting error results. This effect is particularly noticeable for values of approaching to .

**(a)**

**(b)**

**(c)**

**(d)**

Table 2 and Figure 5 show the corresponding results for the side profile section of a car body example. As the reader can see, the value is again a very good choice, although and also yields good results. Once again, increasing the value of towards leads to a strong increase of the fitting errors.

**(a)**

**(b)**

**(c)**

**(d)**

It is worthwhile to recall that the choice of parameters is problem-dependent, meaning that our choice of parameters might not be adequate for other examples. Therefore, it is always advisable to carry out a preliminary analysis and computer simulations to determine a proper parameter tuning for better performance, as it has been done in this section.

##### 6.5. Computational Issues

All computations in this paper have been performed on a 2.9 GHz. Intel Core i7 processor with 8 GB of RAM. The source code has been implemented by the authors in the native programming language of the popular scientific program *Matlab*, version 2010b. Regarding the computation times, all examples we tested can be obtained in less than a second (excluding rendering). Obviously, the CPU time depends on several factors, such as the number of data points and the values of the different parameters of the method, making it hard to determine how long does it take for a given example to be reconstructed.

#### 7. Comparison with Other Approaches

In this section a careful comparison of our method with other alternative methods is presented. In particular, we compare the results of our method with the two most popular and widely used methods: the uniform and the de Boor methods [28, 29, 32]. These methods differ in the way components of the knot vector are chosen. The uniform method yields equally spaced internal knots as This knot vector is simple but does not reflect faithfully the distribution of points. To overcome this problem, a very popular choice is given by the de Boor method as where and . Finally, in our method the knots are determined by the best solution of the optimization problem described in Section 3.2, solved with the firefly-based approach described in Section 5. In this section we have applied these three methods to the three examples of this paper.

Table 3 summarizes our main results. The three examples analyzed are arranged in rows with the results for each example in different columns for the compared methods: uniform method, de Boor method, and our method, respectively. For each combination of example and method, eight different fitting errors are reported, corresponding to the best and average value of the , RMSE, AIC, and BIC error criteria. As mentioned above, the average value is computed from a total of 20 executions for each simulation. The best result for each error criterion is highlighted in bold for the sake of clarity.

It is important to remark here that the error criteria used in this paper for comparison must be carefully analyzed. While the general rule for them is “the lower and the better," this is exclusively true to compare results obtained for experiments carried out *under the same conditions*, meaning that we cannot compare the results for the three examples globally, as they are not done under the same conditions (for instance, they have different number of data points, control points, and even the order of the curve). Therefore, we cannot conclude that the fitting curve for the third example (with the lower global AIC/BIC) fits better than the other ones. The only feasible comparison is among the different methods for the same curve under the same simulation conditions. Having said that, some important conclusions can be derived from the numerical results reported in Table 3.(i) First and foremost, they confirm the good results of our method for the examples discussed in this paper. The fitting errors are very low in all cases, meaning that the fitting curve approximates the given data very accurately. (ii) On the other hand, the numerical results clearly show that our method outperforms the most popular alternative methods in all cases. The proposed approach provides the best results for all error criteria and for all examples, both for the best and for the average values. (iii) Note also that the uniform approach is the worst in all examples, as expected because of its inability to account for the distribution of data points of the given problem. This problem is partially alleviated by the de Boor approach but the errors can still be further improved. This is what our method actually does.

To summarize, the numerical results from our experiments confirms the very accurate fitting to data points obtained with the proposed method as well as its superiority over the most popular alternatives in the field. This superior performance can be explained from a theoretical point of view. Firstly, we notice that the uniform method does not reflect the real distribution of points; in fact, (11) shows that the knots depend exclusively on and , and hence it is independent of the input data. Consequently, it lacks the ability to adapt to the shape of the underlying function of data. Regarding the de Boor method, the choice of knots is based on ensuring that the system of equations is well conditioned and has solution on optimizing the process. Accordingly, it focuses on the numerical properties of the matrices involved. From (4) and (5), we can obtain the following vector equation: where and is the matrix of basis functions sampled at data point parameters. Premultiplying at both sides by , we obtain where is a real symmetric matrix and, therefore, Hermitian, so all its eigenvalues are real. In fact, up to a choice of an orthonormal basis, it is a diagonal matrix, whose diagonal entries are the eigenvalues. It is also a positive semidefinite matrix, so the system has always solution. It can be proved that the knots in the de Boor method are chosen so that it is guaranteed that every knot span contains at least one parametric value [59]. This condition implies that the is positive definite and well conditioned. As a consequence, the system (14) is nonsingular and can be solved by computing the inverse of matrix . As the reader can see, the choice of knots in this method is motivated by the numerical properties of matrix , and never attempted to formulate the fitting process as an optimization problem. In clear contrast, in our approach the problem is formulated as an optimization problem, given by In other words, in our method the emphasis is on optimizing the choice of knots rather than on the numerical properties of the system equation. Note that our method allows consecutive nodes to be as close as necessary; they can even be equal. A simple visual inspection of (11) and (12) shows that this choice is not allowed in these alternative methods. However, near knots are sometimes required for very steep upwards/downwards areas or areas of sharp changes of concavity, such as those in the figures of this paper. These facts justify why our method outperforms both the uniform and the de Boor methods for the computation of the knot vector.

#### 8. Conclusions and Future Work

This paper addresses the problem of curve fitting to noisy data points by using explicit B-spline curves. Given a set of noisy data points, the goal is to compute all parameters of the approximating explicit polynomial B-spline curve that best fits the set of data points in the least-squares sense. This is a very difficult overdetermined, continuous, multimodal, and multivariate nonlinear optimization problem. Our proposed method solves it by applying the firefly algorithm, a powerful metaheuristic nature-inspired algorithm well suited for optimization. The method has been applied to three illustrative real-world engineering examples from different fields. Our experimental results show that the presented method performs very well, being able to fit the data points with a high degree of accuracy. A comparison with the most popular previous approaches to this problem is also carried out. It shows that our method outperforms previous approaches for the examples discussed in this paper for eight different fitting error criteria.

Future work includes the extension of this method to other families of curves, such as the parametric B-spline curves and NURBS. The extension of these results to the case of explicit surfaces is also part of our future work. Finally, we are also interested to apply our methodology to other real-world problems of interest in industrial environments.

#### Acknowledgments

This research has been kindly supported by the Computer Science National Program of the Spanish Ministry of Economy and Competitiveness, Project Ref. no. TIN2012-30768, Toho University (Funabashi, Japan), and the University of Cantabria (Santander, Spain). The authors are particularly grateful to the Department of Information Science of Toho University for all the facilities given to carry out this research work. Special thanks are owed to the Editor and the anonymous reviewers for their useful comments and suggestions that allowed us to improve the final version of this paper.