Applied Mathematics
Vol.05 No.21(2014), Article ID:51992,10 pages

Energy-Minimizing Curve Fitting for High-Order Surface Mesh Generation

Karsten Bock, Jörg Stiller

Institute of Fluid Mechanics (ISM), Technische Universität Dresden, Dresden, Germany


Copyright © 2014 by authors and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

Received 25 September 2014; revised 20 October 2014; accepted 6 November 2014


We investigate different techniques for fitting Bézier curves to surfaces in context of high-order curvilinear mesh generation. Starting from distance-based least-squares fitting we develop an incremental algorithm, which incorporates approximations of stretch and bending energy. In the process, the algorithm reduces the energy weight in favor of accuracy, leading to an optimized set of sampling points. This energy-minimizing fitting strategy is applied to analytically defined as well as triangulated surfaces. The results confirm that the proposed method straightens and shor- tens the curves efficiently. Moreover the method preserves the accuracy and convergence behavior of distance-based fitting. Preliminary application to surface mesh generation shows a remarkable improvement of patch quality in high curvature regions.


Curvilinear Mesh Generation, High-Order Methods, Bézier Curves, Curve Fitting, Energy Minimization

1. Introduction

The present work is motivated by curvilinear mesh generation for high-order numerical methods such as spectral and hp element methods [1] [2] or discontinuous Galerkin methods [3] [4] . Exploiting the superior convergence properties of these methods for problems of practical interest requires an accurate and well behaved piecewise polynomial representation of complex domains including their boundaries. Unfortunately, contemporary state- of-the-art mesh generators are tailored to low-order discretization methods and, hence, provide only piecewise linear or quadratic meshes. Consequently, various approaches were developed for converting straight meshes into curvilinear ones [5] . Recently, considerable effort has been dedicated to assure the validity and to improve the quality of the curved mesh [6] - [8] . Most of this work deals with situations where the mesh spacing is smaller

or of the same size as the radius of curvature, often in conjunction with a moderate polynomial degree. To the best of our knowledge, the generation of suitable meshes with higher order (say) and spacings exceeding the radius of curvature still remains a challenge.

As a first step of curvilinear mesh generation we consider the construction of polynomial curves from a given straight-sided surface mesh. This problem can be regarded as a special case of curve fitting, which, in turn, is a well established research area in computer aided geometric design [9] [10] . Typically the sought curve is not fitted to the surface as such, but to a set of samples extracted from the latter. A widespread approach is to use splines in tension or smoothing splines, which attain a fair shape by minimizing a certain energy functional related to stretching, bending, twist, or a combination thereof [9] . Veltkamp and Wesselink [11] explored these and a variety of related energy functionals in fitting B-splines to a set of given points in or. Higher-order curves such as B-splines are often computed by minimizing an error functional augmented by a suitable energy measure for regularization [12] . Unfortunately, the energy functional tends to act diametrically to the error functional and thus prevents convergence to the exact surface with increasing polynomial order. As an alternative Alhanaty and Bercovier [13] proposed the use of optimal control methods for constructing interpolating splines with minimal energy. Apart from the advantage of guaranteeing optimality, this approach implies that the number of samples does not exceed the degree of freedom available for fitting, which may prove too restrictive for many applications. Flöry and Hofer [14] advocated an incremental strategy for fitting curves on manifolds using a weighted average of error and energy functionals as the objective function. In course of their iterative procedure the weight of energy is successively reduced, thus improving the accuracy of the final curve. It should be remarked, however, that this method is not intended for generating energy-minimized curves. Indeed, the curve energy is essentially determined by the chosen sample points, which remain fixed throughout the fitting procedure. A direct approach to energy-minimizing splines was developed by Hofer and Pottmann [15] [16] . The key idea consists in translating the sample points tangentially on the given manifold until the imposed energy functional attains a minimum. This process is embedded in an iterative procedure which constrains the movement of samples to a trust region. More precisely, the method does not minimize the curve energy, but a finite difference approximation of the latter based on the sample points. The optimality of the resulting curve will therefore depend on the sampling density. Finally we remark that curvature minimizing smoothing offers a powerful alternative to fitting worth considering especially with noisy surface data, see e.g. [17] and references therein.

In this paper we investigate techniques for fitting Bézier curves to surfaces for application with high-order mesh generation. Starting from squared distance minimization we develop an incremental algorithm leading to an accurate, energy-minimizing method that combines the ideas put forward in [14] [15] . The paper is organized as follows: In Section 2 we first revisit least-squares fitting and then elaborate the energy-minimizing curve fitting procedure. Section 3 provides a comparison of both methods covering analytically defined smooth surfaces as well as scattered surface data. Section 4 concludes the paper.

2. Surface Curve Construction

In spectral or hp element methods, each element face coinciding with the domain boundary constitutes a polynomial surface patch. Therefore, the construction of well behaved polynomial patches represents a natural building block in curvilinear mesh generation. Starting from a straight-sided initial mesh, the curvilinear mesh is often built in a hierarchical process consisting of the following steps: i) construction of boundary curves representing the edges of the boundary faces; ii) generation of patches defining the boundary faces; and iii) creation of curved volume elements. Here we focus on the first step, the construction of high-order polynomial boundary curves. Adopting the Bézier form, a curve of order is expressed as


where are the Bernstein polynomials and the corresponding control points. For convenience we assume that the vertices of the initial mesh and, hence, the start and end points of the curves are fixed and given such that they fit to the boundary surface.

2.1. Distance-Based Curve Fitting

Curve fitting is used to construct a boundary edge between vertices and. For this purpose we perform a least squares fitting to a set of sampling points generated from the initial, linear edge.

The procedure starts with the selection of samples in parameter space. Next, a corresponding point distribution is generated on the straight edge by means of linear interpolation. Projecting these points to the boundary surface yields the sampling points


where, ideally, the operator represents the normal projection of a given point to a surface point. In practice the projection is realized by means of an iterative procedure based on the approximate normal vector. In case of scattered data, we use a fine triangulation combined with a smoothing interpolation as the surface definition. Fitting is performed by minimizing the average squared distance


where denotes set of the interior control points of the curve, the set samples in parameter space and the corresponding set of surface points. Since the number of samples typically exceeds the degrees of freedom, the minimization problem is solved using the least-squares method, which yields the interior control points of the curve.

2.2. Energy-Minimizing Curve Fitting

With coarse meshes purely distance-based fitting can lead to severe undulations in regions of high curvature. As a remedy one may look for curves that minimize a certain energy functional. Here we consider the norms of the first and second derivative of:



The norm of the first derivative (4) is related to the elastic stretch energy of a string [11] . A surface curve which minimizes represents the shortest path between between the two endpoints and is called a geodesic. Geodesics are regularly applied in computer vision, image processing and motion design [16] . It is worth noting that the functional also improves the curve towards an isometric parametrization. The norm of the second derivative (5) corresponds to bending or strain energy. Minimizing is a fundamental ingredient in cubic spline interpolation. The combination of both energies leads to the concept of splines in tension (see, e.g. [9] ). We use and for augmenting the error functional (3), leading to the objective function


where overbars indicate the normalization which has been introduced to compensate possible differences in the order of magnitude between the individual terms. In particular we define




where the superscript “0” refers to the curve obtained from distance-based fitting and “lin” to the straight edge. Note that equals the squared length of the latter, whereas is identical to zero. As these two cases are close to the extrema that are to be expected in the fitting process, the normalized distance and error functionals vary typically between zero and one. The parameter represents the relative weight of the energy functionals compared to the error functional, while the value of determines the blend of bending and stretch energy to be used in the former. Hence, applying with is equivalent to distance-based fitting, whereas choosing any and thus taking into account the curve energies leads to smoother and shorter curves but allows for deviation from the surface.

To obtain fair surface curves we employ an incremental approach, which is outlined in Algorithm 1. The basic idea is to start fitting with a high energy weight, which is successively reduced according to a generic shape function, until reaching the minimum value in the final step (line 8). Furthermore, the sampling points are recomputed from the current curve (line 9) before performing the next fitting (line 10). In this manner, curve and surface points move towards each other in course of the process, the former approaching the surface, the latter converging to an energetically improved configuration.

Throughout the present study we used the shape function


which starts with 1 at followed by a gradual transition to 0 at. This choice always recovers the straight edge in the first step. Note that the final curve minimizes the mean quadratic error as well as the energy, since it results from distance-based fitting to optimized sampling points. Various other choices for the shape function are possible, but have not been explored yet.

3. Results

In the following we study the performance of the energy-minimizing fitting method in two different cases. As the first case we consider a coarse, but nearly uniform triangulation of an explicitly defined screw surface. In the second case the “exact” surface is defined as a patchwork of cubic triangles based on a fine mesh derived from CT scans of a rabbit aorta. For assessment we use the error evaluated over all mesh curves


and the curve energy norms


Algorithm 1.Incremental energy-minimizing curve fitting starting from linear edge.

where again corresponds to stretch energy and to bending energy. The integrals in Equation (10) were computed by means of Gauss quadrature using points.

3.1. Analytically Defined Smooth Surface

As a first test case we consider the screw surface defined by the analytical expression


The projection of a given point onto the surface is realized by means of the iterative procedure


starting with. We remark that this scheme results from a linearised solution of exact condition, where is the surface normal vector. Figure 1 shows the exact surface in the interval along with a coarse uniform triangulation. The mesh consists of 32 triangles and has an average spacing of times the minimum curvature radius,.

The energy-minimizing fitting procedure was examined over a wide range of polynomial degrees, ranging from to 20. As an example we consider the case. Figure 2 illustrates the influence of the energy blending parameter and the number of fitting steps on the error and the energy norms and. The graphs correspond to three different energy compositions in the objective function: implies using only stretch energy, only bending energy and an equal mixture of both. Note that the plot does not show stepwise progress, but only the characteristics of the final curve. The case is equivalent to single-step distance-based fitting. Increasing the number of fitting steps reduces the curve energy rapidly until approximately. In the present example, this improvement is accompanied by a slight increase in the error, which lessens with growing. We remark, however, that with finer meshes the energy-minimizing fitting succeeded in reducing the error simultaneously with the curve energy.

Figure 2 indicates that the incremental method improves both energies, even if one is excluded from the objective function by the particular choice of. However, with higher the excluded energy may increase slightly in favor of a further error reduction. In the present example such behavior is observed with for and for, whereas retains both energies close to the minimum.

Figure 3 shows the effect of different parameter choices for an edge crossing the ridge of the screw. In this

Figure 1. Exact definition (left) and a coarse linear grid (right) of the screw surface example. The mesh contains 32 triangles, which results in an average mesh spacing of 10.85 times the minimum curvature radius of the geometry.

Figure 2. error and energy norms for energy-minimizing fitting with order and different blending parameters as a function of the number of fitting steps.

Figure 3. Comparison of different fitting methods for “Curve 54” crossing the ridge of the screw: dotted magenta line distance-based, dashed red energy- minimizing with, solid blue and dash-dotted green. Energy-minimizing fitting was performed with steps.

case, mere distance-based fitting yields a meandering curve, whereas energy-minimizing fitting straightens the path and removes undulations regardless of the chosen energy composition. As expected, the balanced energy mix, , yields a curve nestling in between the extreme cases of minimal stretch, , and minimal strain,. Yet these curves follow a rather similar path.

Figure 4 illustrates the performance of the energy-minimizing fitting procedure with for the edge depicted in Figure 3. Starting with the straight edge, the normalized error functional assumes its maximum, while the normalized energies and vanish after the first step. As a result of tapering both energies increase, thus allowing for error reduction. At the end of the process, the normalized error approaches zero and thus recovers the accuracy of single-step distance-based fitting. The energy functionals finally approach and, respectively. The latter value indicates a significant strain reduction compared to the reference case. Note that is related to the curve length and thus dominated by surface curvature, which explains the comparably lower reduction.

Finally we look at the convergence behavior with respect to the fitting order. Figure 5 depicts the errors of distance-based fitting and energy-minimizing fitting with and. The comparison shows only a negligible difference between the two graphs. In both cases the slope attains an asymptotic regime which indicates a linear dependence of the error exponent on the fitting order, i.e.,. Thus we conclude that the proposed energy-minimizing fitting procedure preserves the spectral convergence of the original, distance-based method.

3.2. Scattered Surface Data

Scanning methods such as computed tomography (CT) or magnetic resonance imaging (MRI) provide scattered data that can be processed to give triangulated volume and surface representations of the investigated object. Here we consider a partition of a rabbit’s aortic arch, given as a fine mesh consisting of 24,644 triangles (Figure 6).

Figure 4. Evolution of the normalized error and energy functional for Curve 54 in course the energy-minimizing fitting procedure with using 100 steps.

Figure 5. Convergence of distance-based and energy-minimizing fitting with respect to the polynomial degree of curves.

