Fast hierarchical fusion model based on least squares B-splines approximation

With manufacturing shifting from traditional products to high value products, the complexity and accuracy of the products are increasing in order to reduce energy costs, create friendly environment and better health care. Structured surfaces, freeform surfaces, and other functional engineering surfaces are becoming the core part of high value manufacturing products. However, measurement of these surfaces is becoming very difficult due to instrumental limitations including measurement range, speed, resolution and accuracy. Multi-instruments/ sensors measurement are now being developed for freeform and structured surface assessment, which requires the fusion of the data into a unified system to achieve larger dynamic measurements with greater reliability. This paper discusses the process of combining data from several information sources (instruments/sensors) into a common representational format and the surface topography can be reconstructed using Gaussian processes and B-spline techniques. In this paper the Gaussian process model is extended in order to take into account the uncertainty propagation and a new data fusion model based on least squares B-splines that drastically reduce the computational time are presented. The results are validated by two for freeform surface measurements.


Introduction
In modern manufacturing technology, manufactured surfaces are characterized by complex features, designed to meet functional specifications. In this competitive environment, the ability to rapidly design, production and inspection of the specifications of these pieces is becoming essential. To this aim, measurement with more than one sensor can improve the accuracy and the precision of the measurement and decreases the acquisition time. The main idea of a multiple sensor measurement is to use different sensors to acquire the same object, the results are then combined to improve the results of each single measurement. A classification of multiple sensor systems according to the way the information contribute to the final objective can be found in Girao et al. [1]. A multiple sensor system can be classified into: (i) complementary: the sensors are not directly dependent but their output can be combined to yield better data, under a pre-defined criterion; (ii) competitive: the sensors provide independent measurements of the same quantity They can be identical or can use different measuring methods and they compete in the sense that when discrepancies occur, a decision on which to believe must be made, sensors are put in competition either to increase the reliability of a system or to provide it with fault tolerance; (iii) cooperative: the sensors' outputs are combined to produce information that is unavailable from each sensor if used alone; iV independent: all the remaining cases, basically when the sensors' outputs are not combined. The process to combine multiple data from different sensors is defined in Weckenmann et al. [2] as combining data from several sources (sensors) in order that the metrological evaluation can benefit from all the available sensor information data.
The literature on data fusion often suggests sequential use of the sensors. In Carbone et al. [3] and Shen et al. [4] optical and contact sensors are sequentially used for reverse engineering applications. First, a point cloud is acquired with an optical system and there is a first digitization of the object. With this information, it is possible to write a CMM part program that acquires the desired point cloud. The final CAD model is thus built by refining the first raw model (based on non contact data) with the more precise CMM data points.
An example of competitive data fusion can be found in Ramasami et al. [5] where a multiscale approach based on wavelet decomposition is presented. The authors show how an appropriate multiscale model can increase the measurement result in micro-scale manufacturing.
Another example where the two sensors are used in a competitive way can be found in Xia et al. [6]. technique in order to align point clouds coming from two different sensors. In this case, the main idea is to reconstruct the information provided by the Lo-Fi sensor with a Bayesian Gaussian Process (GP) model. Then, the alignment with the Hi-Fi data is performed thanks to a local kernel smoothing technique and with a rigid transformation matrix (translation and rotation).
Starting form the approach proposed by Xia et al. [6], Colosimo et al. [7] presented a data fusion model to merge data from non contact and contact devices with a GP model. The model was applied to reconstruct a surface with application to normal and large scale metrology. The fusion process was performed using contact data to locally correct the model of the non contact point cloud.
As an alternative approach, Wang et al. [8] summarized different weighted least squares fusion techniques for homogeneous multi-sensor data fusion, which show remarkable efficiency enhancement than GPbased methods.
Yin et al. [9] proposed a Bayesian multiscale data fusion model by designing composite covariance functions. The authors proposed to use different covariance functions to estimate different features composing the surface: smooth or rough.
In their following work Yin et al. [10] proposed and intelligent sampling plan for the data fusion model based on the uncertainty of the prediction.
Following the work of Colosimo et al. [7], the uncertainty from the Lo-Fi model is taken into account in the fusion model and an equation for the uncertainty propagation is provided. However, a fatal disadvantage of a data fusion model based on GP is the high estimation time, due to the time necessary to solve a large linear system multiple times. If the number of measured point is high, it may be not possible to estimate the model parameters', due to the limited CPU memory. To overcome these issues a data fusion model based on B-splines approximation is proposed. The proposed model allows to estimate the parameters quicker compared to the existing data fusion models based on Gaussian process, achieving similar reconstruction performances.
The paper is structured as follow: in Section 2 the data fusion models are described, in Section 3 the models are applied to real test cases and in Section 5 conclusions and future developments are given.

