Spatio-Temporal Adaptive Sampling for effective coverage measurement planning during quality inspection of free form surfaces using robotic 3D optical scanner

In-line dimensional inspection of free form surfaces using robotic 3D-optical scanners provide an opportunity to reduce the mean-time-to-detection of product quality defects and has thus emerged as a critical enabler in Industry 4.0 to achieve near-zero defects. However, the time needed to inspect large industrial size sheet metal parts by 3D-optical scanners frequently exceeds the production cycle time (CT), consequently, limiting the application of in-line measurement systems for high production volume manufacturing processes such as those used in the automotive industry. This paper addresses the aforementioned challenge by developing the SpatioTemporal Adaptive Sampling (STAS) methodology which has the capability for (i) estimation of whole part deviations based on partial measurement of a free form surface; and, (ii) adaptive selection of the next region to be measured in order to satisfy pre-defined measurement criterion. This is achieved by first, modelling spatiotemporal correlations in the high dimensional Cloud-of-Points measurement data by using a dimension reduced space-time Kalman filter; then, dynamically updating the model parameters during the inspection process by incorporating partial measurement data to predict entire part deviations and adaptively choose the next critical region of the part to be measured. The developed STAS methodology enhances the current free form surface inspection models, which are mostly based on spatial analysis; into spatio-temporal model, which uses (i) the spatial analysis to model part deformation; and, (ii) temporal analysis to model autoregressive behaviour of the manufacturing process for prediction of next part deviations. This provides capability to predict the whole part deviation based on partial measurement information and consequently reduces measurement cycle time. The industrial case study using a robotic 3D-optical scanner for the measurement of an automotive door inner part demonstrates the STAS methodology, which resulted in (i) a 3 Sigma error of prediction of whole part deviations within 0.27mm based on measurement of 33% of the part surface; and, (ii) a corresponding CT reduction of 42.2% from 510.5 s required by current best practice to measure the whole part to 295.18 s required to partially measure the part.


