1 Introduction
Scattered data fitting is nowadays of fundamental importance in a variety of fields, ranging from geographic applications to medical imaging and geometric modeling, see for example [8, 13, 14]
. The topic can be addressed in different approximation spaces, either by using spline spaces or radial basis functions which are particularly attractive in high dimension
[19]. In particular, thanks to their low computational cost and also to their better control of conditioning, two stage schemes relying on local approximations have received a lot of attention, being formulated either as partition of unity interpolation or as two–stage approximation methods, see for example [3, 5, 6, 7, 15] and references therein.In this work we are interested in surface reconstruction of industrial models starting from scattered data obtained by optical scanning acquisitions. In order to reduce the noise, the available data are always preprocessed, since the interest is in highly accurate reconstructions. Anyway, since the available data can not be considered exact, interpolation is not required and a less costly approximation scheme can be used to compute any local fitting. In particular, we are interested in using extensions of tensorproduct Bsplines, see e.g.,
[16], the standard choice for industrial computeraided design applications.As well known, an approximating spline can be obtained by using different approximation approaches and also operating in different kinds of spline spaces. The approximation schemes can be divided into two main classes. The first is the class of global schemes which simultaneously use all the given information and usually require the solution of a linear system whose size is equal to the cardinality of the considered spline space. The final approximation defined with the second approach can be collocated in the field of QuasiInterpolation (QI) which can require the solution of small linear systems whose size does not depend on the cardinality of the entire data set. These local systems descend from the application of local approximation schemes, each one considering a small number of data whose associated parameter values belong to a small portion of the domain intersecting the support of a compactly supported basis function. The computational advantage of the QI approach is evident, since these local linear systems are independent of each other. The quality of the final spline approximation depends in this case not only on the considered spline space but also on the approximation power of the local approximation scheme and on the choice of the local data sets.
Since the data can be characterized by a highly varying distribution, by also including voids, a flexible and reliable approach which automatically (re)constructs the geometric model, by suitably adapting the solution to the shape and configuration of the given point clouds, strongly improves the efficiency of the overall scheme. We then consider adaptive spline approximations based on the hierarchical extension of multivariate Bsplines, usually defined in commercial software tools via the tensorproduct approach, which does not provide the possibility of a full local mesh refinement. In particular, in our construction we exploit the capabilities of truncated hierarchical Bsplines (THBsplines) [11], since they allow an easy extension to hierarchical spline spaces of a quasi–interpolation scheme formulated in a standard spline space [17, 18].
The first proposal based on adaptive THBspline fitting of scattered data for the reconstruction of industrial models was introduced in [12], where a global adaptive smoothed leastsquares scheme was developed. Successively, in order to increase its locality and reduce the computational cost, the same problem has been addressed in different papers by combining a twostage approach with hierarchical spline approximations. The first contribution where this kind of schemes was used by some of the authors was presented in [2] to deal with gridded data of Hermite type. In [3] these kinds of approximants were extended for the first time to scattered data, by using in the first stage of the scheme local polynomial least squares approximations of variable degree. A preliminary application of this scheme to industrial data reconstruction was given in [4], where its theoretical analysis was also presented.
In this paper a variant of the approach considered in [3, 4]
is presented, to further decrease the number of degrees of freedom necessary to reach a certain accuracy and also to reduce the artifacts in the resulting surface. Since the local polynomials used in
[3, 4] need to be converted into Bspline form for being usable by the quasiinterpolation approach adopted in the second stage of scheme, we now work directly in local spline spaces. In this way, the algorithm for the first stage of the scheme is simplified. In order to avoid rankdeficiency problems by simultaneously controlling the smoothness of local approximants, a smoothing term is added to the least squares objective function. To improve the stability of the proposed method for general data configuration, we also inserted an automatic control on the choice of the local data sets.The structure of the paper is as follows. Section 2 presents the model problem, while the first stage of the new scheme, devoted to the computation of local smoothed least squares Bspline approximations, is described in Section 3. Section 4 introduces the construction of the adaptive THBspline surface approximating the scattered data set. A selection of examples on geometric models representing components of aircraft turbine blades is presented in Section 5, and compared with the results obtained in [4]. The numerical experiments include a new scattered data set with voids and the adaptive reconstruction of a surface closed in one parametric direction.
2 The problem
The industrial models here considered are components of aircraft turbine blades which can be suitably represented in parametric form by using just one map, with the possibility of being periodic in one parametric direction. The problem can be mathematically described as follows. Let
be a scattered set of distinct points in the 3D space which can be reasonably associated by a onetoone map to a set of distinct parameter values belonging to a closed planar parametric domain Since the choice of a suitable parametrization method for the definition of the set , which can naturally influence the quality of the final approximation, is not our focus, in this paper we relate to classical parametrization methods based on a preliminary triangulation of the scattered data set see e. g., [9, 10]. Consequently, both and are considered input data for the approximation problem.
Focusing on twostage spline approximation schemes, we can introduce the basic lines referring for simplicity to their formulation in a standard space of tensor product splines of bidegree With this kind of methods, a quasiinterpolation operator is defined so that with
denoting a vector function, possibly periodic in one parametric direction, with components in the spline space
Using a suitable spline basis of such vector spline can be expressed as follows,(1) 
where each coefficient vector is computed in the first stage of the scheme by using a certain local subset of data and the corresponding set of parameter values , so that
When dealing with discrete data, measuring the accuracy of the spline approximation with the maximum of the errors at each parameter site can appear reasonable at the first sight. However, the quality of the approximant is also strictly related to the lack of unwanted artifacts, a feature of fundamental importance for industrial applications of any approximation scheme. In this context, it is then a common practice to require the error under a prescribed tolerance only at a certain percentage of sites in .
3 First stage: local Bspline approximations
For computing each vector coefficient necessary in (1) to define the approximation we consider a local data subset ,
associated to the set
of parameter values in a local subdomain of which has non empty intersection with the support of the basis function , namely . By denoting with the set of Bsplines in not vanishing in (which necessarily includes ), the value of is defined as the vector coefficient associated with in the local spline approximation
minimizing the objective function
(2) 
where is a smoothing coefficient and the thinplate energy,
As well known, the assumption of a positive ensures that this local approximation problem admits always a unique solution, regardless of the distribution of
Since the scheme is locally applied, an automatic (datadependent) selection of the parameter could be considered. For example, the choice may take into account the cardinality of the local sample or the area of , which influence the value of the first and of the second addend in (2), respectively. In view of this influence however, we may observe that a constant value of implies that the balancing between the fitting and the smoothing term in the objective function usually increases when or the area of increases, being this true in the second case because second derivatives are involved in the smoothing term. Both these behaviors seem reasonable and are confirmed by the quality of the results obtained in our experiments, where a constant value for is suitably chosen.
Differently from [3, 4], in order to better avoid overfitting, a lower bound for the cardinality of is now required. Obviously, the value of has to be reasonably chosen with respect to the data distribution. To fulfill this condition, is initialized as and enlarged until . The refinement strategy presented in Section 4 automatically guarantees that and prevents an excessive enlargement of the set , which would compromise the locality of the approximation. For this reason, it is not necessary to set a maximum value for controlling the enlargement of the local data set. Summarizing, the computation of each is done according to the following algorithm.
Inputs

: scattered data set;

: set of parameter values corresponding to the data in ;

: tensorproduct spline space defined on ;

: tensorproduct mesh associated with the space ;

: index of the considered basis function ;

: minimum number of data required for the local approximation;
Local smoothing spline approximant

initialize ;

initialize and ;

while

enlarge with the first ring of cells in surrounding ;

update and ;


compute the local approximation for the data and by minimizing (2);

set
Output

vector coefficient .
Note that exploiting a regularized least square approximation and, as a consequence, being able to directly employ the local spline space has significantly simplified the algorithm originally proposed for the first stage in [3, 4], where a variabledegree local polynomial approximation was considered. In particular, the new scheme does not require the selection of a suitable degree for the computation of any coefficient and eliminates the conversion of the computed approximant from the polynomial to the Bspline basis.
In the following section, after introducing the THBspline basis, the operator is easily extended to hierarchical spline spaces, by also introducing the automatic refinement algorithm here considered. Note that this extension rule ensures that the coefficient associated to a THBspline basis function remains unchanged on a refined hierarchical mesh if this function remains active on the updated hierarchical configuration.
4 Second stage: THBspline approximation
Let us consider a sequence of spaces of tensorproduct splines of bi–degree defined on the closed domain and a sequence of closed domains , with and . Any , for is the union of cells of the tensorproduct grid of level . The hierarchical mesh is the set of active cells at different levels.
We denote by the Bspline basis of . For any , , let
be its representation in terms of Bsplines of the refined space . We define the truncation of with respect to level and its (cumulative) truncation with respect to all finer levels as
and
respectively. For convenience, we also define for . The THBspline basis of the hierarchical space was introduced in [11] and can be defined as
where
is the set of active multiindices, and denotes the intersection of the support of with . The Bspline is called the mother Bspline of the truncated basis function .
By following the general approach introduced in [18], we construct the vector THBspline approximation of the scattered data set in terms of the hierarchical quasiinterpolant as follows,
(3) 
where each vector coefficient is the one of the mother function and is obtained by computing the local regularized Bspline approximation on the data set associated to as described in Section 3.
In order to define an adaptive approximation scheme, we need to specify how to iteratively construct the hierarchical mesh and the corresponding spline space .
At the first iteration we always start from a hierarchical mesh coinciding with a uniform mesh , and therefore the spline space is a tensorproduct space . At each iteration the coefficients of the hierarchical QI defined by (3) are computed, and the approximation is then evaluated at each data point to get the errors
The hierarchical mesh is refined in the areas of the domain where these errors exceed the given tolerance and a refinement criterion is satisfied. The refinement criterion has been motivated by the observation that, when the parameter values corresponding to the local data set are concentrated in a small part of the support of , the quality of the approximation may be affected. More precisely, we consider a splitting of the two sides of in and uniform segments, respectively, and subdivide the support of in the resulting subregions. We then dyadically refine (in the two parametric directions) the support of the Bsplines , , , which contains at least one point where and at least data points in any of the subregions, being a given nonnegative integer. To simplify the usage of the algorithm by default we set but different values can be chosen if suitable, see e.g. the data set with voids considered in Example 3 of Section 5.
Concerning the parameter we may observe that the requirement guarantees that the points needed to compute the coefficients associated with the new functions in the first stage of the next iteration can be found not too far from the support of the functions themselves. Indeed in the algorithm introduced in the previous section, after a few enlargements, will surely include the support of a refined function of the previous level intersecting As a consequence, analogously to , a high value of contributes to the reduction of oscillations deriving from overfitting, but this value should also be low enough to guarantee that the refinement strategy can generate a hierarchical spline space with enough degrees of freedom for satisfying the given tolerance .
Once the hierarchical mesh is refined, the hierarchical space is updated accordingly. This loop is repeated until the tolerance is satisfied for a given percentage of the data points. If this condition is not satisfied within a maximum number of levels, the procedure stops anyway.
The proposed adaptive approximation method also extends to the case, not addressed in the previous works, of surfaces closed in one (or even two) parametric directions. Note that the local nature of the considered approximation approach makes the implementation especially easy, since coefficients associated with a THB present at successive steps of the adaptive refinement procedure (even if possibly further truncated) do not depend on such steps.
5 Examples
We present a selection of tests for the approximation of industrial data obtained by an optical scanning process of different aircraft engine parts. The parameter values are computed in all examples with standard parametrization methods based on a triangulation of the scattered data sets, see e. g., [9, 10]. The results highlight the effects of considering a minimum number of local data points (also) in the first stage of method, as well as the improvements obtained by introducing a regularized Bspline approximation for each local fitting with respect to the scattered data fitting scheme considered in [4]. By combining these two changes, the twostage approximation algorithm is more stable and unwanted oscillations are further reduced.
Example 1 (Tensile). In this example, we consider THBspline approximations to reconstruct a part of a tensile from the set of 9281 scattered data shown in Figure 1 (top). We compare the new local scheme based on local Bspline approximations with the algorithm based on local polynomial approximations of variable degree presented in [4], where this test was originally considered.
As the first test, we ran both algorithms with the same setting considered in [4], namely, by starting with a tensorproduct mesh with , tolerance m, and . The algorithm with local polynomial approximations with the parameter choice considered in [4] () led to an approximation with degrees of freedom, of points below the tolerance and a maximum error of m. The new scheme based on local Bspline approximations with and generated a THBspline surface with degrees of freedom that satisfies the required tolerance in of points with a maximum error of m.
As the second test, we ran both algorithms by starting with a tensorproduct mesh with , tolerance m, and . The algorithm with local polynomial approximations led to an approximation with degrees of freedom, of points below the tolerance and a maximum error of m. The surface and the corresponding hierarchical mesh are shown in Figure 1 (center). This approximation clearly shows strong oscillations on the boundary of the reconstructed surface, due to a lack of available data points for the local fitting in correspondence of high refinement levels. The scheme based on local Bspline approximations, with and , produced a THBspline surface with degrees of freedom that satisfies the required tolerance in all data points ( of points below the tolerance) with a maximum error of m. The surface, free of unwanted oscillations along the boundary, and the corresponding hierarchical mesh are shown in Figure 1 (bottom).
Example 2 (Blade). In this example, we test the second example considered in [4] on the set of 27191 scattered data representing a scanned part of a blade shown in Figure 2 (top). Again, to compare the new local scheme with the algorithm based on local polynomial approximations there considered, we ran both algorithms with the same setting of [4], namely, by starting with a tensorproduct mesh with , tolerance m, and . The algorithm with local polynomial approximations with the parameter choice considered in [4] () led to an approximation with degrees of freedom, of points below the tolerance and a maximum error of m. The new scheme based on local Bspline approximations with and generated a THBspline surface with degrees of freedom that satisfies the required tolerance in of points with a maximum error of m. The surface and the corresponding hierarchical mesh are shown in Figure 2 (bottom).
Example 3 (Endwall). In this example, we illustrate the behavior of the adaptive algorithm on data sets with voids by considering the reconstruction of an endwall part from the scattered data set of 43869 points shown in Figure 3 (top). The figure shows that in this case the data set represents a model with three different holes, where no input data are available. The aim of this reconstruction is to avoid artifacts due to lack of points and obtain a sufficiently regular surface (eg. by avoiding selfintersections), that can be postprocessed with standard geometric software tools to obtain a suitably trimmed model. Consequently, not only the number of points in the local data sets is important to reach this aim, but also their distribution. To properly address this issue, we consider a real density parameter with value between 0 and 1 which determines whether the distribution of the points in the local set is reliable or not for the fitting. The distribution of the local data points is computed as the number of mesh cells of level inside the support of or its enlargement, which contain at least one point, over the total number of mesh cells, either in the support of or its enlargement. If this ratio is below a density threshold, then more data points are required and the function support is enlarged for the computation of the local approximation in the firststage of the method. The approximation is developed by starting from a tensorproduct mesh, with , , and . A choice of the density parameter equals to permits to take care of the distribution of data points in the costruction of the approximation. By considering a tolerance m, the refinement generated a THBspline approximation with degrees of freedom, of points below the threshold and a maximum error of m. The surface and the corresponding hierarchical mesh are shown in Figure 3 (bottom).
Example 4 (Airfoil). This example illustrates the behavior of the new adaptive fitting algorithm with local Bspline approximations for surfaces closed in one parametric direction. We test the scheme to reconstruct a blade airfoil from the set of 19669 scattered data shown in Figure 4 (top). We ran the algorithm by starting with a tensorproduct mesh with , tolerance m, and . The local algorithm was used in this case with , , and produced an approximation with 3446 degrees of freedom, that satisfies the required tolerance in of the data points with maximum error m. The surface and the corresponding hierarchical mesh are shown in Figure 4 (bottom). By trying to force additional refinement, some oscillations appear. In this case, they are consistent with the data distribution since there are clusters of high density points, due to scan noise. Consequently, they do not represent artifacts caused by regions with very low density of data and cannot be prevented by exploiting the bound for cardinality of the local data sample governed by .
Acknowledgements
Cesare Bracco, Carlotta Giannelli and Alessandra Sestini are members of the INdAM Research group GNCS. The INdAM support through GNCS and Finanziamenti Premiali SUNRISE is gratefully acknowledged.
References
 [1] Saboret, L., Alliez, P., Lévy, B., RouxelLabbé, M., Fabri, A., Triangulated Surface Mesh Parameterization. In CGAL User and Reference Manual. CGAL Editorial Board, 5.0 edition, 2019.
 [2] Bracco, C., Giannelli, C., Mazzia, F., Sestini, A. Bivariate hierarchical Hermite spline quasiinterpolation. BIT Numer. Math. 56, 1165–1188 (2016)
 [3] Bracco, C., Giannelli, C., Sestini, A. Adaptive scattered data fitting by extension of local approximations to hierarchical splines. Comput. Aided Geom. Des. 52–53, 90–105 (2017)
 [4] Bracco, C., Giannelli, C., Großmann, D., Sestini, A. Adaptive fitting with THBsplines: Error analysis and industrial applications. Comput. Aided Geom. Des. 62, 239–252 (2018)
 [5] Cavoretto, R., De Rossi, A., Perracchione, E., Optimal selection of local approximants in RBFPU interpolation, J. Sci. Comput. 74, 1–22 (2018)
 [6] Davydov, O., Morandi, R., Sestini, A., Local Hybrid Approximations for Scattered Data Fitting with Bivariate Splines, CAGD 23, 703–721 (2006)
 [7] Davydov, O., Schumaker, L., Interpolation and scattered data fitting on manifolds using projected Powell–Sabin splines, IMA J. of Numer. Anal. 28, 785–805 (2008)
 [8] De Marchi, S., Iske, A., Santin, G., Image Reconstruction from scattered Radon data by weighted positive definite kernel functions, Calcolo 55:2 (2018)
 [9] Floater, M.S.: Parametrization and smooth approximation of surface triangulations. Comput. Aided Geom. Des. 14, 231–250 (1997)
 [10] Floater, M.S., Hormann, K.: Surface parameterization: a tutorial and survey. In: Dodgson N. A., Floater M. S., Sabin M. A. (eds.) Advances in Multiresolution for Geometric Modelling, Mathematics and Visualization, pp. 157–186. Springer, Berlin, Heidelberg (2005)
 [11] Giannelli, C., Jüttler, B., Speleers, H. THBsplines: The truncated basis for hierarchical splines. Comput. Aided Geom. Des. 29, 485–498 (2012)
 [12] Kiss, G., Giannelli, C., Zore, U., Jüttler, B., Großmann, D., Barner, J. Adaptive CAD model (re)construction with THBsplines. Graph. Models 76, 273–288 (2014).
 [13] Le ThiThu, N., NguyenTan, K., NguyenThanh, T., Approximation of triangular Bspline surfaces by local geometric fitting algorithm, NICS 2016, Proceedings of 2016 3rd National Foundation for Science and Technology Development Conference on Information and Computer Science 7725674, 91–96 (2016).
 [14] Majdisova, Z., Skala, V., Smolik, M., Incremental Meshfree Approximation of real geographic data, Applied Physics, System Science and Computers III, Lecture Notes in Electrical Engineering 574, 222–228 (2019).
 [15] Perracchione, E., Rational RBF–based partition of unity method for efficiently and accurately approximating 3D objects, Comp. Appl. Math. 37, 4633–4648 (2018)
 [16] Prautzsch, H., Boehm, W., Paluszny, M., Bézier and Bspline techniques, Springer, Berlin, 2002.

[17]
Speleers, H. Hierarchical spline spaces: quasiinterpolants and local approximation estimates. Adv. Comput. Math.
43, 235–255 (2017)  [18] Speleers, H., Manni, C. Effortless quasiinterpolation in hierarchical spaces. Numer. Math. 132, 155–184 (2016)
 [19] Wendland, H., Scattered data approximation, Cambridge Monographs on Applied and Computational Mathematics 17, Cambridge University Press (2005).
Comments
There are no comments yet.