Fusion model
In this Section, the general scheme of an hierarchical fusion model is presented. The model described in Colosimo et al. [7] is briefly recalled and a formula to estimated the uncertainty of prediction of the fusion model is presented in Section 2.1. In Section 2.2 the proposed B-splines model and its properties are discussed.
Data coming from different measurement systems with different precision and/or accuracy, called Low-Fidelity (Lo-Fi) and High-Fidelity (Hi-Fi) data, are analyzed. Usually the number of Lo-Fi data is greater than the number of Hi-Fi data, but the Hi-Fi data are more accurate. The aim of a fusion model is to link the Lo-Fi to Hi-Fi data in order to achieve better prediction of the real, unknown, surface points. The data acquired using the Lo-Fi device can be expressed by where = s u v ( , ) T is a vector describing the spatial location, s z ( ) l is the value of the measured Lo-Fi point at the spatial location s, s f ( ) is a function of the spatial location that describe the measured surface and s ( ) l represents the randomness of the measurement with the Lo-Fi device.
Once the Lo-Fi data are described, a model to merge the two sets of data can be defined as where s z ( ) h is the value of the measured Hi-Fi point at the spatial location s, s ẑ ( ) l is the prediction with the Lo-Fi model at location same location, s f ( ) is the function that describe the discrepancy between Lo-Fi and the Hi-Fi data and s ( ) is the measurement error of the difference between the Hi-Fi device and the prediction of the Lo-Fi model.
Since the Lo-Fi and the Hi-Fi data can be acquired at different spatial locations, the prediction with the Lo-Fi model has to be used in the linkage model. A further advantage of using the prediction is that the random measurement error should be removed.

Gaussian Process model
Gaussian Process (GP) is a model commonly used in spatial statistics [11] and to describe the output of computer experiments [12]. Many authors, see for example Kennedy and O'Hagan [13], Qian et al. [14], Xia et al. [6], Colosimo et al. [7], have used this model to describe booth Lo-Fi and Hi-Fi data in a data fusion framework. The observed Lo-Fi value, z l , at a generic spatial location s can be described, using a GP model, by the relation is the transpose of f s ( ) l r , which is a vector of known basis functions which act as regressors, T and s s r ( , ) i j is the correlation function modeling the dependence between two distinct points s i and s j as a function of their distance. This type of correlation structure is useful if points closed in the space have similar values and hence a good prediction at a given location can be done by looking at its neighborhood. In this paper the power exponential distance is used When the scale parameters, 1 and 2 , have the same value the correlation is called isotropic, i.e. it does not depend on the specific spatial direction; otherwise the correlation is called anisotropic. The prediction at any new point s 0 can be computed using the best linear unbiased predictor (BLUP) estimator (see for example Santner et al. [15] or Shabenberger and Gotway [16]), defined as where = r s s r ( ( , ), The objective of a data fusion model is to combine the Lo-Fi and Hi-Fi data in order to improve predictions achievable by using the Lo-Fi or the Hi-Fi data sets alone. The core of the data fusion model is a linkage or second-stage model, which represents the connection between Lo-Fi and Hi-Fi data and can be expressed as [7,14] = The aim is to "correct" the Lo-Fi predictions s ẑ ( ) l i using a "shift" effects, represented by + s . Combining the linkage model in Equation (8) and first-stage model in Equation (3), the model obtained via data fusion allows one to predict a process realization at each new location s 0 as: is the variance-covariance matrix of the linkage model, computed as where 0 is the kriging variance of the Lo-Fi model (given by the first stage equation), but computed at the Hi-Fi locations and R il the GP variance-covariance matrix.
While the variance of prediction, or kriging variance, can be computed as: Further details about the parameters estimation and the computation of the prediction variance can be found in Appendix A.

Least squares B-splines approximation fusion model
The least squares B-splines approximation (LSBA) is an algorithm that approximate the observed data using as basis functions using uniform cubic B-splines. The LSBA algorithm is first present, then the fusion model is analyzed.
The evaluation of a point of the surface, at a spacial location s, can be computed as where + + i k j l ( )( ) are the coefficients (control points) to be estimated, t v v and B · are uniform cubic B-spline defined as with < t 0 1, where • represents the floor operator. Since the basis B-splines surfaces are described through the cross product of two B-splines basis (see Equation (13)) the control points ij are defined over a rectangular lattice domain. The control lattice of the model, k , is shown in Fig. 1, it is defined on a rectangular domain and i j , the value of the ij-th control point at i j ( , ) for = … + i m 1,0, , 1 and = … + j n 1,0, , 1. In order to estimate the model parameters', the point evaluation must be rewritten in matrix notation as = + z F (15) where z n is the vector of the observed data, p is the vector of coefficients, × F n p is the sparse matrix of B-splines coefficients where B · correspond to the B-splines basis and is the vector the measurement errors, which are assumed to be independent and identically distributed normal random variables. As it is possible to observe form Equation (14) only 16 coefficients are not zero, so the matrix F is extremely sparse. The unknown coefficients can be computed minimizing the function where z i and ẑ i are the value observed and predicted at location s i while the second term ( 0) is a smoothing penalty term, where × E n n is a symmetric sparse matrix whose entries are where is the divergence operator and is the vector scalar product. The minimum value of (16) respect to can be found by solving the sparse linear system To solve the sparse system in Equation (16) the sparse Cholesky decomposition implemented in the Eigen library [17] was used.
As suggested in Hjelle and Daehlen [18], λ is chosen between the multiples of where · F is the Frobenius matrix norm.
As for the GP model the LSBA algorithm is firstly applied to the Lo-Fi data at the first stage. The observed Lo-Fi data can be described by where z l are the observed Lo-Fi data, F l is the model matrix of the Bsplines basis, l is the vector B-splines coefficients and l is a vector of i.i.d. normally distributed random variables. The LS estimation of l is The estimation of its variance-covariance matrix is Eventually, the estimate of to avoid the computation of the solution of the linear system if the number of measured points is high. The variance of the prediction in a new location s 0 can be computed as It should be noted that, since the inverse of + F F E l T l l l is not explicitly computed a linear system has to be solved for each new location.
The observed Hi-Fi points can be described as Since the prediction with the Lo-Fi model has prediction variance that is a function of the spatial location u v ( , ) this information can be taken into account assigning a weight to each prediction. The weight correspond to inverse of the prediction variance, i.e. a smaller importance is given to predictions with higher variability. The variance of the residual process, = z z ẑ h l , is a multiple of the variance of the predicted values with the Lo-Fi model the variance of the observed Hi-Fi data become The variance of the fusion model is then bigger than the variance of ẑ l because is the sum of two independent random vectors. The parameters can be estimated with the weighted least squares (WLS) procedure, the estimator of is The prediction in a new location s 0 can be computed as where

Q f s F S Q f s f s Q V I Q f s f s Q F S Q Q f s F S Q f s
is a diagonal matrix with entry equal to the mean variance of the Lo-Fi model defined as where ŝ ( ) l 2 is the predicted variance in the spatial location = s u v ( , ), S is the surface domain and S | | is its area. In order to speed up the computation it is assumed that the matrix W is a diagonal matrix. The uncertainty of the Lo-Fi model is taken into account, but not the correlation between the predictions, the prediction variance the matrix W become The data fusion model was coded in C++, the points evaluation and the variance computation was easily parallelized using OpenMP [19]. Although the computation of the prediction variance requires the solution of one (few) linear system for the Lo-Fi (fusion) model it is still

Test case
In this Section the advantages of combining different data sets with the proposed hierarchical data fusion approach are investigated using some simulated data. To this aim, performance of three different models are compared: The performances of the model estimated using both the Lo-Fi and Hi-Fi points as they both come from the same source, called addition, it is not investigated. As observed in Colosimo et al. [7], the addition model does not lead to smaller prediction if the Lo-Fi data have a systematic bias.
The performance indicator used throughout this paper is the root mean square (prediction) error (RMSE), defined as The Hi-Fi data were simulated from the following model while the Lo-Fi data model were randomly generated using the following equation: where a bias term was added. This bias term was assumed to have the following model: 10 10 , 2 2 which describes a Lo-Fi measurement system where the bias increases as the u and v coordinates increase in module. This bias term gives evidence to the need of performing data fusion. As a matter of fact, the Hi-Fi data have to be used to correct the Lo-Fi bias where needed. A total number of 250000 Lo-Fi points were generated from a regular grid of × [500 500] points in u and v directions. Since the number of the available Hi-Fi points is usually lower than the Lo-Fi sample size, only 100 Hi-Fi points were generated on a regular grid of × [10 10] points.
For the estimation of the Lo-Fi model a lattice of × 30 30 control points and a smoothing parameter = 1 were used. A × 10 10 control lattice and a smoothing parameters equal to 0.01 and 0.1 were used to estimate the Hi-Fi and fusion model respectively.
For each of the competitive models, an error map is computed as the difference between the predicted and the nominal value at each location ( s e ( ) i ). The error maps are reported in Fig. 2. The shape of the error map of the Lo-Fi model shows that the simulated bias affected the accuracy of the predictions. The error map of the Hi-Fi model shows that the error is generally smaller (because those data are assumed to be characterized by zero bias). The error map of the fusion model shows that a hierarchical model can correct the error of the Lo-Fi model where needed, using only a small set of Hi-Fi data. It should be also noted that, since the fusion model has to estimate the difference between the prediction of the Lo-Fi model and the Hi-Fi data, a higher smoothing parameter, compared to the Hi-Fi model, can be set to avoid the ondulation in the flat portion of the surface. Considering this peaks model as a toy example, 10 different realizations of the surface were run and the prediction ability of each model was tested using the RMSE.
The average time to estimate the parameters of the LSBA-based fusion model is 0.35s, so with this model is also possible to handle point clouds with a large amount of data.
The confidence interval on the mean of the RMSE computed on these 10 replications are drawn in Fig. 3. It is clear that the fusion model outperforms all the other methods in this case study.

Case study: freeform surface reconstruction
In this Section the advantages of the data fusion model applied to the reconstruction of a real smooth freeform surface are investigated. In this case study, the performance of different models are computed while varying the number of Hi-Fi points. Furthermore, three different sampling plans to locate Hi-Fi points were analyzed: (i) a uniform sampling plan according to a classic regular grid; (ii) Halton [20] and (iii) Hammersley [21] sampling strategies which are widely used in metrology [22].
In their work, Petrò et al. [23] acquired the point clouds of the free form surface shown in Fig. 4 with both a Structured Light (SL) scanner (Lo-Fi) and a Coordinate Measuring Machine (CMM) Zeiss "Prismo 5 HTG VAST" equipped with a analog probe head with maximum probing The number of points acquired with each device was 9635. Clearly, acquiring this sample size via CMM required a long acquisition time, which is almost unfeasible in real practice. As a matter of fact, it is assumed that only a subset of these Hi-Fi (CMM) data were effectively available to perform data fusion, while the remaining set was used as test set to compute the prediction errors of all the competitive methods. In other words, = n h% 9635 h data were used to reconstruct the surface with all the proposed methods, while the remaining set of = n n 9635 test h points was used to represent the "real" surface. Predictions at all the locations of n test data were computed with all the competitive methods and the corresponding prediction errors ( s e ( ) i ) will be computed.
The whole set of 9635 Lo-Fi and Hi-Fi points representing a smooth freeform surface are shown in Fig. 5. The color map in Fig. 5b represents the difference error of the Lo-Fi points, computed as the difference between the measured Hi-Fi and Lo-Fi points. The negative errors of the SL scanner are located near the peak, and on the right corner, while the positive errors are located on the left corner.
In order to compute the performance of each model as a function of the number of Hi-Fi points, the percentage h of the Hi-Fi data used to fit the surface models was changed (from = n 25 h to a maximum value of 4% of the total amount, i.e., = n 400 h Hi-Fi data points). The sampling of 100 Hi-Fi points, for each of the analyzed sampling strategy, are shown in Fig. 6. Since the "real" points were available only on a finite set of locations (the ones measured using the CMM) after computing the theoretical locations the closest point, between the points not already included in the set, were selected as Hi-Fi point.

GP-based fusion model
The surfaces reconstructed with the GP model, using 100 H-Fi points, are shown in Fig. 7. The color represents the difference between the test point and the predicted one. The errors of the Lo-Fi model has an higher magnitude on the left of the surface. On the peak the error is positive, while on the bottom left it is negative. These discrepancies are systematic due to the SL measurement system. On the contrary the systematic error of both the Hi-Fi and the fusion models appear smaller, so the reconstruction result is better. The error map of the Hi-Fi model shows that, due to the small number of Hi-Fi points available, the model has high errors in some locations of the reconstructed surface. On the contrary, the fusion model uses the information of both the Lo-Fi model and the Hi-Fi points to better reconstruct the surface reducing the error magnitude everywhere. The major differences between the Hi-Fi and the fusion models can be observed along the border of the surface, where there aren't Hi-Fi points available, or in the portion of the surface where the distance between the Hi-Fi points is higher. The GP fusion model can therefore "correct" the prediction of the Lo-Fi model using the Hi-Fi points as "attractors". Fig. 8 shows the standard deviation of the predicted mean value. The uncertainty of the points predicted using the Lo-Fi model is lower compared to the other models because the parameters are estimated using a higher number of points. Since the number of Hi-Fi points is low, the uncertainty of the predictions of the Hi-Fi model is the highest. The standard deviation using a grid sampling shows a sinusoidal pattern. This is due to the effect of the anisotropic correlation function, since the estimated correlation along the u direction is higher, the prediction uncertainty is lower. The uncertainty of the fusion model has the same behaviour of the Hi-Fi model one, but with a lower magnitude.

LSBA-based fusion model
In order to use the LSBA algorithm, the first step is to choose the three parameters that characterize the model: The surface reconstructions using the LSBA algorithm and the reconstruction errors are shown in Fig. 10. The behavior of the error is comparable with the one of the GP models, but the error of the Hi-Fi models is higher if there are not points located on the boundary of the surface. The standard deviation of the mean prediction, show in Fig. 11, is lower compared to the GP models for all the model analyzed. Since the spatial correlation of the points is not taken into account the standard deviation of the Hi-Fi model has a check board effect, it is low were there are the measured points and higher in the other locations. Differently from the GP model, the behaviour of the standard deviation of the fusion model is similar to the one of the Lo-Fi model, i.e. the uncertainty due to the second stage model is low.