Introduction
Real-time product geometry assurance plays a vital role in achieving near-zero defects in a smart manufacturing environment [1][2][3][4][5] and is especially important for products with free-form surfaces which have high aesthetical and functional requirements as defined by tight geometric dimensioning and tolerancing (GD&T). Dimensional and geometric deviations of free-form surfaces can lead to numerous quality related problems including (i) high rate of re-work or scrap, (ii) inferior product functional performance, (iii) tooling failures; and, (iv) unexpected production downtime [6] thereby, reducing both product quality and production throughput. Furthermore, around 60-70% of design faults found during production ramp-up in an automotive body assembly process are also linked to dimensional variation [1,2]. These factors together with increasing use of free-form surfaces in automotive and aerospace industries have made inspection of geometric and dimensional characteristics of product essential [7,8].
In order to inspect products with free form surface effectively, measurement gauges need to satisfy a number of requirements as related to their: (i) accuracy, repeatability and reproducibility (R&R); (ii) coverage, i.e., capability to measure sufficient number of points on a single part to estimate geometric features' errors; and, (iii) mean timeto-detection (MTTD) of all required quality defects. of free-form surfaces not only by using contact-type gages such as coordinate measuring machines (CMMs) but also by using non-contact gages such as 3D-optical scanners (also referred to as 3D-scanners in this paper). The contact-type gages have high accuracy, however, they are inherently slow in capturing data points and the inspection process is very time consuming, resulting in relatively small number of points measured on a single part [9,10], i.e., small measurement coverage. As a result, it is difficult and time-intensive to estimate geometric deviations and variation patterns, therefore, contact gages have relatively large MTTD of quality defects.
On the other hand, the accuracy of 3D-optical scanners have improved greatly and can provide a viable solution for inspection of freeform surfaces [11]. 3D-optical scanners have the capability to capture a large number of measurement points on a single part (i.e., point clouds, also called cloud-of-points (CoP), with millions of captured data points) which allows for estimation of geometric deviations and variation patterns, and thus, 3D-optical scanners provide good measurement coverage.
In addition, efforts are also underway to reduce MTTD by changing the configuration and the placement of the robotic-inspection station, from (i) off-line measurements conducted in a stand-alone CMM-roomthis configuration has very high MTTD as the parts need to be taken out of the assembly line and transferred to the CMM-room ( Fig. 1(a)); through to (ii) by-pass (or near-line) measurement conducted next to the production line ( Fig. 1(b)); until the most recent development of (iii) in-line measurements conducted within the flow-through of the production line, with potential MTTD being at production CT level ( Fig. 1(c)).
In-line inspection station provides several key advantages, some of them are summarised below: (i) Significant reduction of MTTD of quality defects by providing realtime process information without removing parts from the manufacturing line. This eliminates the lead-time of off-line measurement system, which includes the time necessary for (a) transfer of part or subassembly from production line to CMM-room for conducting measurements; (b) planning and scheduling measurements; and, (c) transfer of measured part/subassembly from CMMroom back to production line. (ii) Elimination of out-of-sequence production error, which is important in case of multiple product variants being assembled on a single line. The unplanned sequence change leads to potential outof-sequence production error. (iii) Improved quality control due to 100% sampling capability which enables real-time detection of product quality related issues [2,12]. For instance, in-line inspection systems have higher fault detectability, which is at 56% compared to 19% for off-line inspection systems during production ramp-up of an automotive body assembly process when applied methods for root cause analysis of 6-sigma failures (variation reduction approaches) [1,[13][14][15]. Additionally, in case of some manufacturing applications, in-line inspection is the only solution capable of enforcing strict tolerances [16].
In spite of the aforementioned advantages, at present, measurement systems implemented are mostly individual point-based measurement systems instead of 3D scanners needed for free-form surfaces [17][18][19][20][21][22][23]. Currently, 3D scanners are mainly used for reverse engineering and offline inspection [24,25]. With a few studies investigating the application of 3D scanners for off-line process quality monitoring [26,27]. The lack of application of in-line measurement systems for estimation of geometric error patterns is especially critical for high volume manufacturing processes such as the ones used in automotive industry, for example, automotive body assembly.
One key challenge limiting effective application of 3D scanners in high volume production is that of inspection time exceeding production CT. In general, the inspection time conducted using robotic 3D scanner includes: (i) robot travel time; and, (ii) processing time of the scanner head. For example, during the inspection of automotive door inner part using robotic 3D-optical scanner, the ratio between robot travel time and processing time is approximately 35% :65 % of total CT, respectively. In case of high volume automotive body assembly processes, the CT of production is short, which is typically about 50-120 s, as compared to the time required to measure the complete automotive body sheet metal part or subassembly, which can typically take up to 5-15 min. Currently, the in-line 3D scanners are used for (i) estimation of geometric error patterns by measuring 100% of surface, which is feasible for in-line inspection of small stamped parts (less than 25 × 25 cm) moving on a conveyor [28]; and, (ii) estimation of key product characteristics such as gap and flushness between closure panels and automotive body openings at selected locations [35].
This paper addresses the challenge of reducing inspection time for estimation of geometric error patterns in large 3D free-form surfaces by developing the Spatio-Temporal Adaptive Sampling (STAS) methodology. The STAS methodology has the capability to predict deviation error of the whole part based on partial measurement data of a free-form surface. This is achieved by first, modelling the spatio-temporal correlations in the CoP data from a batch of training parts using space-time Kalman filter; then, dynamically updating the model parameters during the inspection process to incorporate information from partial measurement data to predict entire part deviations. Following which the sequence of regions of the part to be measured is chosen adaptively in real-time according to a predefined coverage criterion that can take into consideration (i) part inspection: as defined by the effectiveness of detecting out-of-specification regions on the part, i.e., regions not fulfilling design requirements; or, (ii) quality control: as defined by the effectiveness of detecting out-of-control regions on the part, i.e., regions not fulfilling statistical process control rules for out-of-control regions.
As an outcome, the proposed methodology minimizes necessary part inspection coverage, thereby reducing inspection time needed for a given part produced on a specific manufacturing line. The contributions of this paper are twofold: (i) development of a methodology for in-line prediction of entire part deviation pattern with partial measurement coverage by utilising spatio-temporal correlations in parts' deviations; and, (ii) development of an adaptive region selection strategy to identify the critical regions to be measured. Effective in-line application of 3D free form surface is made possible by the aforementioned contributions which reduce the inspection time of the entire part.
The remainder of the paper is organized as follows: Section 2 provides a detailed review of quality inspection of free-form surfaces; Sections 3 and 4 describe the problem formulation and the STAS methodology, respectively; Section 5 demonstrates, verifies and validates the STAS methodology on an automotive door inner part, it also compares the STAS methodology with state-of-art spatial prediction techniques, and state-of-practice static measurement plans. Finally, conclusions and future research directions are discussed in Section 6.

Literature review
Quality inspection of free-form surfaces as described in Section 1 and illustrated by Fig. 1 can be classified into (i) off-line, (ii) by-pass and near-line, and, (iii) in-line inspection. However, the existing literature can be categorized into two main groupings by considering bypass and near-line inspection as a special category of off-line inspection with shorter MTTD of quality defects. A detailed review of the two afore-mentioned types of inspection is provided below.

Off-line inspection
Existing literature for off-line measurement systems mainly deal with point-based measurements, typically CMMs fitted with tactile probes or non-contact optical scanners. As described in Section 1, point measurements are time consuming and a large number of measurements are required to estimate the shape error of free-form surfaces. Therefore, many sampling strategies have been proposed to plan measurements using CMMs effectively.
Kriging based adaptive sampling plans for point-based measurements using a CMM was proposed in [17,29], where sequential plans for points to be measured were identified using Kriging variance. The methodologies used the information up to the latest measurement to adaptively choose the next point to measure so as to evaluate straightness and roundness. A Gaussian process based methodology was further developed to obtain adaptive inspection plan for CMM's using prediction uncertainty and Euclidean distance criteria in [30]. Similarly, AK-ILS, a Kriging based adaptive inspection methodology for large aeronautical surfaces was proposed in [18]. The study optimised points to be measured on the surface according to a probability of conformance to tolerance criteria.
An adaptive point measurement planning methodology that uses a small number of sampled points to find the minimum deviation zone using neighbourhood search approach was developed in [20]. The methodology however is applied only to planar surface to evaluate flatness tolerance. A B-spline based regression methodology was used to optimise the sampling of points to ascertain tolerance requirements of surfaces with sharp-edged features in [21] and it was shown to have similar accuracy as Kriging based approaches at a lower computational cost.
Presently 3D-scanners are predominantly used for reverse engineering and off-line inspection [24,25]. The main aim of reverse engineering is to generate the CAD file of the part from the measured CoPs; a detailed review of existing techniques can be found in [25].
While off-line inspection by 3D-scanners can be performed using commercially available systems, research focuses mainly on (i) improving the quality of obtained measurements, and (ii) optimising the path for scanning the part [27].
The quality of CoP output from the 3D-surface scanner depends on shadows, occlusion and surface reflection, factors that have been considered in [26] to develop a dynamic off-line measurement methodology with a focus on obtaining high-quality CoP. Robot path optimisation can be performed to improve CoP quality by considering scanner parameters such as viewpoint and field of view, as studied in [31,32].

In-line inspection methodologies
In-line point-based measurements with optical CMM have been used in the industry using commercially available devices for 100% dimensional inspections [1,33]. A Kalman filter based measurement uncertainty reduction for in-line point-based measurements from CMMs was studied in [19,34].
While in-line point-based measurements are widely used, there is very little literature dealing with in-line CoP-based free-form surface inspection using 3D-scanners. In-line CoP-based inspection of shape of stamped sheet-metal parts utilising photogrammetry has been proposed in [28], the approach, however, is applicable for small parts up to 25 cm by 25 cm. Additionally, in-line 3D surface scanners are used in industry for partially measuring key product characteristics such as gap and flushness between closure panels and body at selected locations [35].
While research is lacking on in-line adaptive measurement technique for CoP-based measurements of free-form surfaces in sheet metal assembly systems, it has nonetheless been applied in the semiconductor fabrication industry [36,37] where it focussed on optimum selection of the batches to be sampled or the frequency of sampling in each batch to minimise semiconductor wafers at risk [38]. Whereas, the presented research focuses on optimal selection of regions of a given part to be measured while measuring all parts passing through the assembly system. Additionally, the size of semiconductor wafers typically ranges between 200-450 mm [39] and use measurement sensors that are accurate in the sub-micrometer scale [40], whereas the present study focuses on large free-form surfaces between 1-3 m.

Summary of the state-of-art methodologies for quality inspection of free-form surfaces
A table summarising the afore-discussed literature review for quality inspection of free-form surfaces is presented in Table 1. It shows that existing methodologies deal with (i) off-line and in-line point-based inspection, (ii) off-line CoP-based inspection, and (iii) in-line CoP-based inspection for small parts.
To enable effective in-line inspection of large free-form surfaces, this paper develops the STAS methodology, which, in contrast to existing methodologies that utilise only spatial correlations, utilises the spatio-temporal correlations in part deviations to predict entire part deviations from partial measurements. To date, methodologies modelling spatio-temporal correlations have been utilised for effective weather prediction, disease mapping, and modelling energy production [41,42]. The necessity and advantages of modelling spatio-temporal correlations in part deviations is explained in detail in the following section.

Problem formulation
The state-of-art modelling and analysis methods for inspection of free-form surface deviations are based only on spatial correlations assuming that surface deviations of any two consecutive parts are independent of each other during elapsed time of production. However, in many manufacturing processes, observations are not independent, but are serially (temporally) correlated, such as in multi-station assembly processes [43] and in machining and forging operations [44]. In case of a multi-station assembly process, temporal correlation is induced by process dynamics as all the parts and intermediate subassemblies go through the same Place-Clamp-Fasten-Release (PCFR) operations, i.e., are: (i) located and clamped by the same fixtures, (ii) fastened by the same equipment (for example, robotic welding cells); and, (iii) undergo the same sequence of processes and operations in the assembly line. The induced autocorrelation has a significant impact on various manufacturing performance measures as reflected by prediction errors incurred in case of ignoring dependence [45,46].
Considering both spatial and temporal correlations while modelling free-form surface deviations provides the opportunity to use incomplete measurement data of a given part to predict the whole part deviation errors by augmenting missing data using spatio-temporal correlations obtained from the training data. This spatio-temporal correlation is schematically illustrated in Fig. 2. Where, y(i, T) and y(i, T + 1) represent perpendicular to surface deviations from design nominal of points i = 1, …, 4 of the part manufactured at time T and T + 1, respectively; C s ab and C ab represent the spatial and temporal correlations between points a, b, respectively. This paper focusses on modelling and predicting free-form surface deviation pattern from its design nominal as it passes through a given station in an assembly line. Surface normal deviation at a point on the surface of the part is represented by y(s, T), where, s = 1, 2, …, N is the spatial index of a point in the mesh representation of the part, N is the total number of mesh nodes, and T is a part number in the sequence of entering the assembly line as used in this paper.
The surface normal deviations of the entire T th part is represented by Y(s, T) = {y(s, T)|s = 1, 2, …, N}. In contrast to continuous surface representation with infinite number of points, discrete mesh representation is conducive for assembly process simulation through Finite Element Method (FEM) and can easily be integrated into existing simulation tools.
The measured part deviations Z(s, T) is the sum of true deviations Y (s, T) and error Γ(s, T) as in Eq. (1). The error term Γ(s, T) represents the sum of measurement and input part uncertainty and is assumed to be Gaussian. Measurement uncertainty is a characteristic of the 3D scanner and is specified by the manufacturer. Whereas, input part uncertainty refers to the variation in shape of the part due to upstream assembly or manufacturing process variation. In this study, input part is assumed to be ideal and shape error is due to variation in the present assembly stage, thus input part uncertainty is zero.
The temporal correlation of part deviations Y(s, T) is modelled as the sum of a temporally evolving component W(s)Y(s, T − 1) and a Gaussian noise process Θ(T) as in Eq. (2).
where W(s) is the N × N transition matrix quantifying temporal dependence. In the presented study first order temporal correlations are considered, however the formulation can be extended to include higher order temporal correlations [47]. Eqs. (1) and (2) together form the state space representation of the part deviations and are called measurement and state transition equations, respectively. Modelling the part deviations using Eqs. (1) and (2) is not efficient as it requires inverse of N × N matrices which needs O(N 3 ) computations and becomes infeasible for real-time analysis when N is large. This limitation can be overcome by utilising the fact that deviations Y(s, t) are spatially correlated due to geometric covariance [48] and thus, can be effectively represented as a linear combination of K orthogonal bases and an error term as described in Eq. (3). Typically, K < < N thus achieving dimension reduction and enabling to model large number of points with K columns (basis vectors) of Φ(s), which enables part deviations prediction with the inverse of K × K matrices.
where Φ(s) represents an orthogonal basis set, and Ψ(T) represents the Table 1 Summary of literature review for quality inspection of free-form surfaces. The proposed STAS methodology* *The proposed STAS methodology is most beneficial for in-line measurement of large part  (1) and (2) we obtain the dimension reduced measurement and state transition equations as in Eqs. (4) and (5), respectively.

Point-based inspection CoP-based inspection
where Z(s, T) is an N × 1 vector of deviations of all mesh nodes of the part at time T, Φ(s) is the N × K system matrix with K orthogonal basis vectors as its columns, Ψ(T) is the vector of state variables of length K which evolve in time, H is the K × K modified state transition matrix, E (s, T) is the error vector of length N that combines uncertainty due to measurement (Γ(s, T)) and dimension reduction s T (˜( , )), and T ( ) is the state variable error vector of length K. Φ(s) models the spatial correlation, whereas H models the temporal correlation. The covariance matrices of the error terms E(s, T) and T ( ) are represented by R and Q, respectively, they are independent of each other and constant. The dimension reduced state space formulation of part deviations in Eqs. (4) and (5) is a modified version of the formulation by [49] without the spatially correlated error term, as the part deviations patterns were explained by the combination of basis vectors with no spatial structure in the error E(s, T).
The state variables Ψ(T), are low dimensional latent variables which indirectly influence the entire part deviations. The key idea in the proposed methodology is to predict the optimum values of Ψ(T) with partial measurements s T Z ( , ). The information obtained from partial part measurement is utilised to update Ψ(T) from which the entire part deviations is predicted and the new region to be measured is adaptively chosen based on a predefined coverage criterion described in Section 4.3.
The minimum mean squared error estimate of Ψ(T) is obtained using Kalman filter prediction and update recursions. The Kalman filter is a Bayesian on-line regression model where the latent and observed variables have a Gaussian distribution [47,50], and has been used extensively for forecasting, control applications and fault detection [51,52].
The application of Kalman filter in this paper is motivated by the need to dynamically update predictions in-line as new measurements arrive, and the ability of the Kalman filter to deal with missing values [47]. Which enables to predict part deviations with partial measurements. The computational effort to model spatio-temporal correlations scales linearly with time in the Kalman filter as opposed to cubic complexity when using Kriging or Gaussian process regression approaches [53]. Finally, the dimension reduction using orthogonal bases reduces redundancies in measurement data and enables effective part deviations pattern prediction with partial measurements. The scope of the proposed STAS methodology is limited to shape or form variations of the free form surface and does not include high frequency deviations such as roughness and waviness.

Methodology
The proposed methodology models the spatio-temporal correlations imparted to part deviations by the assembly system and utilises these spatio-temporal correlations to make prediction of complete part deviations pattern based on partial measurements. The state variables in the prediction model are updated after each measurement to incorporate new information and provide updated predictions of deviations at unmeasured regions. The new region to be measured is chosen based on the coverage criterion detailed in Section 4.3. The adaptive selection of regions in a part continues until the stopping criteria detailed in Section 4.4 is satisfied. The methodology is schematically illustrated in Fig. 4 and explained in detail below.

Initial model development
Modelling the spatio-temporal correlations for a given part deviation in an assembly system is achieved through the estimation of the system matrix Φ, the state transition matrix H, and the covariance matrix of error terms E(s, T) and Θ(s, T), i.e., R and Q, respectively.
The modified state transition matrix H, and covariance matrices R and Q are calculated by method of moment estimators [49] utilising training part deviations data. The required training data can be obtained from two sources, namely, (i) complete part measurements from the assembly line, or (ii) part deviations data obtained from simulation of part fabrication or assembly process. The system matrix Φ as described in Section 3 is a matrix whose columns are orthogonal bases obtained from decomposition of the centred training part deviations data, i.e., the global mean is subtracted from the deviations data.
The discrete surface representation of large parts have a high number of mesh nodes (typically greater than 30000), notwithstanding the application of dimension reduced Kalman filter, the methodology requires handling of large error covariance matrices. Which to be handled in real-time need a large system memory (RAM). To overcome this limitation a subset of all nodes called key points (ζ k ) are chosen and used as surface measurement points, and are utilised to update the state variables Ψ during partial measurement. The key points are chosen uniformly from all the mesh nodes, and the distance between two key points is chosen based on the Mean Absolute Error (MAE) of prediction and computational time of the resulting model.
Additionally, a large part cannot be measured completely by a single snapshot using the 3D-optical scanner. The robotic measurement plan manually programmed to measure the entire part requires a large number of snapshots. For instance, the automotive door inner part used in the case study requires 83 snapshots (Fig. A.15). The size and number of the regions depends on; (i) field of view of the scanner; (ii) orientation of scanner with respect to the surface of the part; and, (iii) surface reflectivity of the part. Ideally, a region is the surface area of the part which can be measured by the 3D-optical scanner in one snapshot. However, to demonstrate the ability to predict entire part deviations with partial measurements, in this study, the part is divided into larger non-overlapping rectangular regions. The division of the part into larger regions also enables the demonstration of effectiveness of the developed adaptive region selection methodology against a deterministic measurement plan, as described in Section 5.2.3 by eliminating the high computational requirements of evaluating a large number of deterministic sequences (for instance, the part with 83 regions requires evaluation of 83 C 1 + 83 C 2 + ⋯ + 83 C 83 = 9.67 × 10 24 sequences).
The assumption of rectangular regions necessitates more than one snapshot to cover a given region due to the smaller field-of-view of the sensor, and, visibility and surface reflectivity constraints. In this study, the CT for partial measurements was calculated using the manually programmed measurement plan with 83 snapshots, however considering only snapshots that covered at least one point in the regions measured. Following the identification of snapshots corresponding to the regions to be measured, the CT for adaptive partial measurement of the test parts was estimated using ABB RobotStudio® software for all adaptive partial measurement sequences.
In the following discussion, synonymous with time, a part with index T − 1 precedes the part with index T. Each part T is measured iteratively with a region being measured during each partial measurement. The notation s T t Z ( , | ) represents all points measured up to t th partial measurement of the part T. Whereas, s T t Ŷ ( , | ) represents the prediction of Y for the whole part T given measured points s T t Z ( , | ) (equivalent to filtering step in traditional Kalman filter terminology, but with observations missing for unmeasured regions). Similarly, + s T T Ŷ ( , 1| ) represents the prediction of Y for the whole part T + 1, given all the partial measurements of part T (equivalent to prediction step in traditional Kalman filter terminology). The above representations are also applicable to the state variables (Ψ) and the covariance matrix of state variables (P).

State variables' update and part deviations prediction
After the initial model development, the update of state variables to incorporate partial measurement data, and the prediction of complete part deviations is carried out in the following two stages:

Update of state variables (Ψ) and prediction of part deviations for the current part (T)
The update reflects the assimilation of partial measurement data into the state variables Ψ. Whereas, prediction refers to the prediction of deviations at unmeasured points on the surface of the part. The updated values of state variables T t ( | ) and the covariance matrix of state variables P(T|t) are estimated by first setting T T ( | 1) and P (T|T − 1) to their initial estimates Ψ 0 (obtained from orthogonal decomposition of training data) and P 0 (a diagonal matrix with large principal diagonal elements), respectively in Eqs. (6) and (7), the Kalman filter recursion equations which are derived in [47,50].
where G(T) is the Kalman gain which is estimated using Eq. (8).
During these updates, the ability of Kalman filter to deal with missing values is utilised [47], i.e., points of the part which are not measured are considered to be missing values. Accordingly, s T t Z ( , | ) in Eq. (6) contains deviation information from regions of the part measured up to iteration t, rather than all the points in the part. The size of s T t Z ( , | ) increases by the number of points in each region measured in each iteration.
represents the system matrix Φ modified to take into account the nodes with missing deviation information. Similarly, in Eq.
where s 0 = 1, 2, … # of unmeasured nodes, and Φ(s 0 ) is the orthogonal bases matrix corresponding to unmeasured points. The variance of predicted values is estimated by Eq. (10).
where P(T|t) is obtained from (7). The region to be measured in each iteration is selected based on coverage criterion described in Section 4.3. Following the measurement of a new region, T t ( | ) and P(T|t) are updated utilising Eqs. are made. This iterative measurement is represented by the region loop in Fig. 4 and continues until the stopping criterion described in Section 4.4 is satisfied, following which the prediction for part T + 1 is carried out as described below.

Update of state variables and prediction of part deviations for the next part (T + 1)
When the stopping criterion for adaptive region selection of part T is attained, index T is incremented to T + 1, and the measurements of all regions up to t th partial measurement together represent the measurements acquired for part T. The final estimates for T t ( | ) and P(T|t) for part T obtained from Eqs. (6) and (7) are substituted in Eqs. (11) and (12) to obtain the prediction of state vector and its covariance matrix for the part T + 1.
Finally, the prediction of part deviations for the unmeasured part T + 1 (called zeroth prediction) and variance matrix of predicted deviations for the part T + 1 are obtained using Eqs. (13) and (14).
The afore-discussed sequence of update and predictions is represented by the part loop in Fig. 4. The selection of first region of the part T + 1 to be measured is based on the predictions for part T + 1 obtained using Eqs. (13) and (14).
The state variables + T T ( 1| ) and the covariance matrix P (T + 1|T) for the part T + 1 obtained using Eqs. (11) and (12) are used in Eqs. (6) and (7) as T T ( | 1) and P(T|T − 1) during subsequent partial measurements of part T + 1 to update the state variables and their variance matrix as described in Section 4.2.1.

Adaptive region selection
A critical aspect of the methodology is to adaptively identify the next region of the part to be measured according to a specified coverage criterion. In this paper, the coverage criterion intends to improve the effectiveness of part inspection, i.e., improve the detectability of out-ofspecification points in a region of the part with respect to given form tolerance specification for profile of a surface. The proposed criterion chooses the region with maximum expected information gain with-respect-to the probability of unmeasured points in the region being outof-specification, and thus enables automatic selection of critical region to be measured without user intervention.
The selection of the region is adaptive because the regions to be measured are estimated utilising the complete part deviations predicted after assimilation of partial measurement data, causing the measurement sequence to be adapted for each part by considering its deviation pattern. The proposed adaptive measurement sequence as illustrated in Section 5.2.3, performs similar to the best deterministic measurement sequence which can only be identified retrospectively when all part deviation patterns are known. The mathematical formulation of the proposed criterion is explained in detail in the paragraphs below.
The probability of an unmeasured point y(s, T) being out of specification, P y(s,T) , can be calculated by taking into account that points on the surface are (i) independent from each other, or (ii) correlated with each other. In this paper, each point in the unmeasured regions is considered independent and P y(s,T) is calculated by using Eq. (15).
where F represents the CDF of a Gaussian distribution, USL and LSL are Lower and Upper Specification Limits, respectively. Pr 1 and Pr 2 are non-conformance probabilities calculated based on the predicted deviation at a given point which follows a Gaussian distribution with mean and variance given by Eqs. (9) and (10) or Eqs. (13) and (14), according to the Kalman filter model assumptions.
The assumption of independence causes the part regions with points such as γ in Fig. 5 to be classified as within specification limits, while, they are actually out-of-specification in the correlated case (i.e., such points are unlikely to occur since they lie outside the elliptical equal probability contour). However, since it is most important to find all regions which line outside the specifications, i.e., outside square ABCD in Fig. 5, the inability to detect points such as γ is critical only when no points lie outside the square ABCD, or when two unmeasured regions are equivalent in terms of estimated P y(s,T) . During the extensive simulation, verification and validation studies carried out in this paper the aforementioned situation did not occur.
The probability that each point in a region r is out of design specification is used to calculate the entropy I r [54] using Eq. (16).
Sr (16) where N Sr represents the total number of points in region r and P y(s,T) is calculated using Eq. (15). The entropy I r represents the expected information gain which measures the information obtained by measuring a given region r ∈ {R|R = 1, 2, …, Number of unmeasured regions} in terms of improving the detection of out-of-specification points on the whole part. Expected information gain calculated using Eq. (16) represents the gain obtained for a system with independent random variables and can be considered as an upper bound to the expected information gain of a system with correlated independent variable [54]. The unmeasured region with maximum expected information gain is chosen as the candidate for next measurement using Eq. (17).
= I Next region argmax( ) r R r (17) Following the measurement of the selected region, updated complete part deviations are obtained using Eq. (9) for the current part or using Eq. (13) for the next part.

Measurement stopping criteria
The measurement stopping criteria enable to ascertain the number of partial measurements for a part in a given time step. The stopping criteria is a function of both the assembly system CT and the measurement coverage criteria. For instance, CT can constrain the number of partial measurements depending on the time available for measurements, or, the measurement can be stopped before the maximum number partial measurements possible within the CT if (i) a given level of confidence in deviations pattern prediction is attained, or (ii) a part geometric and dimensional fault and its root cause have been identified. The number of measurements can also be predetermined based on the part tolerance requirements and error of prediction obtained in the training data. In this paper to facilitate CT reduction estimation and model sensitivity analysis, the part is measured adaptively until all regions are covered.

Industrial case study: automotive door inner
The proposed STAS methodology for prediction of part surface deviation from partial measurements is demonstrated, verified and validated on a sport utility vehicle door inner sheet-metal part. The door inner is a part of the door subassembly consisting of window channel; halo; and hinge, latch and seat belt reinforcements. The quality of joining/fastening of door inner with other parts belonging to the subassembly is affected by geometric and dimensional variations in the input parts. For instance, remote laser welding is an emerging joining process with stringent part-to-part gap requirements of (i) 0.05 −−0.3 mm for galvanized steel sheet metal parts [55,56] or (ii) from 0 mm gap up to 50% of top parts for aluminium components [57,58], which are affected by part geometric variations. The faults generated by geometric and dimensional variations propagate in the assembly line causing quality defects such as gap and flushness issues between the door closure panel and automotive body-in-white openings. Application of the proposed STAS methodology enables the reduction of MTTD of such faults during the part fabrication process at stamping-press line (near line placement as shown in Fig. 1) or during the door assembly process.
The true manufactured part deviations data obtained from complete part measurements displayed limited deviation patterns, therefore to test the effectiveness of the proposed STAS methodology and to perform model sensitivity analysis, additional part deviation patterns were simulated using the VRM software [59,60]. A detailed explanation of the data generation methodology is provided in A.1. Seven batches of data were simulated in total with each batch consisting of 50 part deviation patterns, one batch was utilised for testing and six batches were utilised for training data sensitivity analysis (further increase in training batch size did not cause a significant change in average MAE of prediction) as described in Section 5.2.1.
The case study in this paper is conducted in three phases, firstly, the STAS methodology is demonstrated on the door inner part using the simulated part deviations data. Secondly, the following model performance and comparison analysis are conducted: (i) sensitivity analysis of model parameters, (ii) comparison of the proposed methodology with state-of-art spatial kriging methodologies [17,18,30] modified to be applicable to CoP measurements, and (iii) demonstration of the effectiveness of adaptive region selection methodology proposed in this paper against deterministic or fixed measurement sequences. Finally, the proposed methodology is verified and validated on the true part deviations data obtained from measurement of 33 manufactured door inner parts. A detailed analysis of the aforementioned studies beginning with the illustration of methodology on door inner part is presented below.

Initial model development
The simulated part deviations data in the case study are generated by imparting auto-correlated deviations to part fixture components, a few specimen part deviations is illustrated in Fig. 6. The optimum number of batches of training deviation data and the number of orthogonal bases are estimated by sensitivity study described in Section 5.2.1. Accordingly, in this case study, 100 part deviation patterns corresponding to two batches of the simulation data was utilised for training the model and K was set to 45 eigen bases. Since the part deviation data are obtained from simulation, the measurement error is assumed to be zero. The system matrix Φ is estimated by Singular Value Decomposition (SVD) of the training data. The state transition matrix H, and covariance matrices R and Q are calculated by method of moment estimators [49]. Key points ζ k are chosen uniformly from all nodes, by finding the points of intersection of the part and its voxalised bounding box with cubic voxels. For the case study presented in this section, a voxel size of 40 mm was chosen. This measure was decided after conducting a sensitivity study based on MAE of prediction and computational time of the model, for voxel sizes ranging between 20 to 50 mm in 5 mm increments. The part, is divided into 9 rectangular regions with each region covering approximately 11% of part surface area, as illustrated in Fig. 7. Following the initial model development, the proposed STAS methodology is applied to the test data consisting of 50 part deviation patterns.

State variables' update and part deviations prediction
The update of state variables to incorporate partial measurement data, and the prediction of complete part deviations is carried out in the following two stages: (a) Update of state variables (Ψ) and prediction of part deviations for the current part (T).
Following the measurement of a region of part T, data from a partial measurement Z (˜) are utilised to update state variablesˆusing Eq. (6) to reflect the data assimilation. The updated state variablesˆare used to predict deviation of the part surface from nominal using Eq. (9), and variance of the predicted deviations are obtained using Eq. (10). The updated predictions and their variance in turn are utilised to select the next region of part T to be measured, according to the criterion described in Section 4.3.
(b) Update of state variables and prediction of part deviations for the next part (T + 1).
When all regions of a part T are measured, the STAS methodology by utilizing temporal correlations present in the system can predict deviations of part T + 1 before any region of part T + 1 is measured (called zeroth prediction). The state variables for part T + 1 are predicted using Eq. (11) as described in Section 4.2.2, which in turn are utilised to make zeroth prediction using Eq. (13). The variance of the predicted deviations are obtained using Eq. (14). The first region of part T + 1 to be measured is selected using the zeroth prediction and its variance according to the criterion described in Section 4.3. An illustration of zeroth prediction along with with (T) th and (T + 1) th part deviation pattern is shown in Fig. 8.

Adaptive region selection
A key capability of the STAS methodology developed in this paper is the ability to adaptively select critical regions of the part to be measured without user intervention. In this paper, the criterion for adaptive region selection is chosen to be the identification of out-of-specification areas with-respect-to form tolerance requirements for profile of a surface with LSL and USL of a deviation of −1 and +1 mm from nominal, respectively.
The complete part deviations Ŷ and their variances estimated after each partial measurement as described in Section 5.1.2 are utilised to estimate P y(s,T) , the probability each mesh node of the part being out of specification using Eq. (15). Following the estimation of P y(s,T) , the expected information gain I r achieved by measuring a given region of the part is calculated using Eq. (16), and the region with maximum information gain is selected to be measured according to Eq. (17). The adaptive region selection for a given part is carried out until a stopping criterion is satisfied as described in Section 4.4.

Measurement stopping criteria
In this paper, a part is adaptively measured using the criterion described in Section 4.3 until all regions are covered to facilitate model sensitivity analysis studies. The Average MAE of prediction and 3-Sigma bound on error of prediction using the proposed STAS methodology for the test batch after each partial measurement is summarised in Table 2. The manually programmed robotic measurement plan to measure the entire part required 83 snapshots (Fig. A.15) and a CT of 510.5 s. The average CT for each partial measurement and percentage reduction compared to the CT for complete measurement is summarised in Table 2.   The average MAE and CT reduction presented in Table 2 can be utilised as guide to the stopping criteria of the STAS methodology for a given system. The number of partial measurements can be chosen based on the requirements of the process, for instance, if a MAE of 0.005 mm is acceptable in determining the component's conformance to a given form tolerance and the available CT for measurement larger than 296 s, the number of regions to be measured can be reduced to 3 from the total 9 required to measure the complete component, resulting in a measurement CT reduction of 42.2% compared to that required for the measurement of the complete part.

Model sensitivity analysis for number of batches of training data and eigen basis shapes utilised for model development
The developed adaptive spatio-temporal sampling model was analysed for sensitivity to two key model parameters: (1) number of batches of training data utilised for model development, and (2) number of eigen basis K for a given number of training batches. The first parameter gives the optimal training data required to develop the model, whereas the second determines the length of the state variable Ψ.
In this paper, different models are evaluated on the basis of MAE of prediction of the true deviation pattern. The model MAE as shown by Fig. 9 is unaffected by the number of eigen basis after two batches. Each training batch consists of 50 part deviation patterns, the change in MAE between two batches represents the change due to the increase in available training deviation patterns. The MAE plotted in Fig. 9 is 95% upper confidence bound on the mean of MAE after each partial measurement averaged for 50 part instances in the test data. The change in average MAE after a minimum of two batches of training data and 45 basis vectors, is insignificant for practical purposes. Based on this sensitivity analysis, the first 100 deviation patterns obtained from two batches of simulation are utilised to obtain model parameters with K representing the number of orthogonal basis chosen to be 45.

Comparison between state-of-art spatial prediction prediction and STAS methodologies
The state-of-art part deviation prediction methodologies typically applied to point data [17,18,30] are applied in the present study to CoP data by making the following changes: (1) using Full Independent Training Conditional (FITC) [61] to optimise the covariance function parameters for large data set of CoP data, and (2) selecting key points as measurement points. Additionally, to facilitate comparison, the measurement sequence obtained from the STAS methodology is used to predict part deviation using state-of-art spatial prediction techniques.
The MAE of predictions from both methodologies after each partial measurement averaged for 50 simulated deviation patterns in the test data is shown in Table 3, the STAS methodology is shown to be consistently better compared to the state-of-art spatial predictions with statistically significant lower average MAE values. A one tailed two sample t-test for each of the eight measurement cases (the case of complete measurement is same for both methodologies as the complete part is measured) rejected the null hypothesis that the average MAE of STAS methodology is greater than the average MAE of the state-of-art spatial prediction (p value ≈0).
The STAS approach, after a single measurement, is able to predict deviation pattern of regions far from the initial measurements. This ability can be attributed to the information carried by the state variable Ψ, which influences the entire part and helps in deviation pattern prediction with information gained from partial measurements. In contrast, the state-of-art spatial prediction has no information regarding deviations in regions far from the measurement area. An illustration of this is seen in Fig. 10 where the prediction of the true pattern (Fig. 8c) after three partial measurements (regions 3, 4 and 9) along with the error in prediction for both methodologies is shown.

Comparison between deterministic and adaptive measurement sequence
The effectiveness of the proposed adaptive region selection criterion in comparison with a deterministic measurement sequence is analysed in this section. Deterministic measurement sequence refers to measuring the same regions for all parts passing through the assembly line. The best and worst deterministic measurement sequence for the test batch consisting of 50 simulated deviation patterns is found based on the MAE of prediction. For instance, if three regions of a part can be measured due to constraints in cycle time, the fixed regions to be measured for all parts in a batch to obtain the best and worst performance are found through an exhaustive search of all the possible three region combinations by comparing the average MAE in predicting the original pattern after three measurements. This exhaustive search is performed individually for all 8 partial measurement options, resulting in 16 sequences, a best and a worst sequence for each of the 8 partial measurements, which are compared with the adaptive measurement sequence obtained through the proposed methodology.
However, these best and worst sequences can only be found retrospectively when deviation pattern for all parts in the batch are known, such sequences cannot be estimated during actual production. Thus, a deterministic sequence in real implementation could be anywhere between the best and worst case scenarios. The averaged MAE after each partial measurement using the proposed methodology along with the same metric for the best and worst deterministic sequences is plotted in Fig. 11.
A two sample t-test with the null hypothesis (H 0 ) of equal means  was performed to test a significant difference between average MAE for best deterministic and the adaptive sequences. The test failed to reject H 0 for partial measurements where one and two regions are measured.
For partial measurements with number of regions measured from 3 to 8 the null of hypothesis equal means was rejected, however the maximum difference in average MAE was 2 μm which is acceptable for deviation pattern recognition and small compared to the measurement error of 3D-optical scanners which is approximately 25 μm. The ability to obtain MAE similar to the best deterministic case, in-line adaptively makes the adaptive region selection methodology desirable.

Verification and Validation (V&V) of the proposed methodology
The methodology was applied to measurement data from a batch consisting of 33 manufactured door inner parts by using 15 measurements to train the model and 18 to validate the methodology. The experimental setup consisted of a Hexagon WLS400A 3D-optical scanner mounted on an ABB industrial robot, with the door inner part loaded on a reconfigurable alufix fixture as shown in Fig. 12. The average MAE after each partial measurement and its 95% confidence bound for the V&V test cases is shown Fig. 13. Higher MAE is consistent with simulation studies in Section 5.2.1 and should decrease with more training data.

Conclusions and future research
A novel STAS methodology to enable efficient in-line CoP measurements for industrial sized sheet metal parts is developed in this paper. It utilises spatio-temporal correlations present in the assembly system to enable prediction of deviation pattern of complete part from partial measurements. The necessity of partial measurements arises due to the high cycle time of 3D-surface scanners compared to the assembly system cycle time. The methodology effectively utilises the historical measurements or simulation data to model spatio-temporal correlations and optimise measurement coverage by predicting deviation pattern of the entire part from partial measurements. As illustrated by the industrial case study, the methodology has the ability to: (i) work with little data and has an optimum requirement of approximately 100 complete parts deviations obtained through measurement, simulation or historical data; (ii) reduce the measurement CT by 42.2% compared to that required for complete part measurement; and, (iii) achieve MAE comparable to the best deterministic measurement sequences. The part deviations predicted at each time step can be used to enable in-line inspection to ensure conformance to GD&T specifications in production systems with CT greater than that presented in Table 2. For production systems with CT lower than that presented in Table 2, the STAS methodology enables higher inspection frequency thus reducing MTTD of quality defects.
However, methodology has the following limitations: (i) it can only predict global deviation patterns and cannot detect local deformations, (ii) model parameters and the sensitivity study results are system specific and need to be performed for each new system configuration, and (iii) effective prediction of deviation patterns depends on the training data used to model the system parameters. Thus, the data utilised to  model the system should closely represent the true system. Future research can (i) incorporate the following aspects for the adaptive region selection: robot kinematics, and, part regions considering the quality of snapshot which depends on ambient conditions, surface reflectivity and sensor position [32]; (ii) mathematically characterise the quantity of training data required for optimal performance of the STAS methodology for a given system; and (iii) quantify the effect of non-ideal part variation on the performance of the STAS methodology [62].

Data statement
The source-code for the developed STAS methodology can be accessed online at [63].

Declaration of Competing Interest
The authors declare that there is no conflict of interest.
The simulated part deviation data are generated by imparting auto-correlated deviations to part fixture components, i.e., clamps, in the surface normal direction. A schematic representation of the door inner along with its fixture layout utilised to generate the part deviation data is shown in Fig. A.14. Temporal correlation in the deviation pattern is simulated using a first order Auto-Regressive process (AR(1)) defined according to Eq. (A.1).
where d(T) and d(T − 1) represent deviation of a clamp in time T and T − 1 respectively; γ is a zero mean Gaussian random variable with variance 0.25; constants a 0 and a 1 determine the mean and autocorrelation strength of the AR(1) process and are set to 0 and 0.7, respectively. Spatial correlation is exhibited by the part through physical deformation due to the imparted locator variation. Higher order temporal correlations can be modelled by modifying the state variables and the system matrices in Eqs. (5) and (4). Clamp deviations generated using Eq. (A.1) are applied as boundary conditions to the door inner in the surface normal direction. Part deviation pattern corresponding to each of the clamp's deviation from ideal position is generated through Finite Element Analysis (FEM) using VRM software [59]. Each batch of simulated data consists of 50 complete part deviations corresponding to 50 auto-correlated locator deviations. The entire simulated data consists of one test batch and six training batches, as further increase in training batch size did not cause a significant change in average MAE of prediction.   15 illustrates the 83 regions (i.e., the snapshots of the 3D-optical scanner) of the door inner part which takes into account the visibility and surface reflectivity constraints. However, since the aim of this paper is to predict complete part deviations with partial measurements, we divide the part into large regions with no overlaps (Fig. 7) without considering the visibility and surface reflectivity constraints.