Figure 6. Fine mesh (blue) of a rabbit aorta and coarse mesh used for fitting (red). Fine mesh courtesy of Spencer Sherwin, Imperial College London.

This representation was enhanced in two steps, yielding the “exact” surface: First, computing the vertex normals using the method of Max [18] . Second, constructing a cubic interpolant based on the point-normal vertex data in terms of the PN Triangles proposed by Vlachos et al. [19] . The projection of a point onto the surface is defined as the shortest projection along the Phong normal [20] , which is based on the linear mesh. For details we refer to [21] .

To evaluate the curve fitting methods we generated a curvature-dependent coarse mesh comprising 532 triangles (Figure 6). The ratio between the local mesh spacing and the radius of curvature was limited in order to allow for reasonable accuracy with moderate fitting order.

Starting from the linear mesh we constructed curves of order. Figure 7 presents the curves after 50 fitting steps in comparison with the distance-based fitted curves in the high curvature bifurcation region. For clarity, the “exact” surface is also shown. In accordance with the previous test case the energy-minimizing fitting method yields straighter and shorter curves. In low curvature regions both methods result in nearly straight curves. This is not surprising, since a straight line is energetically optimal. Therefore, the computational cost can be reduced by applying energy-minimizing fitting only to those curves, which can be expected to profit from it. We explored a decision criterion based on the relative change in length


where is the length of the curve resulting from distance-based fitting and the length of the straight edge. Only when exceeds a certain threshold energy-minimizing fitting is performed. Table 1 compiles the results of distance-based fitting with energy-minimizing fitting using either no threshold (equivalent to) or. Note that the energy norm is adjusted by the value, which represents a mesh-dependent lower bound corresponding to the straight edges.

Remarkably, the energy-minimizing approach not only succeeds in reducing the energy norms, but also improves the approximation accuracy. Assuming a threshold of energy-minimized fitting is applied only to about ten percent of edges, which saves almost 90 percent of the computational cost. It is important to note that the improvements in energy and accuracy are not affected by this measure.

In addition we constructed triangular Bézier patches in two steps [21] : First, we inject the previously generated boundary curves. Second, we determine the inner control points by least-squares minimization of a distance functional. To assess the patch quality we consider the distortion measure, which is based on the Jacobian of the mapping from parameter to physical space [5] . This measure has an upper bound of corresponding to an isometric mapping, which is the preferred case. Smaller values indicate distorted patches and with the mapping becomes invalid. As can be seen from Table 1, the energy-minimized curves lead to improved patches without further optimization. The quality measure increases by a factor of, yielding for all patches, which represents an acceptable quality for a high curvature geometry.

4. Conclusion

We investigated different techniques for fitting Bézier curves to surfaces as a first step of curvilinear mesh generation for high-order discretization methods. As a starting point we examined a distance-based least-squares fitting method. This method achieves a high accuracy, but tends to produce distorted curves where the mesh spacing is large compared to the radius of curvature. As remedy, we included approximations of stretch and bending energy into an incremental algorithm, resulting in an energy-minimizing fitting method. Both approaches were evaluated using two examples: an analytically defined screw surface and a surface triangulation of a rabbit aorta. The results confirm that the energy-minimizing method straightens and shortens the curves efficiently. Moreover the method preserves the accuracy and convergence behavior of distance-based fitting. In accordance with previous work (see e.g. [9] [11] [16] ), our study indicates that combining stretch and bending energy yields better results than using only one of those. Additionally, we analyzed the influence of the curve fitting method on

Figure 7. Comparison of fitting methods in the bifurcation of the rabbit aorta: dotted magenta lines distance-based, solid blue energy-minimizing with and steps.

Table 1. Results obtained with different fitting approaches for the rabbit aorta test case with polynomial order. Listed are: the error, the energy norms and and the Jacobian distortion. The offset corresponds to the straight edges. Energy-minimizing fitting was performed with and.

patch construction. This investigation shows a clear improvement in patch quality when using energy-minimized curves. Nonetheless distortion remains an issue in the patch interior. Therefore, future work should address the extension of energy-minimizing approach to surface patch fitting.


The authors gratefully acknowledge the funding of this project by the German Research Foundation (DFG, STI 157/4-1).

Cite this paper

Karsten Bock,Jörg Stiller, (2014) Energy-Minimizing Curve Fitting for High-Order Surface Mesh Generation. Applied Mathematics,05,3318-3327. doi: 10.4236/am.2014.521309


  1. 1. Karniadakis, G.E. and Sherwin, S. (2005) Spectral/hp Element Methods for Computational Fluid Dynamics. Numerical Mathematics and Scientific Computation. Oxford University Press, Oxford.

  2. 2. Deville, M., Fischer, P.F. and Mund, E.H. (2002) High-Order Methods for Incompressible Fluid Flow. Cambridge University Press, Cambridge.

  3. 3. Cockburn, B.B., Karniadakis, G.E. and Shu, C.-W. (2000) Discontinuous Galerkin Methods: Theory, Computation and Applications. Lecture Notes in Computational Science and Engineering. Springer, Berlin.

  4. 4. Hesthaven, J.S. and Warburton, T. (2008) Nodal Discontinuous Galerkin Methods/Algorithms, Analysis, and Applications. Springer, Berlin.

  5. 5. Dey, S., O’Bara, R.M. and Shephard. M.S. (1999) Curvilinear Mesh Generation in 3d. Proceedings of the Eighth International Meshing Roundtable, John Wiley & Sons, Hoboken, 407-417.

  6. 6. Sherwin, S.J. and Peiró, J. (2002) Mesh Generation in Curvilinear Domains Using Highorder Elements. International Journal for Numerical Methods in Engineering, 53, 207-223.

  7. 7. Persson, P.-O. and Peraire, J. (2009) Curved Mesh Generation and Mesh Refinement Using Lagrangian Solid Mechanics. Proceedings of the 47th AIAA Aerospace Sciences Meeting and Exhibit.

  8. 8. Toulorge, T., Geuzaine, C., Remacle, J.-F. and Lambrechts, J. (2013) Robust Untangling of Curvilinear Meshes. Journal of Computational Physics, 254, 8-26.

  9. 9. Hoschek, J. and Lasser, D. (1996) Fundamentals of Computer Aided Geometric Design. Wellesley, Massachusetts.

  10. 10. Farin, G. (2002) Curves and Surfaces for CAGD—A Practical Guide. 5th Edition, Academic Press, Waltham.

  11. 11. Veltkamp, R.C. and Wesselink, W. (1995) Modeling 3D Curves of Minimal Energy. Blackwell Science Ltd., Computer Graphics Forum, 14, 97-110.

  12. 12. Wang, W., Pottmann, H. and Liu, Y. (2006) Fitting B-Spline Curves to Point Clouds by Curvature-Based Squared Distance Minimization. ACM Transactions on Graphics, 25, 214-238.

  13. 13. Alhanaty, M. and Bercovier, M. (2001) Curve and Surface Fitting and Design by Optimal Control Methods. Computer-Aided Design, 33, 167-182.

  14. 14. Flöry, S. and Hofer, M. (2008) Constrained Curve Fitting on Manifolds. Computer-Aided Design, 40, 25-34.

  15. 15. Hofer, M. and Pottmann, H. (2004) Energy-Minimizing Splines in Manifolds. ACM Transactions on Graphics, 23, 284-293.

  16. 16. Hofer, M. (2007) Constrained Optimization with Energy-Minimizing Curves and Curve Networks: A Survey. Proceedings of the 23rd Spring Conference on Computer Graphics, Comenius University, Bratislava, 27-35.

  17. 17. Lawonn, K., Gasteiger, R., R?ssl, C. and Preim, B. (2014) Adaptive and Robust Curve Smoothing on Surface Meshes. Computers & Graphics, 40, 22-35.

  18. 18. Max, N. (1999) Weights for Computing Vertex Normals from Facet Normals. Journal of Graphics Tools, 4, 1-6.

  19. 19. Vlachos, A., Peters, J., Boyd, C. and Mitchell, J.L. (2001) Curved PN Triangles. Proceedings of the 2001 Symposium on Interactive 3D Graphics, I3D’01, New York, 159-166.

  20. 20. Phong, B.T. (1975) Illumination for Computer Generated Pictures. Communications of the ACM, 18, 311-317.

  21. 21. Bock, K. and Stiller, J. (2014) Generation of High-Order Polynomial Patches from Scattered Data. In: Azaïez, M., El Fekih, H. and Hesthaven, J.S., Eds., Spectral and High Order Methods for Partial Differential Equations-ICOSAHOM 2012, Lecture Notes in Computational Science and Engineering, Vol. 95, Springer International Publishing, 157-167.