IGS: an IsoGeometric approach for Smoothing on surfaces
We propose a novel approach for smoothing on surfaces, namely estimating a function starting from noisy and discrete measurements. More precisely, we aim at estimating functions lying on a surface represented by NURBS, which are geometrical representations commonly used in industrial applications. The estimation is based on the minimization of a penalized least-square functional. The latter is equivalent to solve a 4th-order Partial Differential Equation (PDE). In this context, we use Isogeometric Analysis (IGA) for the numerical approximation of such surface PDE, leading to an IsoGeometric Smoothing (IGS) method for fitting data spatially distributed on a surface. Indeed, IGA facilitates encapsulating the exact geometrical representation of the surface in the analysis and also allows the use of at least globally continuous NURBS basis functions for which the 4th-order PDE can be solved using the standard Galerkin method. We show the performance of the proposed IGS method by means of numerical simulations and we apply it to the estimation of the pressure coefficient, and associated aerodynamic force on a winglet of the SOAR space shuttle.
Keywords: Functional Data Analysis; Isogeometric Analysis; Smoothing on Surfaces.
The estimation of a function from a set of noisy data is a very common task, which is often tackled by minimization of a penalized least-square functional, where the penalty involves a differential operator, commonly based on second order derivatives, and occasionally on derivatives of a different order. Classical examples are offered by smoothing splines (see, e.g., Ramsay and Silverman (2005) and references therein) for estimating functions defined over real intervals, thin-plate splines (see, e.g., Wahba (1990), Wood (2006) and references therein), Soap film smoothing (Wood et al., 2008), splines over triangulations (see, e.g., Lai and Schumaker (2007)) for estimating functions defined over regions of , and spherical splines (see, e.g., Wahba (1981); Baramidze et al. (2006); Alfeld et al. (1996)) for estimating functions defined over spheres or spheres-like surfaces. In this context, one of the main challenges in minimizing such penalized least-square functional consists in determining a suitable finite dimensional space representing, at the discrete level, the infinite dimensional space to which the function belongs. In other words, the challenge is to find a finite dimensional problem which is tractable and whose solution is close to the solution in the infinite dimensional space. The same challenge arises when looking for the numerical solution of a PDE. This common goal has been recently exploited to develop statistical tools to deal with spatially distributed data. In this respect, Ramsay (2002) considered planar smoothing optimizing a penalized least-square functional with a regularization term involving the Laplacian, and used the Finite Element Method (FEM) (see, e.g., Quarteroni and Valli (1994)) to solve the estimation problem. Sangalli et al. (2013) generalized the method proposed by Ramsay (2002) to include space-varying covariates and to account for boundary conditions, while Wilhelm (2013) explored a generalized linear version of the method. Azzimonti et al. (2014, 2015) further extended the model of Sangalli et al. (2013) to account for any elliptic differential penalization, not necessarily based on the Laplacian operator. Ettinger et al. (2015) and Dassi et al. (2015) dealt with the case of functions defined over two dimensional Riemannian manifolds, by considering a regularizing term based on the Laplace-Beltrami operator and exploited a conformal parametrization of the manifold to solve an equivalent estimation problem on . Duchamp and Stuetzle (2003) also explored smoothing on complex surfaces using a penalization based on the Laplace-Beltrami operator. A common feature of all these contributions is the use of FEM to estimate the underlying function corresponding to the observed data.
In this paper, we consider the estimation of functions defined over surfaces in starting from a discrete set of noisy observed data in points distributed on such surfaces. More precisely, we refer to sufficiently smooth functions defined over surfaces that can be represented by NURBS (Non-Uniform Rational Basis Splines). Indeed, NURBS are commonly used in Computer Aided Design (CAD) to represent most of the geometries of Engineering and industrial interest. From a more general point of view, the proposed model is a particular case of a Generalized Additive Model (GAM) (Hastie and Tibshirani, 1990) where the smooth component is defined on a surface. Hence, most of the theoretical results of GAMs can be directly applied to this model, e.g., for the quantification of uncertainty.
IGA has been first introduced by Hughes et al. (2005) with the main idea of using the same basis functions to represent the geometry and then to approximate the solution of the PDEs defined in such computational domains. This facilitates encapsulating the exact geometrical representation when performing the analysis of PDEs defined in the computational domain; see Cottrell et al. (2009). In this respect, the isogeometric paradigm facilitates the use of an exact representation of the surface, while most of the current methodologies, as FEM, generally only handle its approximation. This may induce an error on the solution due to the approximation of the geometry and may require complex meshing procedures. As mentioned, to estimate the function over the surface, starting from its discrete and noisy measurements, we minimize a least-square functional where the penalty involves the Laplace-Beltrami operator associated to the surface, analogously to Ettinger et al. (2015), Duchamp and Stuetzle (2003) and Wahba (1981). The estimation problem is tackled directly on the surface and using IGA to numerically solve the associated surface PDE. For this reason, we name the resulting method IsoGeometric Smoothing (IGS). High-order PDEs can be straightforwardly solved by NURBS-based IGA; indeed, globally continuous NURBS basis functions can be easily defined on surfaces for some , where is the polynomial degree. This allows to use standard Galerkin method for numerically solving the PDE. In this respect, Dedè and Quarteroni (2015) studied the use of IGA for 2nd-order PDEs defined on surfaces, Tagliabue et al. (2014) analysed IGA for high-order PDEs, while Bartezzaghi et al. (2015) for high-order PDEs defined on surfaces. All these works showed the efficiency and accuracy of IGA for the spatial approximation of surface PDEs. Partially related to our work, NURBS-based and IGA approaches are used in Beaubier et al. (2014) and in Dufour et al. (2015) as calibration tools.
In this paper, we introduce the geometrical background to define the application framework of IGS. Then, we show how one can minimize a least-square functional involving a penalization term and hence estimate a function starting from a set of discrete and noisy measurements. We assess and compare IGS to the Thin Plate Splines (TPS) (Duchon, 1977; Wahba, 1990) by means of numerical simulations. Finally, as industrial application, we estimate the pressure coefficient and the aerodynamic force acting on the inboard winglet of the space shuttle SOAR, designed by S3, Swiss Space Systems Holding SA. S3 is a Swiss company currently developping, manufacturing, certifying, and operating a launch system for small satellites of weight inferior to 250 kg. A geometry of the SOAR suborbital shuttle for a preliminary study of aerodynamic forces is shown in Figure 1. In this application, the data are pointwise measurements of the pressure coefficient and the associated quantity of interest is the aerodynamic force. We propose a method to estimate functions defined over a NURBS surface starting from scattered noisy observations.
This paper is organized as follows. In Section 2, we briefly recall B-splines and NURBS. In Section 3, we introduce some analytical tools and the IGS method. We show results for two numerical simulations in Section 4: the first one allows to compare IGS to TPS, while the second one corresponds to a complex NURBS surface, for which planar smoothers are not suited. Finally, in Section 5, we show the results of the estimation of the pressure coefficient and aerodynamic force over the winglet of the SOAR space shuttle. Conclusions follow.
2 B-splines and NURBS
NURBS are widely used in CAD for geometrical representation of surfaces (Piegl and Tiller, 1997). We start by defining an open knot vector of degree as a sequence of values , with knots , where the first and the last knots are repeated times. The interior knots can be repeated at most times and, if a knot is repeated times, we say that its multiplicity is . The B-spline basis functions are defined recursively using the Cox-De Boor formula. We summarize some properties of B-spline basis functions of index and degree , for .
The basis function possesses continuous derivatives across a knot of multiplicity and is continuous between the knots.
The support of the basis function is compact and contained in knots spans .
The basis functions are pointwise nonnegative, i.e. for all .
They form a partition of the unity, that is , for all .
Starting from the basis function , the NURBS basis functions are defined as projective transformations of B-spline basis functions. Let be positive weights, then, NURBS basis functions are defined as:
A NURBS curve , with and control points is defined as:
To define surfaces with NURBS, we resort to the tensor product scheme. Given two open knot vectors and , , for , , for , the corresponding univariate basis functions, the control points and positive weights for we define the bivariate NURBS basis functions as:
and a NURBS surface as:
where . For simplicity, we consider henceforth the same polynomial degree along both the parametric directions, for which we rewrite the bivariate basis functions simply as , for where . For an exhaustive description of NURBS, we refer the reader to Piegl and Tiller (1997).
We remember that in IGA, one uses the same basis functions to represent the surface and to approximate the solution of the PDE defined in such computational domain. In general, it is possible to enrich the NURBS basis without changing the geometry, i.e. by preserving the geometrical mapping, with the goal of obtaining a more accurate solution of the PDEs. In this respect, refinement indicates a uniform knot insertion, which adds new basis functions while preserving the geometrical mapping, while refinement refers to an elevation of the polynomial degree of the basis functions, similarly to FEM. In addition, the so-called refinement is peculiar of NURBS and refers to a consecutive order elevation and a knot insertion which allow to increase the polynomial degree and continuity of basis functions. All the refinement procedures are discussed in details in Cottrell et al. (2009), while for an introduction to NURBS in the context of IGA, we refer the interested reader to Hughes et al. (2005).
3 The IGS method
We describe the mathematical framework of the IGS method. First, we introduce the surface differential operators and then we introduce the IGS method for smoothing functions on surfaces.
3.1 Geometrical framework
Following Dedè and Quarteroni (2015), let be an open and bounded parametric domain of finite measure with respect to the topology of . Then, let be a compact, connected, and oriented surface, defined by a NURBS geometrical mapping such that:
We assume that is sufficiently smooth, e.g. . Then, we define the Jacobian of the mapping , denoted by as:
We denote by the -th column of the matrix . The metric tensor of the mapping is represented by:
and we denote with the square root of the determinant of the metric tensor of the mapping :
The metric tensor of the mapping is assumed to be invertible almost everywhere in .
3.2 Differential operators
Let be a smooth function defined over the surface , i.e.:
such that we can define the differential operators on the surface . Then, the projection operator on the tangent plane to the surface at the point is given by:
where is the unit normal vector to the surface at the point and is the identity matrix
The Laplace-Beltrami operator, which is the surface Laplacian operator, is expressed as:
where , for all , is the surface divergence and
is the Hessian matrix of . These operators can be rewritten in the parametric domain using the geometrical mapping (3.1) as:
3.3 Mathematical model
Let us consider points located on for which the observed values are respectively. We assume that:
where is a sufficiently smooth function defined on that we aim at estimating and are independent observational errors of zero mean and constant variance.
Given a positive smoothing parameter , we aim at minimizing the following parameter dependent functional:
where and is the vector of evaluations of the general function at the points . This functional is the same considered by Ettinger et al. (2015) and by Duchamp and Stuetzle (2003) that uses a finite element representation, and by Wahba (1981), that considers functions on spheres and proposes spherical splines. The functional (3.4) involves the surface Laplace-Beltrami operator, which is, roughly speaking, a measure of the curvature of the function related to the surface. The use of the Laplace-Beltrami operator ensures that the regularization is invariant to rigid transformations of the coordinate system, which is desirable both for theoretical and practical reasons. If the observational errors are normally distributed, the functional (3.4) can be viewed as a negative rescaled Gaussian penalized log-likelihood. Hence, in such case, minimizing (3.4) is equivalent to maximizing a penalized log-likelihood. Then, for a given positive parameter , the estimation problem is:
where is a suitable functional space to be defined to ensure the well-posedness of the problem and is the estimate of in the functional space , which should lie at least in , i.e. . Indeed, since , the evaluation of at the points is well defined. Since the problem (3.5) will be later associated to a PDE, some essential boundary conditions (Brezis, 1999) should be specified in relation with the choice of . In the simpler context of functions defined over planar domains, and with the penalizing terms involving linear second order elliptic operators, the well-posedness of problem (3.4) is extensively discussed in Azzimonti et al. (2014) under different kind of boundary conditions (Dirichlet, Neumann, or Robin (see, e.g., Brezis (1999)). In the following, we prove the well-posedness of the estimation problem (3.5) in the particular case of homogeneous boundary conditions using the Lax-Milgram theorem, following Ramsay (2002) and Sangalli et al. (2013). Let be defined as:
where denotes the outward directed unit vector normal to , the boundary of . Then, one can characterize the solution of the minimization problem (3.5) and ensure the existence and the uniqueness of the solution in the case where is assumed to lie in . With this aim, we recall the Lax-Milgram Theorem (see, e.g., Quarteroni and Valli (1994)):
Theorem 3.1 (Lax-Milgram).
Let be a Hilbert space, a continuous and coercive bilinear form and a linear and continuous functional. Then, there exists a unique solution of the following problem:
Moreover, if is symmetric, then is the unique minimizer in of the functional , defined as
Let be two functions and a positive smoothing parameter. Let and be defined as:
where and for some distinct points on . Then, the bilinear form is coercive, continuous and symmetric and is linear and continuous in .
First, we note that , defined as , is equivalent to the norm in . This implies that there exists a constant such that (see, e.g., Quarteroni (2015)). Then, we have:
and so is coercive. We now show the continuity of . Since , there exists a constant such that . We have:
Then, the bilinear form is also continuous and its symmetry is obvious.
Finally, we have:
which proves the continuity of the linear form and concludes the proof. ∎
The functional (3.4) can be rewritten as:
By defining as
we have that
From the definitions of and of (3.7), the functional can be written as:
Thanks to the Lax-Milgram Theorem, it is then sufficient to use Lemma 3.1 to establish the well-posedness of problem (3.5). Moreover, since is symmetric, problem (3.5) is equivalent to problem (3.8). ∎
Setting a priori the essential boundary conditions and following (3.6) may be an inadequate choice in several applications, especially when the behaviour of the function at the boundaries is not known a priori. In such cases, it may be more convenient to consider instead natural boundary conditions, that is:
In the case that the boundary conditions (3.9) are set, we are unable to show the well-posedness of problem (3.8) with . Nevertheless, numerical experience indicates that it still yields a numerically stable problem. We can also note that the usual planar smoothers, such as TPS, also implicitly use natural boundary conditions.
3.4 Numerical approximation: IGS
Let be a finite dimensional approximation of obtained by means of IGS. Let be a basis of a discrete function space of dimension . In the finite dimensional space , problem (3.8) reads:
where and . Since is finite dimensional, problem (3.10) is equivalent to:
where . Let us define the matrix as
and the matrix as Since belongs to , it can be written as a linear combination of the basis functions:
or compactly as
. Then, we have:
Problem (3.11) in matrix form reads as:
where . Then, the explicit form of is given by:
We see from (3.12) that the estimator has the typical form of a penalized least-square estimator. Since we finally get the evaluation of the function in the points as:
The smoothing matrix, that maps the observed data values in the fitted data values , is given by:
The trace of is a measure of the equivalent degrees of freedom of the estimator (see Buja et al. (1989)). If , the number of degrees of freedom is equivalent to the number of basis functions . However, the two notions differ for . While different definitions of equivalent degrees of freedom can be considered, these can be assumed as a consistent measure of the number of degrees of freedom that takes into account the harmonic penalization. In this respect, if the smoothing parameter is strictly positive, the number of equivalent degrees of freedom is smaller than the number of basis functions used in .
We use IGA to solve the minimization problem (3.5), for which we define:
where are the NURBS basis functions used to build , eventually after the application of some , or refinement procedure, as described in Section 2. is the number of basis functions, which is the dimension of . NURBS allow to define basis functions which are globally continuous on . As consequence, one can approximate problem (3.5) with the standard Galerkin method, since ; see Bartezzaghi et al. (2015) and Tagliabue et al. (2014). In this manner, we obtain a method, which we name IGS, allowing to perform smoothing on surfaces by means of NURBS-based IGA. This also allows encapsulating the original description of the geometry of the surface in the analysis.
The smoothing parameter may be chosen by minimization of a generalized cross-validation criterion (GCV), defined as:
see Craven and Wahba (1978). Here, we use as a measure of the equivalent degrees of freedom (EDF) of the model (Buja et al., 1989). In order to solve the optimization problem corresponding to GCV minimization, we use a BFGS quasi-Newton method (see, e.g., Nocedal and Wright (1999)) with a sufficiently small tolerance. One can observe that the computation of the GCV criterion involves the inversion of the matrix of size . In our implementation we use a direct method to compute the matrix because of the moderate size of this matrix. Methods based on matrix decompositions can improve both the efficiency and the stability of the optimization procedure used for the GCV criterion, (see, e.g., Wood, 2006, pp. 178–181). Other methods for the choice of the smoothing parameter are also available, e.g. the restricted maximum likelihood estimation (Wood, 2011).
We remark that our model only considers a deterministic location of the measurement points, according to (3.3). While it is easy to account for random location of points at the implementation level of the IGS method, the uncertainty quantification in this setting is not straightforward.
3.5 Distributional properties and quantification of uncertainty
Here, we denote by and the expectation and the variance of respectively. Moreover, let be the constant variance of the noise introduced in (3.3). For a given smoothing parameter , the estimate is a linear transformation of the observations , as shown in (3.12). Moreover, we have and . Then, from (3.12), we get:
Similarly, we can directly express the variance as:
where we used the fact that the matrix is symmetric and hence . In particular, under the assumption of normality of the errors, we have:
Then, we can also express the expectation and the variance of the fitted values explicitly as:
respectively, since the smoothing matrix is symmetric. In practice, the error variance must be estimated from the data. Following Hastie and Tibshirani (1990) we can estimate by:
Given an additional point on , the predicted value of the function is given by:
while its variance is given by
An estimate of , say , reads:
These results fully characterize the estimates in the case of Gaussian noise. In such case, one can derive confidence bands on the estimated functions and thus quantifying the uncertainty of the estimations of any predicted value. If the Gaussian assumption does not hold, this confidence interval should be used with caution and the confidence level is only approximated. More generally, the quantification of uncertainty has been widely studied in the context of generalized additive models (see, e.g., Wood (2006) and references therein). A Bayesian approach to uncertainty quantification for this class of models is also possible; see Marra and Wood (2012).
4 Numerical simulations
In order to assess the IGS methodology, we consider two simulations on surfaces represented by NURBS for which the function to be estimated is given a priori. We aim at showing different properties of IGS in different settings. In the first simulation, the configuration is such that any method used for planar smoothing such as Thin Plate Splines (TPS) (Duchon, 1977; Wahba, 1990) can be efficiently used and thus compared to IGS. In the second simulation, the setting is such that methods for two dimensional smoothing are not appropriate, while IGS can be straightforwardly used.
4.1 Simulation 1
We consider a quarter of cylinder defined as . Using cylindrical coordinates, it is parametrized by the following mapping:
This is an isometric mapping (Stoker, 1989). That means that the mapping preserves lengths and angles, and thus the area. In other words, the parametric domain is not distorted by the mapping. This kind of mapping allows to indifferently work on the parametric domain or directly on the surface . Thus, the smoothing can be performed on the parametric domain, which is planar, namely using any traditional method for planar smoothing. In particular, in the following, we shall compare the proposed technique to TPS.
The surface is exactly representable using NURBS basis functions of degree 2 or higher which are at least globally continuous. We are interested in recovering the function:
from noisy and discrete observations. The quarter of cylinder and the function are shown in Figure 3. The function is evaluated in points located on a 1010 grid of equally spaced points in the parametric domain and is affected by independent Gaussian observational errors of zero mean and standard deviation . We generate the data for and the simulation is repeated times. We use random scalars drawn from the standard normal distribution to generate the noise (specifically, we used the MATLAB function randn). The dimension of the IGA space varies between and and we use different NURBS basis functions, namely globally , , and continuous of degrees , and respectively, obtained from refinement. The smoothing parameter is chosen at each simulation repetition and for each basis setting considered, by minimizing GCV criterion in (3.15).
The first two estimated functions over the ones are displayed in Figure 4. We observe that these are qualitatively good estimates of . In general, the number of basis functions must be chosen carefully. Indeed, this choice should depend both on the complexity of the function to be estimated and of the number of data points available. When the number of basis functions is small, IGS is not able to capture the behaviour of the function , as it would be the case with any other smoother. On the contrary, when the number of basis functions is high, we see that there is a larger variability in the estimated functions . Indeed, when the number of basis used is too high, GCV criterion can lead to overfitting, that is the estimated function incorporates noise. However, we see in Figure 4 that the estimates are not very sensitive to this choice. Finally, we remark that the minimum number of basis functions is basically dictated by the number of functions used to represent the surface with NURBS.
We notice, following Section 3.3, that we used natural boundary conditions, for which we have not formally proved the well-posedness of the problem. We report in Table 1 the mean condition number for the matrix of (3.12), with chosen with the GCV criterion and for only, since results for and are similar. As a matter of fact, the system (3.12) results to be well-conditioned in all our numerical experiments.
In order to better assess the IGS method, we use the empirical mean function and the associated empirical variance function where denotes the -th estimated function. These are both shown in Figure 5, in the case of globally continuous basis functions of degree and with basis functions. The estimates appears to have a negligible bias and a relatively small variance.
We compare our methodology with a widely used smoothing technique, TPS, using cylindrical coordinates. This method is implemented for instance in the R package mgcv (Wood, 2015). We use different number of basis functions for a comparison, selecting the smoothing parameter at each simulation repetition and for each number of basis considered via GCV. As a criterion for the comparison, we use the mean squared error (MSE) of the estimator, computed as:
where is a lattice of evaluation points on . We compute the MSE for , and and degree and for IGS and for for TPS. The comparisons of MSE in terms of the equivalent degree of freedom is shown in Figure 6. We also compare in Figure 7 the best setting of each methodology and for each number of basis functions, that is the lowest median MSE. Finally, Figures 6 and 7 show that IGS is comparable to TPS in terms of performance, in a setting where the latter technique may be applied. Figure 6 illustrates a key feature of the IGS method and more generally of all smoothing techniques. Specifically, one can notice that increasing the number of basis functions does not improve the estimated function . Indeed, although we increase the number of basis functions to build , the number of measurement points remains the same, i.e. is still built from the same set of data, but only through a richer finite dimensional space . Therefore, the convergence of to should be simultaneously regarded through the number of data and the quality of the NURBS space . We highlight in Figure 6 that the same happens with TPS-based smoothing. A key role in the convergence of the method may be played by the NURBS space when the number of measurement points is clustered in a region of , for which mesh refinement techniques can be used.
4.2 Simulation 2
We consider the surface reported in Figure 8 which is represented in terms of NURBS basis functions of degree starting from the knot vector along both the parametric directions and control points , , , , , , , , and , and with the corresponding weights vector . As reported in Figure 8, we consider the exact function to be estimated:
which is evaluated in points located on ; more specifically, they are located on a 1010 grid of equally spaced points in the parametric domain. The observations are affected by independent Gaussian observational errors of zero mean and standard deviation . The sampling of the data , , is repeated times. The sampling locations remain the same for all the repetitions. We consider three NURBS spaces of dimensions , of globally , and continuous NURBS basis functions of degrees , and , respectively. These are obtained by refinement of the original NURBS basis.
In Figure 9 we report the estimated functions for the last two simulation repetitions, by considering NURBS spaces of basis functions of degree and globally continuous on . Figure 10 highlights the empirical mean function and the corresponding empirical variance function over the simulation repetitions, obtained for the same NURBS basis functions. The results obtained for the degrees and are very similar and are not reported here for the sake of brevity. Our estimation has a negligible bias and a small variance, as shown in Figure 10. The behaviour of the function seems to be very well captured by our estimator , as shown in Figure 9. In this setting, we increase the regularity of the basis functions without changing the number of basis functions . The estimated functions are thus globally continuous. We then compare the mean integrated squared errors (MISE) of IGS for different regularity and degrees of NURBS basis functions. The MISE of any estimated function is defined as:
As we observe in Figure 11, the quality of the estimation is not significantly affected by the regularity of the basis functions.
5 Estimation of aerodynamic force on the SOAR’s winglet
We aim at estimating the pressure coefficient field and the corresponding aerodynamic force on the inboard winglet of the SOAR shuttle shown in Figure 1. The pressure coefficient is a dimensionless field related to the pressure field. It describes the relative pressure over the surface of the winglet and is given by:
where is the pressure at point , and are the far field wind speed and pressure, with the air density. The flow regime of the free stream is subsonic, namely the Mach number is , for which we can reasonably assume that the pressure coefficient field over the winglet remains sufficiently smooth, since transonic effects are marginal. The force acting on the winglet and due to the contribution of the pressure (see, e.g., Anderson (2010)) is given by: