Abstract
A classical issue in many applied fields is to obtain an approximating surface to a given set of data points. This problem arises in Computer-Aided Design and Manufacturing (CAD/CAM), virtual reality, medical imaging, computer graphics, computer animation, and many others. Very often, the preferred approximating surface is polynomial, usually described in parametric form. This leads to the problem of determining suitable parametric values for the data points, the so-called surface parameterization. In real-world settings, data points are generally irregularly sampled and subjected to measurement noise, leading to a very difficult nonlinear continuous optimization problem, unsolvable with standard optimization techniques. This paper solves the parameterization problem for polynomial Bézier surfaces by applying the firefly algorithm, a powerful nature-inspired metaheuristic algorithm introduced recently to address difficult optimization problems. The method has been successfully applied to some illustrative examples of open and closed surfaces, including shapes with singularities. Our results show that the method performs very well, being able to yield the best approximating surface with a high degree of accuracy.
1. Introduction
Obtaining a curve or surface that approximates a given cloud of data points is a classical problem in several scientific and technological domains such as computer-aided design and manufacturing (CAD/CAM), virtual reality, medical imaging, computer graphics, computer animation, and many others. In real-world settings, data points come from real measurements of an existing geometric entity, as it typically happens in the construction of car bodies, ship hulls, airplane fuselage, and other free-form objects [1–8]. This process is also applied in the shoes industry, in archeology (reconstruction of archeological assets), in medicine (computed tomography), and in many other fields. The primary goal is to convert the real data from a physical object into a fully usable digital model, a process called reverse engineering. Such digital models are usually easier and cheaper to modify than their real counterparts, leading to a significant reduction of the costs associated with the processing and manufacturing time of the real goods they represent. Furthermore, due to their inherent digital nature, they become available anytime and anywhere, a very valuable feature in our current digital-world era.
Data points in reverse engineering are usually acquired through laser scanning and other digitizing methods (light digitizers, coordinate measuring machines, CT scanners, and tactile scanners) and are, therefore, subjected to measurement noise, irregular sampling, and other artifacts [7, 9]. Consequently, a good fitting of data is generally based on approximation schemes (where the curve/surface is expected to pass near the data points) rather than on interpolation (where the curve/surface is constrained to pass through all input data points). 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 irregularly sampled noisy data points.
There are two key components for a good approximation of data points: a proper choice of the approximating function and a suitable parameter tuning. The usual models for data fitting in CAD/CAM and other industrial fields are free-form parametric entities, such as Bézier, B-spline, and NURBS, as they have a great flexibility and can represent well any smooth shape with only a few parameters, thus leading to substantial savings in terms of computer memory and storage capacity [10–17].
In this paper we focus particularly on the case of polynomial Bézier surfaces, a kind of free-form splines very popular in fields such as CAD/CAM and computer graphics. Bézier splines were developed independently in the early 60s by Paul de Casteljau and Pierre Bézier for the CAD systems of the French automotive companies Citröen and Renault, respectively. Mathematically, they are based on the Bernstein polynomials (see Section 4 for details), developed as early as 1912 but whose applicability to engineering design was unknown until the 60s. A Bézier curve is a linear combination of the Bernstein polynomials and vector coefficients called control points. The curve follows approximately the shape of its control polygon (the collection of segments joining the control points), and hence, it reacts to the movement of its control points by following a push-pull effect. This powerful feature was fundamental for the popularization of free-form curves and surfaces for interactive design. The generalization of this idea to surfaces leads to the Bézier surfaces, which are linear combinations of the control points (now arranged in a three-dimensional net) and the so-called tensor-product basis functions (given by the products of all possible combinations of univariate Bernstein polynomials in surface parameters and , resp.).
Although nowadays Bézier splines have been overtaken by the B-splines (developed during the 70s and of which the Bézier splines are a particular case), they played a key role in the current development of computer design. In addition to their historical value, they are still widely used today for different purposes, such as computer fonts (e.g., TrueType fonts, PostScript), computer animation (for simple movements of objects in programs such as Adobe Flash), and computer design (Adobe Photoshop, Corel Draw, Adobe illustrator). The reader is referred to [18–20] for further details about the subject. See also [21] for a nice historical approach written by some of the most prominent figures in the field.
Best approximation methods make commonly use of least-squares techniques [1, 8, 10, 13, 14, 22–28], where the goal is to obtain the relevant parameters of the polynomial approximating surface that fits the data points better in the least-squares sense. This problem is far from being trivial: because the surface is parametric, we are confronted with the problem of obtaining a suitable parameterization of the data points [18, 20]. As remarked in [29], the selection of an appropriate parameterization is essential for a good fitting. Unfortunately, it also becomes a very hard problem, specially for the cases of irregularly sampled noisy data points. In fact, it is well known that it leads to a very difficult overdetermined continuous nonlinear 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. Finally, it is usually a multimodal problem as well, because of the potential existence of several (global or local) optima of the objective function.
In this context, the present paper describes a new method to solve this challenging parameterization problem for free-form polynomial Bézier surfaces. Our method applies a powerful nature-inspired metaheuristic algorithm, called firefly algorithm, introduced recently by Professor Yang (Cambridge University) to solve difficult 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 an optimal approximating Bézier surface to a given set of noisy data points, provided that an adequate representation of the problem and a proper selection of the parameters are carried out. To check the performance of our approach, it has been applied to some illustrative examples of open and closed surfaces, including shapes with singularities. Our results show that the method performs very well, being able to yield the best approximating surface with a high degree of accuracy.
The structure of this paper is as follows: in Section 2 the previous work in the field is briefly reported. Then, the fundamentals and main ideas of the firefly algorithm, the method used in this paper, are briefly explained in Section 3. Our proposed firefly-based method for data fitting with Bézier surfaces is described in Section 4. The section begins with the description of the problem to be solved. Then, the application of the firefly algorithm to solve it is explained in detail. Some illustrative examples of its application to open and closed surfaces, including shapes with singularities, along with some implementation details are reported in Section 5. The paper closes with the main conclusions of this contribution and our plans for future work in the field.
2. Previous Work
The problem of data fitting through free-form parametric surfaces has been the subject of research for many years [1, 20, 30–35]. One of the most important problems regarding this issue is the surface parameterization, that is, the computation of suitable parametric values for the fitting surface to data points. In many practical situations, it is advisable to obtain a parameterization as similar as possible to the arc-length parameterization. The ultimate reason for this is that a constant step on the parametric domain automatically translates into a constant distance along an arc-length parameterized curve on the surface. In other words, for constant parameter intervals, the curve on the surface exhibits a point spacing that is as uniform as possible. Therefore, this parameterization is very convenient for surface interrogation issues, such as surface intersections or measuring distances on a surface [36, 37]. For instance, it has been traditionally applied in metrology for design and manufacturing, to collect measurement data from industrial parts of the designed and manufactured products. Many other industrial operations also require a uniform parameterization. For example, in computer controlled milling operations, the curve path followed by the milling machine must be parameterized such that the cutter neither speeds up nor slows down along the path [9]. This property is only guaranteed when the curve path is parameterized with the arc-length parameterization. Consequently, this has been the preferred and most classical choice for surface parameterization.
Some recent papers have shown that the application of Artificial Intelligence techniques can achieve remarkable results regarding this parameterization problem [2, 5, 6, 38–40]. Most of these methods rely on some kind of neural networks, either standard neural networks [38], Kohonen’s SOM (Self-Organizing Maps) nets [29, 39], or the Bernstein Basis Function (BBF) network [40]. In the case of surfaces, the network is used exclusively to order the data and create a grid of control vertices with quadrilateral topology [39]. After this preprocessing step, any standard surface reconstruction method (such as those referenced in the bibliography) has to be applied. In some other cases, the neural network approach is combined with partial differential equations [29] or other approaches. The generalization to functional networks (an extension of neural networks where the weights are replaced by functions) is also analyzed in [2, 5, 6, 41].
Due to their good behavior for complex optimization problems involving ambiguous and noisy data, there has recently been an increasing interest in applying nature-inspired optimization techniques (such as metaheuristics and evolutionary methods) to this problem. However, there are still few works reported in the literature. A previous paper in [42] describes the application of genetic algorithms and functional networks yielding pretty good results for both curves and surfaces. 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 are described in [4, 10] for particle swarm optimization (PSO), [3, 27, 28, 43] for genetic algorithms (GA), [44, 45] for artificial immune systems, [46] for estimation of distribution algorithms, and [11] for hybrid GA-PSO techniques. The method used in this paper also belongs to this category, as described in next section.
3. The Firefly Algorithm
The firefly algorithm is a nature-inspired metaheuristic algorithm introduced in 2008 by Yang to solve optimization problems [47, 48] (see also [49] for a recent modified version of this algorithm). 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. The reader is kindly referred to [50] for a comprehensive review of the firefly algorithm and other nature-inspired metaheuristic approaches. See also [51] for a gentle introduction to metaheuristic applications in engineering optimization.
In the firefly algorithm, there are three particular idealized rules, which are based on some of the major flashing characteristics of real fireflies [47]. They are(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: where is the distance defined as in (1), 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 .
The method described in previous paragraphs corresponds to the original version of the firefly algorithm (FFA), as originally developed by its inventor. Since then, many different modifications and improvements on the original version have been developed, including the discrete FFA, multiobjective FFA, chaotic FFA, parallel FFA, elitist FFA, Lagrangian FFA, and many others, including its hybridization with other techniques. The interested reader is referred to the nice paper in [52] for a comprehensive, updated review and taxonomic classification of the firefly algorithms and all its variants and applications.
4. The Proposed Method
A free-form polynomial parametric surface is defined as [18, 19]: where are vector coefficients in (usually referred to as the control points as they roughly control the shape of the surface), are the tensor-product functions obtained from two sets of basis functions (or blending functions) , and , and are the surface parameters, usually defined on a bounded rectangular domain . Note that in this paper vectors are denoted in bold.
In this work we will focus on the particular case of free-form polynomial Bézier surfaces. In this case, (4) becomes where the blending functions are the Bernstein polynomials of index and degree , given by where and the surface parameters , are defined on the unit square . Note that, by convention, .
Let us suppose now that we are given a set of data points in an -dimensional space (usually or ). Our goal is to obtain the free-form polynomial Bézier surface that fits the data points better in the discrete least-squares sense. To do so, we have to compute the control points of the approximating surface by minimizing the least-squares error, , defined as the sum of squares of the residuals:
In the case of irregularly sampled data points , our method will work in a similar way by simply replacing the previous expression (8) by
The least-squares minimization of either (8) or (9) leads to the system of equations: where corresponds to the vectorization of the set of data points (alternatively, ), corresponds to the vectorization of the set of control points , and is a matrix given by , with , , for any , and represents the tensor product of vectors. The indices in (10) vary in the ranges of values indicated throughout the section.
The algebraic solution of (10) is given by, , where denotes the Moore-Penrose pseudoinverse of . Due to the fact that the blending functions are nonlinear in and , the least-squares minimization of the errors is a strongly nonlinear problem, with a large number of unknowns for large sets of data points. Our strategy for solving the problem consists of applying the firefly algorithm to determine suitable parameter values for the least-squares minimization of functional according to either (8) or (9). However, in order to do it, some previous steps must be carefully carried out.(1)First of all, we need an adequate representation of the unknowns of the problem. Because of the tensor-product structure of the free-form Bézier surfaces, the fireflies in our method can be encoded as either strings of two sorted real-coded vectors on the interval of length and , respectively, for organized data points, or as sorted real-coded vectors of length for the case of irregularly sampled data points. All fireflies are initialized with sorted uniformly distributed random numbers on the coordinate parametric domain.(2)The objective function corresponds to the evaluation of the least-squares function given by either (8) or (9). Since this error function does not consider the number of data points, we also compute the RMSE (root-mean squared error), given by for (8) or, alternatively by: for (9) and report our results by using these error criteria.(3)We also need to choose the degree of the approximating surface, which in turn depends on the number of control points. This value is chosen according to the complexity of the shape of the underlying function of data. In general, a small amount of control points is needed for simple, smooth shapes, while a large number of control points must be selected for complicated, twisted, or irregular shapes. Since this number is unknown a priori, it is advisable to start with a low number of control points for each parametric coordinate and increase it until the error reaches values below a prescribed threshold, which generally depends on both the underlying surface and the application domain.(4)Regarding the firefly algorithm, some control parameters should be set up. As usual when working with metaheuristic techniques, the choice of suitable control parameters is very important as it determines the performance of the method at large extent. It is also challenging, because it is strongly problem dependent. In this paper, our choice is based on a large collection of empirical results. These control parameters are(a)the number of fireflies, : this value is set up to fireflies 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;(b)the number of iterations, : this number is another parameter of the method that has to be determined in order to run the algorithm until the convergence of the minimization of the error is achieved. 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 simulations, we found that = 10 is a suitable value, as larger values for this parameter does not improve our results;(c)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;(d)the absorption coefficient, : it is set up to in this paper, as this value provides a quick convergence of the algorithm to the optimal solution;(e)the potential coefficient, : although any positive value can be used for this parameter, the light intensity varies according to the inverse square law. Therefore, we choose accordingly;(f)the randomization parameter, . This parameter varies on the interval and allows us to determine the degree of randomization introduced in the algorithm. This stochastic component is necessary in order to allow new solutions appear and avoid getting stuck in a local minimum. However, larger values introduce large perturbations on the evolution of the firefly and, therefore, delay convergence to the global optima. Consequently, it is advisable to select values in between. In this work, we take .
After the selection of those parameters, the firefly algorithm is performed iteratively for the given number of iterations. To remove the stochastic effects and avoid premature convergence, 20 independent executions have been carried out for each choice of the surface degree. Then, the firefly with the best (i.e., minimum) fitness value is selected as the best solution to the problem.
5. Experimental Results
To check the performance of our method described previously, it has been tested with a large collection of examples with excellent results in all cases. To keep the paper at manageable size, in this section we consider only three of them. They have been primarily chosen to reflect the diversity of situations to which the method can be applied. The examples correspond to both open and closed surfaces, including shapes with singularities. As the reader will see, they clearly show the good performance of our approach.
Examples in this paper are shown in Figures 1, 2, and 3. For each example, two different pictures are displayed: on the left, we show the original cloud of input data points, represented as small red points; on the right, the best approximating Bézier surface, as obtained with our firefly-based method, is displayed. Our input consists of sets of irregularly sampled data points (this fact can readily be seen from simple visual inspection of the point clouds on the left), which are also affected by measurement noise of low to medium intensity (signal-to-noise ratio of 15 : 1, 25 : 1, and 10 : 1, resp.). In all examples, no information about the data points parameterization is available at all. In fact, no information about the structure and properties of the underlying surface of data is either assumed or known beyond the data points.
(a)
(b)
(a)
(b)
(a)
(b)
Table 1 summarizes the main results of our computer simulations. The different examples are arranged in rows. For each example, the following data are arranged in columns: number of data points, error value (according to (8) and (9)), the maximum of the error (denoted by and that provides a useful upper bound for that error), and RMSE error value (according to (11) and (12)). The error values are reported for each coordinate in all cases.
First observation is that, although our data points are irregularly sampled and affected by noise, the method yields very good fitting results in all cases. The RMSE is of order in all cases, while the order of the least-squares error is within the range 10−3–10−2 and so is its maximum. Furthermore, these very small fitting errors are obtained for surfaces that are more complicated than it may seem at first sight. For instance, the surfaces of the first and third examples are apparently simple, flat, and height-map surfaces. However, a careful observation reveals that they oscillate several times, and hence, they exhibit a rich variety of hills and valleys, which have been highlighted by using an illumination model for the sake of clarity. On the other hand, the second surface is a closed surface with a strong singularity at its uppermost part, where many data points concentrate in a very small volume. This is usually a very challenging feature for free-form parametric surfaces, which typically tend to distribute the control points by following a rectangular topology. Clearly, such a distribution is not adequate for this surface. To our delight, the proposed method identifies this situation automatically and rearranges the control points by itself to adapt to the underlying structure of data points. In opinion of the authors, this is a striking and very remarkable feature of this method and shows its ability to capture the real behavior of data points even under unfavorable conditions.
To summarize, a visual inspection of the three figures clearly shows that our method yields a very good approximating surface to data points in all cases. This fact is validated by the numerical results reported in Table 1, which confirm the good behavior of the method. From these examples and many other not reported here for the sake of brevity, we conclude that the presented method performs very well, with remarkable capability to provide a satisfactory, accurate solution to our parameterization problem with polynomial Bézier surfaces.
Regarding the implementation 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 for Windows 8 operating system.
6. Conclusions and Future Work
This paper introduces a new method to address the surface parameterization problem, that is, to compute a suitable parameterization of a set of data points in order to construct the free-form parametric surface approximating such data points better in the least-squares sense. This is a challenging problem that appears recurrently in reverse engineering for computer design and manufacturing and in many other industrial fields. Very often, data points in real-world settings are irregularly sampled and subjected to measurement noise, leading to a very difficult nonlinear continuous optimization problem, which cannot be solved by using standard optimization techniques. To overcome this limitation, this paper proposes a new method based on a powerful nature-inspired metaheuristic algorithm called firefly algorithm, introduced recently to solve difficult optimization problems. The method has been successfully applied to solve the parameterization problem for polynomial Bézier surfaces. The paper discusses the main issues in this problem, such as the solution representation and the selection of suitable control parameters. To check the performance of our approach, it has been applied to some illustrative examples of open and closed surfaces, including shapes with singularities. Our results show that the method performs very well, being able to yield the best approximating surface with a high degree of accuracy.
As mentioned in Section 3, the original firefly algorithm has been improved and modified in many different ways. Some of its variants have shown to be more efficient than the original version, meaning that the presented approach can arguably be improved with new, optimized features for better performance. An illustrative example is given by a very recent version called memetic self-adaptive firefly algorithm [53], whose new capabilities (the use of self-adaptation strategies on the control parameters, a new population model based on elitism, and the hybridization with a local search heuristics) improve the original firefly algorithm significantly. The application of many of these variants to our parameterization problem along with a comparative analysis of their performance is part of our future work. We are also interested to extend this method to other families of surfaces, such as the B-splines and NURBS, where the existence of additional parameters (such as knots and weights) can modify our procedure significantly. The application of this method to some interesting real-world problems in industrial settings is also part of our plans for future work.
Conflict of Interests
The authors of this paper have no current or past, direct or indirect financial relationship with any commercial identity mentioned in this paper that might lead to any conflict of interests. They are solely mentioned here for scientific purposes.
Acknowledgments
This research has been kindly supported by the Computer Science National Program of the Spanish Ministry of Economy and Competitiveness, Project Reference 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 due to the Editor and the anonymous reviewers for their useful comments and suggestions that allowed us to improve the final version of this paper.