Case study: portion of impeller blade
A portion of an impeller blade was measured using both a structured light scanner and a CMM. The SL measurement was carried out using Creaform HandyScan 7000, with a + L 20 60 /1000 μm uncertainty, L expressed in mm. Before the measurement, the blade was painted to matt and a set of anchor points was arranged in the exterior region of the measuring blade as the coordinate reference as shown in Fig. 15a. The painting was performed according to the ISO 3452 [24], the sensitivity of detection was equal or less than 0.5 μm. The red-bounded area was then densely scanned with the sampling resolution of 0.2 mm (0.05 mm with the refinement algorithm). The employed CMM was a Hexagon Global Classis SR 05.07.05 with the composite uncertainty in + L 2.3 3 /1000 μm, L expressed in mm. The measuring probe was a Renishaw SP25 with a 1.5 mm diameter ruby ball. Without using a CAD model, the blade was automatically measured by manually setting the initial scanning points, directions and scanning bounds as shown in Fig. 15b. Once the center coordinates of the ruby ball were obtained after scanning, the surface points were calculated by compensating a ball radius along the surface normal, which were fitted using a set of neighboring points.
A total number of 291 293 Lo-Fi (with the SL scanner) and 1616 Hi-Fi points (with the touching probe) were measured. The two point clouds were first aligned using the rigid iterative point closest algorithm implemented in CloudCompare [25]. The measured points are shown in Fig. 16 along with the estimated systematic error of the Lo-Fi data. Since, due the sampling size, is not possible to see the measured points, a decimated version is shown in Fig. 16b.   Fig. 11. Predicted surface and mean standard deviation using the LSBA model. In this second test case the systematic error appears to have non linear shape, it is positive in the top and on the bottom left portion of the surface and negative elsewhere. It is not possible to estimate the parameters using the GP because the number of the Lo-Fi points is too high, with a 64 GB RAM is not possible to compute the correlation matrix. To check the prediction performance of the GP model, the SL point cloud was decimated using a uniform sampling in the space with a minimum sampling distance distance of 0.5 mm. The decimated data, shown in Fig. 16b, has a total number of 3828 points.
In order to check the goodness of the analyzed models (Lo-Fi, Hi-Fi and fusion) a subset of Hi-Fi points was used to estimated the models' parameters, while the remaining points were used to compute the prediction error. A total number of 25, 100, 225 and 400 points were sub-sampled from the CMM data, using the grid, Halton and Hammersley samples strategies. For each sampling strategy the nominal points were first generated, the closest available Hi-Fi points were then selected to perform the analysis. Since the duplicate points were deleted, the total number of Hi-Fi points is smaller than the designed one.    The analyzed Hi-Fi points (nominal size of 100 samples) are shown in Fig. 17. The RMSEs were then computed using the remaining Hi-Fi points inside the domain of the training points. Fig. 18 shows the results of the test case. A lattice grid of × 30 30 control points and a λ parameter of 0.1 were used for the Lo-Fi model, while a grid of × 10 10 control points and λ equal to 0.01 were set for the Hi-Fi and fusion models. The execution time of the Hi-Fi model is the lowest due to the limited number of measured points, but the best RMSE can be obtained using the fusion model. It should be noted that the estimation times of both the Lo-Fi and fusion model are fast, both the models can estimate the parameters in less than 1 s.
The results of using the decimated version of the Lo-Fi data are shown in Fig. 19. Due to the limited complexity of the analyzed mesh, it is possible to achieve the same prediction performances. As in the previous test case, the results of the GP and LSBA methods are comparable, but the time required to estimate the model parameters of the LSBA method is one order of magnitude faster.

Conclusions
In this paper a LSBA-based fusion model and a formula for the propagation of the uncertainty of prediction from the Lo-Fi to the linkage model were presented. A formula to compute the prediction variance for the GP-based fusion model proposed in Ref. [7] was also proposed.
The advantages of the proposed LSBA-based fusion model were firstly explored with artificially generated data, then the prediction performances of the proposed model were then compared to the GPbased fusion model with a real test case. It was shown the RMSE of the two models are comparable, but the time required for the parameters' estimation of LSBA-based fusion model is orders of magnitude lower compared to the GP-based one. It was further observed that the uncertainty of the prediction of the LSBA model is lower compared to the GP one.
In the test case three sampling strategies were analyzed, the results shown that the sampling had no effect on the RMSE values. Other sampling plans based both on the curvature [26,27] of the Lo-Fi model or on the uncertainty prediction uncertainty [28] has to be investigated.
A second test case on real data showed that the proposed LSBA model can deal with a larger amount of points in a lower execution time compared to the GP model. The proposed LSBA model uses a B-splines basis function that  correspond to the cross product of the two B-splines basis function in u and v directions, a data fusion model based on local B-splines function, such as the truncated hierarchical B-splines [29] or the locally refined B-splines [30] has to be developed. The reconstruction of shapes from large data set can be performed by combining segmentation (to identify surface patches) and surface fitting via data fusion, new approaches for data fusion of multi-patches surfaces have to be developed.