A Registration-free approach for Statistical Process Control of 3D scanned objects via FEM
AA Registration-free approach for Statistical Process Control of 3Dscanned objects via FEM
Xueqi ZhaoEnrique del Castillo ∗ Engineering Statistics and Machine Learning LaboratoryDepartment of Industrial and Manufacturing Engineering and Dept. of StatisticsThe Pennsylvania State University, University Park, PA 16802, USA
January 8, 2021
Abstract
Recent work in on-line Statistical Process Control (SPC) of manufactured 3-dimensional (3-D) objectshas been proposed based on the estimation of the spectrum of the Laplace-Beltrami (LB) operator, adifferential operator that encodes the geometrical features of a manifold and is widely used in MachineLearning (i.e., Manifold Learning). The resulting spectra are an intrinsic geometrical feature of each part,and thus can be compared between parts avoiding the part to part registration (or “part localization”)pre-processing or the need for equal size meshes, characteristics which are required in previous approachesfor SPC of 3D parts. The recent spectral SPC methods, however, are limited to monitoring surface datafrom objects such that the scanned meshes have no boundaries, holes or missing portions. In this paper weextend spectral methods by first considering a more accurate and general estimator of the LB spectrumthat is obtained by application of Finite Element Methods (FEM) to the solution of Helmholtz’s equationwith boundaries. It is shown how the new spectral FEM approach, while it retains the advantages of notrequiring part localization/registration or equal size datasets scanned from each part, it provides moreaccurate spectrum estimates, which results in faster detection of out of control conditions than earliermethods, can be applied to both mesh or volumetric (solid) scans, and furthermore, it is shown howit can be applied to partial scans that result in open meshes (surface or volumetric) with boundaries,increasing the practical applicability of the methods. The present work brings SPC methods closer tocontemporary research in Computer Graphics and Manifold Learning. MATLAB code that reproducesthe examples of this paper is provided in the supplementary materials.Keywords: Manifold Learning; Part localization; Noncontact sensor; Spectral methods; Helmholtz equation ∗ Corresponding author. Dr. Castillo is Distinguished Professor of Industrial & Manufacturing Engineering and Professor ofStatistics. e-mail: [email protected] a r X i v : . [ s t a t . A P ] J a n Introduction
Modern digital manufacturing deals not only with larger metrology data sets but also more complex datatypes that have various structures. Point cloud, mesh, and voxel datasets are some of the most common typesof data acquired with non-contact sensors. In the area of Quality Control of manufactured parts, methods forthe assessment of the quality in a sequence of manufactured parts have relied on point measurements of eachpart, typically acquired with a coordinate measurement machine (CMM), resulting in point cloud or meshdatasets with the exact same number of points from part to part, so that a point to point correspondencebetween parts can be established via registration or superposition. While there is considerable work on3D part inspection, relatively little work has taken place in on-line inspection or on-line Statistical ProcessControl (SPC) of manufactured parts based on non-contact sensor data (Babu et al., 2017), and both fieldshave relied on different types of registration of the parts as a first step or “pre-processing” of the scans ormeasurements obtained from a sequence of parts or between a scanned part and its CAD model (Wells et al.,2013; Huang et al., 2018; Zang & Qiu, 2018).In recent work, Zhao & Del Castillo (2020) (hereafter, ZD) introduce a novel SPC method based onmonitoring the spectrum of the Laplace-Beltrami (LB) operator estimated from the scans of each part. Asdiscussed below, the LB operator, present in both the heat and wave partial differential equations (PDEs)and widely used in Machine Learning (more specifically, in the area of Manifold Learning), codifies thegeometrical information of a manifold (in this case, surfaces or solids). The spectrum of the LB operatoris intrinsic, that is, it does not depend on the coordinates of the ambient space in which the object isembedded, and hence it can be compared between parts without any registration, totally avoiding the partlocalization problem. ZD’s spectral SPC method is however restricted to meshes modeling the surface of3-dimensional objects that must be “closed”, that is, must have no boundaries due to missing parts or holes.More precisely, we define a closed mesh as a mesh of discrete elements (flat triangles for surface data andcuboidal voxels for volumetric data) such that all the boundaries of all elements in the mesh are in contactwith neighboring elements. Then, we simply define an open mesh as a mesh that is not closed, hence thereare some boundary elements in it.In the present paper we extend the spectral SPC method in ZD in three major directions: 1) we estimatethe spectrum of the LB operator with Finite Element Methods, resulting in a considerably more accurateestimator of the analytical LB spectrum than the Localized mesh Laplacian used by ZD; 2) we show how thenew FEM spectral SPC method can be applied to both surface data (triangulation meshes) and volumetricdata (i.e., voxel data) of parts, as acquired by either a range sensor or Computed Tomography (CT) scanner,respectively; 3) we demonstrate how by solving Helmholtz equation as a boundary value problem, the newspectral FEM method can be applied to scans that result in open surface or volumetric meshes. Given thatscans are often open due to part regions that are inaccessible to the scanner (a problem which we will refer toas occlusion), the ability to deal with incomplete or open meshes greatly increases the practical applicabilityof the proposed spectral methods. The focus therefore is on the more accurate and versatile estimation ofthe LB spectrum and the consideration of voxel and open meshes; for the specific SPC methods once thespectrum is estimated we follow those used by ZD.As a preview of the power of the spectral FEM SPC methods presented below to detect defects onmanufactured parts relative to earlier methods, consider the prototype part studied by ZD, shown in Figure1, which displays the CAD model followed by three similar parts with different defects, two with small“chipped” errors a corner, and one with a small “protrusion” in one of the top “teeth” of the part. Theseare triangulation meshes with no boundaries, but as mentioned, the methods presented in this paper applyequally to open triangulation meshes or open volumetric (solid) meshes. We can compare the power of thedifferent methods to distinguish between the four types of parts (the non-defective and the three defective)by performing Multidimensional Scaling (Borg & Groenen, 2005) on the first 15 eigenvalues of the LBspectrum computed, for 10 different parts of each type, with the FEM methods presented in this paperand the method used earlier by ZD (based on the localized mesh Laplacian of Li et al. (2015)). Figure 2 isthe 3D multidimensional scaling plots of the first 15 LB eigenvalues for the four different part types whendifferent Laplacian discretizations are used. Each part type has 10 simulated realizations, resulting in 10points of the same color in the abstract 3D Euclidean space. Since multidimensional scaling tries to preservethe between-object distances, this visualization reveals the reason why the FEM methods we propose are2 igure 1: Prototype parts used by Zhao & Del Castillo (2020) to demonstrate their spectral SPC scheme. The left most partis a perfect part (CAD model), while the other three parts on the right contain different types of defects (“chipped” corners ora “protrusion”). All meshes are triangulations without boundaries. This prototype type is used further below for the SPC runlength analysis of the new methods in Table 1.Figure 2: 3D multidimensional scaling plot of the first 15 eigenvalues of the LB spectrum for the four part types in Figure 1estimated with the proposed FEM methods in this paper (linear and cubic FEM) and with the Localized Laplacian methoddue to Li et al. (2015) used by Zhao & Del Castillo (2020) in their spectral SPC approach. Ten realizations with noise of eachpart were simulated. The tight clustering of the FEM LB spectra in 4 groups compared to the lack of separation in the spectraof the Localized Laplacian explains why FEM LB-based Statistical Process Control methods can detect defects much faster,while maintaining a low false detection rate. more powerful in differentiating subtle shape changes than the Li et al. (2015) Localized Laplacian, as theleading spectra calculated by the FEM methods are tightly clustered basing on the part type (shown asdifferent colors in the plot), while the LB spectra computed based on the Localized Laplacian is unable toseparate the different parts. Hence, if a SPC chart is to monitor the LB spectra as a feature or “profile”from each part, it will detect the defective parts considerably faster with the FEM LB methods presentedin the present paper.The rest of the paper is organized as follows. Section 2 reviews Differential Geometric notions neededin the sequence, in particular, the definition of the Laplace-Beltrami operator and its properties, which weuse in the rest of the paper. Section 3 presents the Finite Element Method (Garlekin approach) for solvinga boundary value Helmholtz PDE problem (which contains the LB operator and results from the spatialpart of both wave and heat PDEs) and estimates the spectrum of the LB operator in this way. WhileFinite Element Methods are well known in engineering and science, we present the peculiarities behind bothFEM formulations, for surface meshes and voxel data, in enough detail for readers to be able to reproduceour results. Finally, in section 4 we show the performance of an on-line or “Phase II” nonparametric SPCchart based on a sequence of estimated FEM LB spectra of parts, and demonstrate its performance againstprevious Laplacian methods and SPC methods that are based on registration (superposition) of the parts. Asis common in the field of SPC, we are concerned with both quick detection of out of control states (defects)and avoidance of false positives for as long as possible. We close with conclusions and further research insection 5. The supplementary materials provide MATLAB code that implement our FEM methods andreproduce the examples in the paper. 3
The Laplace-Beltrami operator and its spectrum
To introduce the Laplace-Beltrami operator, consider a parametric surface p ( u, v ) = ( x ( u, v ) , y ( u, v ) , z ( u, v )) (cid:48) ,( u, v ) ∈ D ⊂ R , where u = x , and v = x are local coordinates, defining M , a surface or Riemannian2-manifold (similar definitions will apply in the case of 3-manifolds or solids). Define the surface differentialvectors at p ( u, v ) as: p u = ∂ p ( u, v ) ∂u = (cid:18) ∂x ( u, v ) ∂u , ∂y ( u, v ) ∂u , ∂z ( u, v ) ∂u (cid:19) (cid:48) and p v = ∂ p ( u, v ) ∂v = (cid:18) ∂x ( u, v ) ∂v , ∂y ( u, v ) ∂v , ∂z ( u, v ) ∂v (cid:19) (cid:48) . Define next g = (cid:104) p u , p u (cid:105) g = (cid:104) p u , p v (cid:105) g = (cid:104) p v , p v (cid:105) (where (cid:104) , (cid:105) denotes the standard inner productin Euclidean space) and define the Riemannian metric tensor associated with the surface M , which definesan inner product on vectors tangent to M : (cid:104) w , w (cid:105) M = w T g w , g = (cid:18) g g g g (cid:19) .The metric g is induced by the ambient (Euclidean) space on the surface M , but note it is intrinsic , i.e.,it does not rely on the coordinates of the ambient space. Intrinsic geometrical properties, those exclusivelybased on the metric tensor, are invariant with respect to rigid transformations. Hence, Zhao & Del Castillo(2020)’s idea was to compute an intrinsic differential operator based on scanner data that models geometricalfeatures of an object, because being invariant with respect to rigid transformations it could be used forinspection or statistical quality control without having to register the scanned parts.The Laplace-Beltrami operator extends the notion of the Laplacian of a function defined on flat (Eu-clidean) space to functions defined on curved space or Manifolds. Recall the Laplacian of a twice differentiablefunction f : R n → R is minus the divergence of its gradient field: ∆ f = − div ∇ f = − (cid:80) ni =1 ∂ f∂x i and isevidently a measure of curvature of f at a point x . Similarly, for a function f : M → R , the Laplace-Beltrami (LB) operator is defined as ∆ M f = − div M ∇ M f , where div M is the divergence taken on M . This is indeedan intrinsic measure of curvature of f defined at a point on the manifold, and in contrast to the Laplacianof a function defined on flat space, it encodes the curvature of the manifold itself as well. In general, appliedto a function f ( x , ..., x k ) ∈ C defined on a k -manifold M , the LB operator is:∆ M f = − (cid:112) det( g ) k (cid:88) j =1 ∂∂x j (cid:32)(cid:112) det( g ) k (cid:88) i =1 g ij ∂f∂x i (cid:33) (1)where g ij are the elements of g − and det( g ) is the determinant of the metric tensor. One importantproperty of the LB operator which aids in its interpretation is that, for a surface (2-manifold):∆ M p ( u, v ) = 2 H n ( u, v ) (2)and similarly for higher dimensional manifolds, where n ( u, v ) is the normal at the point p ( u, v ) on M and H is the mean curvature of M at p , which is the average of the maximum and minimum curvatures in anydirection on M from point p . The LB operator appears in both the heat and wave partial differential equations where the space of interestis a Riemannian manifold. In either case, simple separation of variables and consideration only of the spatialvariables results in the eigenproblem:∆ M f = λf, (3)called the Helmholtz partial differential equation, with an infinite number of eigenfunctions f : M → R (providing spatial or static solutions to both the heat and wave equations) and corresponding eigenvalues4 ∈ R . We consider solving Helmholtz’ equation subject to either Dirichlet boundary conditions ( f = 0 onthe boundary Γ of M ) or Neumann boundary conditions ( ∂f∂ n = ∇ f · n = 0 on Γ, where n is the outwardnormal vector to M and ∇ is the gradient operator). We point out that considering boundary conditionswill permit us to address the case of open or incomplete meshes or volumes, not considered by Zhao &Del Castillo (2020), whose spectral methods were limited to closed objects without a boundary.The collection of eigenvalues { λ i } ∞ i =0 (0 = λ ≤ λ ≤ λ , .... ) obtained from solving (3) is called the spectrum of the LB operator, which Zhao & Del Castillo (2020) proposed to compute numerically to monitorthe quality of discrete parts in manufacturing. In the particular case M ⊂ R , the eigenfunctions f ( u, v )satisfying the Helmholtz equation are often interpreted as the modes of vibration of a membrane or “drum”with resonances at frequencies λ i (Kac, 1966).Even though the analytical LB spectrum of very few 3D objects is known, one of them being the sphere(see below), the estimated LB spectra can be monitored via a multivariate SPC chart for any object scan(surface or volumetric), comparing the estimated spectrum of a new part produced under regular production(what in SPC is called “Phase II”) against the spectra estimated from parts obtained while the monitoringscheme was started up (what in SPC is called “Phase I”). As it will be shown, only the lower part ofthe spectrum is needed for part to part comparisons. Following the classical SPC paradigm (see, e.g.Montgomery (2020)) we assume the Phase I spectra were obtained while the process was in a state ofstatistical control. In contrast with Differential Geometry, we do not have an analytical expression for the object being modeledas a surface (2-manifold) or solid (3-manifold) M , and hence, the first task is to estimate the LB operatorfrom metrology data. In the surface case, we assume in this paper we have available triangulation mesh data(including the possibility of incomplete meshes with holes due to regions on the object that are unreachable tothe scanner, see below), consisting of a sample of points from the surface of the object and their adjacencyinformation, usually generated by built-in algorithms used by the sensor mechanism. For solid data, weassume we have available volumetric data acquired by a Computed Tomography (CT) scanner and pre-processed so the result is a set of 3D set of voxels obtained after application of reconstruction and edge-detection algorithms to find the boundaries of the object from the voxel attenuation values (Kruth et al.,2011).The analytic LB operator is a differential operator acting on a continuous function. Therefore, the firsttask in practice is to estimate a discrete version of the LB operator in the form of a matrix. Several suchdiscretizations exist in the literature, for instance, del Castillo & Zhao (2020) evaluated the performanceof the heat kernel based approximations (Belkin et al., 2008; Li et al., 2015) for surface data and latersuggested to use the Localized LB estimator of Li et al. (2015) for SPC applications using mesh data, due toits sparseness (Zhao & Del Castillo, 2020). See Wardetzky et al. (2007); Patan´e (2016) for a comprehensivereview of discrete LB estimators and their convergence properties. In this article, we use Finite ElementMethods proposed by Reuter et al. (2006, 2009) who use them in medical applications. FEM methods canbe used with either surface or volumetric data (Reuter et al., 2007; Niethammer et al., 2007), and, as willbe shown below, provide a more accurate estimation of the true analytical LB spectrum and can easilyincorporate boundary conditions permitting the estimation of the LB spectrum on open meshes or solids, akey advantage over the Li et al. (2015) method used in Zhao & Del Castillo (2020).Given the very large literature on Finite Element Methods, we will only present next the details of theirapplication to the specific solution of the Helmholtz equation (3) using the classical Galerkin variationalformulation, from which the spectrum can be obtained, for both surface mesh and volumetric metrologydata. Additional details can be found in the Appendix A.5 .1 Galerkin Variational Formulation for the Helmholtz Boundary Value Prob-lem The classical presentation of FEM methods for the solution of a partial differential equation (PDE) startswith the so-called weak, variational, or Galerkin formulation of the problem (see e.g. Le Dret & Lucquin(2016)) which for the Helmholtz equation (3) we are concerned with in this paper consists in finding afunction f ∈ V defined on the manifold M that satisfies the equation: − (cid:90) M ∇ f · ∇ φ dV = λ (cid:90) M φf dV ∀ φ ∈ V (4)where dV is either the surface element on a 2-dimensional manifold (surface) or a volume element in a 3-manifold M . The weak form (4) is arrived at by considering functions f that satisfy the boundary conditions(see Appendix A for some notes about the derivation of the variational form). The space of functions V is called the trial space and V is called the test space, and for the Helmholtz equation, they are naturallydefined to be both Sobolev spaces H ( M ) = { u : ||∇ u || + || u || < ∞} where || · || is the L norm. Inorder to account for the curvature of the manifold, a subtlety about the dot product ∇ φ · ∇ f in the firstintegral, called the first differential parameter of Beltrami (Kreyszig (1991), p.230) and sometimes denotedby ∇ ( φ, f ), is that it must be defined as: ∇ φ · ∇ f = ∂φ (cid:48) g − ∂f = (cid:88) i,j g ij ∂ i φ∂ j f which is an inner product between the gradients of φ and f each expressed in local basis form (see Lee(2018), p. 27) that is, ∇ f = (cid:80) i,j g ij ∂f∂x j e i = (cid:80) i,j g ij ∂ j f e i (where { e i } di =1 is a local Euclidean orthonormalbasis) and likewise for ∇ φ , a consideration necessary given that gradient vectors are covariant .The customary way to summarize the weak or variational form is based on defining the inner products: (cid:104)∇ f, ∇ φ (cid:105) = (cid:90) M ∇ f · ∇ φ dV and (cid:104) f, φ (cid:105) = (cid:90) M φf dV so that (4) is usually written in compact form as: − (cid:104)∇ f, ∇ φ (cid:105) = λ (cid:104) f, φ (cid:105) . (5)It can be shown that if f ∈ V is a solution of the Helmholtz equation (3) then it must satisfy equation(5) for all φ ∈ V . The reverse implication is also true, it can be shown that if f ∈ H satisfies (5) for all φ ∈ V then it is a solution of the Helmholtz PDE (3), see (Le Dret & Lucquin, 2016, propositions 4.1 and4.2 respectively). Equation (4) can be motivated also as being the Euler-Lagrange equation of either anenergy or a least squares error functional, see Appendix A. Solutions f so obtained are eigenfunctions forthe Helmholtz equation for the corresponding eigenvalue λ .Rather than directly solving an infinite dimensional problem in functional space, the FEM strategyconsists in solving a finite dimensional problem by approximating the solution by f N , a linear combinationof N known “shape” functions h i , whose coefficients u i must be determined: f N = u h + u h + · · · + u N h N = N (cid:88) m =1 u m h m (6)The N linearly independent shape functions h , h , · · · , h N are selected to form a basis for the space ofapproximate solutions and to be such that each h i function has local support only over a single discretefinite element in which the space M (surface or volume) is then partitioned. The local form of the gradient can be expressed in vector form as g − ∇ f and hence Beltrami’s first differential parametercan be written in terms of the inner product defined on M as (cid:10) g − ∇ φ, g − ∇ f (cid:11) M . N different test functions φ to solve for the N coefficients. In the Garlekin method, the test functions are exactly the same as the shape functions, thuswe also substitute φ for each h i ( i = 1 , ..., N ) in (4), which results in the N equations: − (cid:90) M (cid:88) g ij ∂ i h l ∂ j (cid:32) N (cid:88) m =1 u m h m (cid:33) dV = λ (cid:90) M h l (cid:32) N (cid:88) m =1 u m h m (cid:33) dV, l = 1 , , · · · , N N (cid:88) m =1 u m − (cid:90) M (cid:88) i,j g ij ∂ i h l ∂ j h m dV = λ N (cid:88) m =1 u m (cid:18)(cid:90) M h l h m dV (cid:19) , l = 1 , , · · · , N (7)This system of equations can be written as a generalized eigenvalue problem in the matrix form: AU = λ BU (8)where A and B are N -by- N Gram matrices with entries: A lm = − (cid:104)∇ h l , ∇ h m (cid:105) = − (cid:88) i,j g ij (cid:90) M ∂ i h l ∂ j h m dV and B lm = (cid:104) h l , h m (cid:105) (9)and U is the vector ( u , u , · · · , u n ) T . Once the shape functions are chosen, both eigenvalues, which arethe Laplace-Beltrami eigenvalues, and eigenvectors, which give the Laplace-Beltrami eigenfunctions, can beeasily computed by solving the eigenproblem (8). Though f N is theoretically only an approximation of f ,given that a finite basis can not span the whole space for infinite dimensional functions, it can become anexact solution when the manifold M is discretized based on a mesh of finite elements in which case functionson M are reduced to vectors of dimension N , the mesh size. This justifies the choice of the basis size in theprevious step. Now we discuss how to choose the shape functions when the scan of a part has generated surface datain the form of a triangular mesh. As can be seen from (6), the N shape functions compose a basis of thesolution space, R N , with N being the mesh size. The simplest way to ensure linear independency is to use N indicator functions, one for each nodal point. Thus, for l = 1 , , ..., N , the l th form function, h l , takes value1 at the l th nodal point and 0 at the other points. Its function values elsewhere on M will be determinedlater based on the properties of h l . A popular choice for the shape functions are piecewise polynomials, forexample: • h l is linear over each finite element h l ( u, v ) = c l, + c l, u + c l, v (10) • h l is quadratic over each finite element h l ( u, v ) = c l, + c l, u + c l, v + c l, u + c l, uv + c l, v (11) • h l is cubic over each finite element h l ( u, v ) = c l, + c l, u + c l, v + c l, u + c l, uv + c l, v + c l, u + c l, u v + c l, uv + c l, v (12)where we have used x = u, x = v as the local coordinates on the manifold M , and the finite elements aresimply triangles as we focus on surface triangulations. Note h l is an indicator function, so for each triangle, h l = 1 only when the triangle has point l as one of its vertices, and h l = 0 otherwise.Since each triangle is associated with three nodal points at its vertices where the function values h l ( u, v )are known, either 1 or 0 depending on the indices of the shape function under consideration and the vertex,7 igure 3: Nodal points used to construct linear shape functions (left), quadratic shape functions (middle), and cubic shapefunctions (right). no additional information is needed to uniquely determine the three coefficients c l,j , j = 1 , , A and B matrices for surface data (triangu-lations) We take the linear FEM as an example to show the construction of the A and B matrices. First let usconsider the simplest case, when the mesh only contains one flat triangle as shown in Figure 3 (left) andthe Euclidean coordinates are P (0 , , , P (0 , , , P (1 , , p ( u, v ) = ( u, v, , ≤ u ≤ , ≤ v ≤ , u + v ≤ u, v ) are the surface coordinates which coincide with ( x, y ) in this case. Using this local parametriza-tion, P = p (0 , , P = p (0 , , P = p (1 , (cid:40) ∂ p ( u,v ) ∂u = (1 , , ∂ p ( u,v ) ∂v = (0 , , ⇒ g = ∂ p ( u,v ) ∂u · ∂ p ( u,v ) ∂u = 1 g = ∂ p ( u,v ) ∂u · ∂ p ( u,v ) ∂v = 0 g = ∂ p ( u,v ) ∂v · ∂ p ( u,v ) ∂v = 1 ⇒ g = (cid:18) (cid:19) = I (14)which is not surprising since the triangle is flat. Since there are three points, we need three shape functions h , h , and h . Taking the linear case, h ( u, v ) = c , + c , u + c , v as an example, which takes value 1only at P and 0 elsewhere, we obtain the following system of equations h ( P ) = h (0 ,
0) = c , = 1 h ( P ) = h (0 ,
1) = c , + c , = 0 h ( P ) = h (1 ,
0) = c , + c , = 0 ⇒ c , = 1 c , = − c , = − ⇒ h ( u, v ) = 1 − u − v (15)Similarly, we can solve for h and h : h ( u, v ) = v, h ( u, v ) = u (16)Figure 4 plots some of the shape functions obtained in this way. As the degree of the polynomials increases,the shape functions become more flexible, as expected. Once the analytical expressions for the shapefunctions are known, entries in matrices A and B can be easily calculated.8 igure 4: Plots of the shape functions on a unit triangle: h in the linear case (left), h in the quadratic case (middle), and h in the cubic case (right). The z-axis represents function values and lighter colors indicate higher values. Consider now the case of three points with nontrivial coordinates: P ( x , y , z ) , P ( x , y , z ) , P ( x , y , z ).This can be converted to the previous case by modifying only the surface parametrization: p ( u, v ) = ( x + u ( x − x ) + v ( x − x ) , y + u ( y − y ) + v ( y − y ) , z + u ( z − z ) + v ( z − z )) , ≤ u ≤ , ≤ v ≤ , u + v ≤ P = p (0 , P = p (0 , P = p (1 , h , h , and h remain the same as in (15) and (16), so there is no need to reevaluatethe integrals in (9). On the other hand, the metric tensor g changes with the parametrization: (cid:40) ∂ p ( u,v ) ∂u = ( x − x , y − y , z − z ) ∂ p ( u,v ) ∂v = ( x − x , y − y , z − z ) ⇒ g = (cid:18) (cid:107) P − P (cid:107) ( P − P ) · ( P − P )( P − P ) · ( P − P ) (cid:107) P − P (cid:107) (cid:19) (18)This affects the entries of the A and B matrices in two ways. First, obviously g ij differs from triangle totriangle. Secondly, the surface area element in the integrals changes too, since dV = (cid:112) det( g ) dudv . So eq(9) becomes A lm = − (cid:90) M (cid:88) i,j g ij ∂ i h l ∂ j h m (cid:112) det( g ) dudv = − (cid:112) det( g ) (cid:18) g (cid:90) M ∂ u h l ∂ u h m dudv + g (cid:90) M ∂ u h l ∂ v h m dudv + g (cid:90) M ∂ v h l ∂ u h m dudv + g (cid:90) M ∂ v h l ∂ v h m dudv (cid:19) B lm = (cid:112) det( g ) (cid:90) M h l h m dudv. (19)Finally we consider the case of meshes consisting of an arbitrary number of connected triangles. Note thatthe ( l, m )th elements in matrices A and B require ∂ i h l ∂ j h m , and h l h m respectively, which are nonzeroonly when points l and m are connected by an edge and therefore appear in the same triangle(s). Thisimplies that we can process the mesh triangle by triangle, and only fill in the entries of A and B as needed.Furthermore, as we discussed above, for each triangle, only the metric tensor g (18) needs to be recalculated.When an edge connecting two points, say l and m , is not on the boundary, it will be included in two adjacenttriangles, each of which gives a value for a lm and b lm . In this case, the two different values of a lm (or b lm ) arethe integrals of ∂ i h l ∂ j h m (or h l h m ) evaluated in the two individual triangles, respectively, and thus shouldbe added up. In addition to their use for computing the LB spectrum for 2-dimensional surface data, the FEM methodscan be easily extended to the voxel or volumetric data case. The shape functions have three variables now9 igure 5: Nodal points and plots of the shape functions on a voxel. From left to right: a) Nodal points with two points labeledin their Euclidean coordinates, P (1 , ,
1) and P (2 / , , h corresponding to P . c) Cubic h corresponding to P . d) Cubic h corresponding to P . Lighter colors indicates higher function values. The graphs show how shape functions ensure h l ( P m ) = 0for l (cid:54) = m . with the increased dimension, and as suggested by Reuter et al. (2007) we use the trilinear function andthe cubic function of the serendipity family (Arnold & Awanou, 2011) for the linear FEM and cubic FEM,respectively: h l ( u, v, w ) = c l, + c l, u + c l, v + c l, w + c l, uv + c l, uw + c l, vw + c l, uvwh l ( u, v, w ) = c l, + c l, u + c l, v + c l, w + c l, uv + c l, uw + c l, vw + c l, uvw + c l, u + c l, u v + c l, u w + c l, u vw + c l, v + c l, v u + c l, v w + c l, v uw + c l, w + c l, w u + c l, w v + c l, w uv + c l, u + c l, u v + c l, u w + c l, u vw + c l, v + c l, v u + c l, v w + c l, v uw + c l, w + c l, w u + c l, w v + c l, w uv (20)Similar to the mesh case, a linear shape function uses the 8 vertices of each finite element, a voxel in thiscase, to uniquely determine its 8 coefficients, while a cubic shape function needs 24 more nodal points, 2trisection nodes on each of the 12 edges, together with the original 8 vertices to provide a total of 32 degreesof freedom (first graph in Figure 5). Each shape function is the indicator function of a corresponding node,and each is again a piecewise polynomial as before. Figure 5 shows a plot of fitted linear and cubic shapefunctions in a voxel, where lighter colors indicate higher function values. Since the voxel and the nodalpoints we use are highly symmetric, the linear shape functions have only one pattern, as in the secondgraph, while the cubic shape functions have two patterns, depending on whether the node is on a corner oran edge, both of which are shown in the last two graphs. Note the color scales differ in different graphs for amore detailed illustration of how function values vary within each object, but it can inferred that whichevercolor that occurs on the edges corresponds to a function value of 0 on that particular graph.The entries of the metric tensor g for voxels are simpler than for triangles. A voxel is essentially a cuboidof fixed dimensions and can be parametrized by: p ( u, v, w ) = ( x + s u, y + s v, z + s w ) , ≤ u ≤ , ≤ v ≤ , ≤ w ≤ x, y, z ) are the coordinates of the cuboid vertex that is closest to the origin, and s , s , and s arethe three edge lengths of the cuboid, respectively. Then the metric tensor is ∂ p ( u,v,w ) ∂u = ( s , , ∂ p ( u,v,w ) ∂v = (0 , s , ∂ p ( u,v,w ) ∂w = (0 , , s ) ⇒ g = ∂ p ( u,v,w ) ∂u · ∂ p ( u,v,w ) ∂u = s g = ∂ p ( u,v,w ) ∂u · ∂ p ( u,v,w ) ∂v = 0 g = ∂ p ( u,v,w ) ∂u · ∂ p ( u,v,w ) ∂w = 0 g = ∂ p ( u,v,w ) ∂v · ∂ p ( u,v,w ) ∂v = s g = ∂ p ( u,v,w ) ∂v · ∂ p ( u,v,w ) ∂w = 0 g = ∂ p ( u,v,w ) ∂w · ∂ p ( u,v,w ) ∂w = s ⇒ g = s s
00 0 s (22)10 igure 6: The computational time in seconds versus the mesh size N for computing the first 50 LB eigenvalues of a unit sphere.Experiments run on a machine with 3 GHz Quad-Core Intel Core i5 and 8G RAM. Note g is independent of ( x, y, z ), that is, independent of the location of the voxel, and depends only on thesize of the voxel, which is constant. Therefore, we can calculate only once the matrix elements a lm and b lm in a voxel and fill in all the other entries in A and B computing (19) by simply looking up the global nodalpoint indices (each nodal point has an index local within each voxel and a global index within the wholevolume of voxels). As can be seen in (9), b lm is nonzero if and only if both h l and h m are not constantly zero over at leastone finite element (a triangle or a voxel), which happens only when point l and m appear in the samefinite element. This indicates that matrix B is sparse. For example, in the surface linear FEM method,the number of nonzero elements along row l or column l equals the degree, or number of neighbors, ofpoint l . Similarly, a lm is nonzero if and only if both h l and h m are not constant over at least one finiteelement. Since h l is an indicator function and piecewise polynomial over each finite element, it is zero forelements not associated with point l , and it is never a constant function for elements associated with point l , l = 1 , , · · · , N . Therefore, when h l is non-constant (has a gradient) it is when it is nonzero, thus a lm and b lm are zero or nonzero at the same time. In other words, matrix A is sparse as well and has the samenonzero structure as matrix B .Another immediate property drawn from (9) is that both matrices A and B , being Gram matrices,are symmetric, which assures the estimated Laplace-Beltrami spectrum is real. There are computationalbenefits of solving the generalized eigenvalue problem (8) for symmetric and sparse matrices. For example,the Arnoldi algorithm has a typical computational complexity of O ( kN ) to solve for the first k eigenvalues,where N is the matrix size (Zhao & Del Castillo, 2020). This can also be seen from Figure 6, where thecomputational time for finding a fixed number of LB eigenvalues is linear in N , the mesh size, in both linearand cubic FEM methods. As expected, the cubic FEM takes longer as it works with larger matrices byadding additional nodes. Reuter et al. (2006) illustrate how the FEM methods appear to yield accurate results for the LB spectrum,with higher order shape functions resulting in more accurate estimations compared to the known analyticspectrum of some 3D objects. To verify this claim, we can compare the FEM LB spectrum with the methods11 igure 7: The leading spectrum of the FEM Laplacians versus that obtained from the heat kernel based Laplacian (Li et al.,2015) for a unit sphere under different conditions. The true (analytical) LB spectrum is also plotted for comparison. From leftto right: a-c) are noise-free spheres with increasing mesh size. Though all estimated spectra converge to the true spectrum asthe density of the mesh increases and approximates the manifold better, the FEM spectra are already very accurate with amesh size as small as 300 points. c-e) are the same sphere with increasing noise level, which show how all estimated spectraare affected by surface noise. A more accurately estimated LB spectrum will result in better detection properties when thespectrum is used for SPC or inspection purposes. used in Zhao & Del Castillo (2020) to compute the LB spectrum. Figure 7 displays the first 10 eigenvalues inthe LB spectrum of a unit sphere (one of the few 3D objects for which the spectrum is known analytically)using different LB estimation methods. Note how much more accurate the FEM methods already are fora mesh size of only 300 points compared with the Localized Mesh Laplacian of Li et al. (2015). As it canbe seen from the same figure, the FEM LB methods are robust to small surface noise and still sensitive toreflect surface changes caused by noise. Figure 8 also shows how both linear and cubic FEM LB spectraclosely approximate the analytical LB spectrum (Reuter, 2006; Helffer & Sundqvist, 2016) of several 3Dobjects in the voxel case. In the last two cases, since both the cube and cuboid can be represented exactlydespite the small voxel numbers, the cubic FEM spectra are extremely close to the corresponding analyticalLB spectrum. The difference between the first 50 LB eigenvalues obtained with the cubic FEM and thecorresponding eigenvalues from the analytical LB operator has a norm of only 0.2391 for the cube and 0.0317for the cuboid. A thicker line is used for the cubic spectrum so that it does not overlap with the red linerepresenting the true spectrum.In practice, convergence of the FEM method is achieved either by refining the mesh, i.e., decreasing themesh elements relative to the scanned object, or increasing the degree of the polynomial approximation.Reuter (2006) mentions how the convergence of the FEM method with shape functions of order p behavesasymptotically with an error of order O ( s p +1 ) as the largest mesh element size s goes to zero (here, s = max s i where s i is a measure of each element size, e.g., the length of the largest side or the radius of the largestinscribed circle in element i , see Ihlenburg (2006)). This fast convergence rate can be observed in Figure12 igure 8: The first 50 LB eigenvalues of the FEM methods versus the analytical solutions (assuming the Dirichlet boundarycondition). The solid ball has 10 voxels in the radius direction. The cylinder has 10 voxels in both the radius direction andheight. The cube and the cuboid have 10 × ×
10 and 10 × ×
20 voxels, respectively. A thicker line is used for the cubicspectrum to avoid overlapping. λs <
1. Hence, determination of the eigenvaluesin a relatively upper part of the spectrum (when sorted by magnitude) requires very large meshes with smalland not wildly variable in size elements. Fortunately, the SPC spectral methods we develop utilize the lower,or leading, part of the spectrum only.
With respect to the computational and storage costs for creating and storing the A and B matrices, theshape functions and their integrals can be calculated beforehand, as discussed at the end of section 3.3. Forthe case of surface meshes, we only need to calculate the metric tensor g for each triangle, which has thecomputational complexity of O ( T ), with T being the number of triangles in a mesh. To construct the A and B matrices given the metric tensors, the computational complexity is linear in the number of non-zeroelements in A and B , which is (cid:0) (cid:1) per triangle for linear FEM and (cid:0) (cid:1) per triangle for cubic FEM. Overall,the construction of the FEM Laplacians has order O ( T ) for both linear and cubic methods, with a largercoefficient for the cubic case. The constructions of the A and B matrices in the voxel case is much simplerthan in the surface case, since the metric tensor g remains constant and does not need to be recalculatedfor each finite element. As a result, the construction of the A and B matrices only depends on the numberof non-zero elements and has the computational complexity of O ( v ), with v being the number of voxelsdenoted as “active”. Again the cubic FEM is expected to have a larger coefficient for the computationalcomplexity compared to the linear FEM due to the increased number of nodes per voxel as the matrices areless sparse. The computational cost of solving for the first eigenvalues is mentioned in section 3.5.1.The storage cost is similar to the computational cost of constructing the A and B matrices given themetric tensor g , as it depends on the number of non-zero elements of each matrix. However, the storagecost is expected to be smaller because two adjacent triangles can share the same pair of points that counts13 igure 9: The first 50 eigenvalues of the FEM Laplacians converge to the true analytical LB spectrum (assuming the Dirichletboundary condition) as we have more voxels in the radius direction, resulting in a more refined voxel approximation of theperfect solid ball with radius 1. as one non-zero element but calculated twice in the A and B matrices, once for each triangle they appearin. Take the linear FEM for a 2D closed mesh as an example, let d i be the degree of point i , which isthe number of edges connected with point i in a mesh, then row i in matrix A or B has ( d i + 1) non-zeroelements. Furthermore, since both A and B matrices are symmetric, only the lower or upper triangularpart needs to be stored, so the smallest storage needed for matrix A or B is (cid:80) i ( d i / O ( N ), with N being the mesh size. We want to point out that thelinear FEM Laplacians are the most sparse that a discretized Laplacian could be, because only interactionsbetween directly connected points are taken into account. The exact storage cost for cubic FEMs and 3Dvoxel cases is more complicated, but it is easy to see it will again be linear in the number of nodes, wherethe number of additional nodes for the cubic FEM depends on the number of edges (in both 2D and 3Dcases) and the number of triangles (in the 2D case only). Another advantage of the FEM methods over the methods used in ZD to estimate the LB spectrum isthat the boundary conditions in the Helmholtz equation can be conveniently implemented, and this permitsthe control and inspection of parts from partial or open meshes that can easily result due to “occlusion”(unreachable areas to the scanner), in contrast to the spectral SPC method in Zhao & Del Castillo (2020)which can only be applied to closed objects with no boundaries. As discussed in Appendix A, the Dirichlet( f ≡
0) and Neumann ( ∂f∂n ≡
0) boundary conditions simplify the weak form from (24) to (4). Further,for the Dirichlet condition, we can just omit the boundary points and their corresponding shape functions,which are constantly zero over the whole surface and do not contribute to the A and B matrices at all.For the Neumann condition, Reuter (2006) suggests simply treating the boundary points as inner points.Either boundary case allows for the statistical process monitoring of “open” meshes, as opposed to “closed”meshes using the FEM spectra. Figure 10 plots the leading spectra of different LB approximations as wellas the true LB spectra. The sequence of noise-free unit spheres have larger holes on their surface from leftto right in the figure, and all estimated spectra are affected. The Dirichlet LB spectrum is more sensitiveto the presence of the holes compared to the Neumann LB spectrum which changes less in the presence ofthe increasing holes, staying closer to the true (analytic) spectrum of the whole (hole-free) surface. Underthe same boundary condition, there is not much difference between the linear FEM and cubic FEM, sincethey both converge fast as shown in Figure 7. Both boundary conditions can be implemented for the voxelFEM methods, discussed above, in the same way. 14 igure 10: LB spectrum estimates for an open surface mesh with a hole. The graphs show the leading spectrum of the FEMLaplacians versus that of the localized Laplacian (Li et al., 2015) for a unit sphere with differently sized boundaries (holes).The FEM Laplacian spectra are closer to the true LB spectrum of a complete unit sphere than that of the localized Laplacian(Li et al., 2015), especially under Neumann boundary conditions. In contrast, the Dirichlet FEM LB spectrum is more sensitiveto the presence of the hole. Run length behavior
A standard performance metric of any SPC chart is the out of control run length , defined as the numberof parts sampled between a defect (or out of control condition) occurs in a sequence of measured parts andwhen this is detected by the chart mechanism (Montgomery, 2020). Also important is the in-control runlength, defined as the number of parts sampled between false detections, when the process is actually in astate of control. One seeks short out of control run lengths and long in-control run lengths. Closed formexpressions for the in and out of control run length distributions of most SPC charts are intractable, andit is customary to estimate the Average Run Length (ARL) and the standard deviation of the run length(SDRL) using Monte Carlo simulation. In this section we adopt the nonparametric control chart in Chenet al. (2016), called the distribution-free multivariate exponentially-weighted moving average (“DFEWMA”)chart, and apply it to monitor changes in the first 15 eigenvalues of the estimated FEM LB spectra of partsvia simulation. It is crucial to use a nonparametric chart, given that the non-normality (non-gaussianity)of the LB spectra of measured parts has been observed even in simulated surface meshes with isotropic,uncorrelated normal-distributed noise (Zhao & Del Castillo, 2020). In Appendix B we provide a briefoverview of the DFEWMA chart operation and its tuning parameters, which will be referred to in thissection. To gain a more complete sense of the effectivity of the SPC chart using the estimated FEM spectra,we conduct the run length analysis based on simulated objects in different scenarios and compare the runlength performance against some other previously proposed methods for SPC of 3D objects. We considerdifferent practical cases, from varying the noise structure to generalizing the type of data, from 2D meshesto 3D voxels, in both cases including cases where the mesh is open.
The simplest case is where the coordinate measurements contain noise, due to the combination of measure-ment error and manufacturing error, which is uncorrelated and isotropic in space. We show two examplesunder these conditions, a prototype part with three types of local defects as the out-of-control scenarios, anda series of cylindrical parts with a “barrel-like” shape controlled by δ >
The first example is a prototype part used in Zhao & Del Castillo (2020), which is typical in an additivemanufacturing process. Three types of defects are considered, namely two types of “chipped” corner partsand a part with a “protrusion” in one of the top edges, shown earlier in Figure 1. To simulate manufacturingand measuring noise, isotropic N (0 . I ) noise is added to the coordinate of each point and to simulatethe case of unregistered meshes with unequal number of (non-corresponding) points, between 0 and 5 pointsare randomly deleted from each simulated part, resulting in mesh sizes of 1675-1680 points. Table 1 showsthe average run lengths (ARL) and the standard deviation of the run lengths (SDRL) for the FEM methodscompared to the Li et al. (2015) Laplacian and a registration based method proposed by Zhao & Del Castillo(2020) as a benchmark. The in-control case uses a nominal ARL of 20 to avoid long simulation times, andit is achieved by all methods thanks to the DFEWMA control chart, which is easy to tune for a desiredin-control ARL (Chen et al., 2016). Changing the DFEWMA chart design parameters such that the nominalin-control ARL is 200, most of the methods are able to signal within the first 5 defective parts when theprocess is out of control. Furthermore, both FEM methods are able to detect the defects quickly withminimal requirements on mesh sizes and mesh qualities, while the Li et al. (2015) Laplacian needs a largermesh with a higher mesh quality, obtained from the Loop subdivision algorithm (see Loop (1987) and Zhao& Del Castillo (2020)) as a pre-processing step, to better capture the local defects. This indicates that theFEM methods are far more sensitive than the Li et al. (2015) Laplacian to reflect local shape changes.We also simulated cases with “open” meshes, as they show one of the greatest advantages of the FEMLB method. In this case, all the bottom of the mesh and a portion of the interior of the cylindrical regionof the prototype part are deliberately omitted to better represent the practical case when a range scanner16n-control Part Chipped Table 1: Run length performance of the DFEWMA SPC charts applied to the FEM LB spectrum for the test part exampleswith closed meshes. Results are obtained from 10,000 replications. Chart parameters were set at m = 100 , w min = 1 , w max =10 , λ = 0 . , and α = 0 .
005 for the out-of-control cases and α = 0 .
05 for the in-control case (resulting in shorter in-control runlengths to reduce simulation expense). First 15 LB operator eigenvalues were used. All methods achieve the nominal in-controlrun length, but there are notable differences in our of control performance, with the FEM methods and the use of the IteratedClosest Point (ICP) method (proposed in Zhao & Del Castillo (2020) as a benchmark) performing the best. The computationalexpense and non-convexity of the ICP formulation makes this method impractical. Due to the nature of the DFEWMA chart,a run length of 2.0 is the minimum possible that can be achieved.Figure 11: Prototype parts with open meshes used for the run length analysis in Table 2. The left most figure is the CADmodel, while the other three figures show the different types of defects (chipped corners or a protrusion). All meshes are openand have boundaries, where the bottom of the part and a portion of the interior of the cylindrical region are missing, a situationcommon due to regions that are unreachable to a scanner. encounters unreachable regions of an object during a scan. The resulting meshes are shown in Figure 11 andthe corresponding estimated run length results using the FEM spectra are listed in Table 2. The Dirichletboundary condition is applied because it is more sensitive to the “holes” in the mesh, as shown in Figure10, and also results in smaller sparse matrices after deleting the rows and columns corresponding to theboundary points. From the table, both the linear FEM and cubic FEM are able to detect the out-of-controlscenarios immediately, while achieving the nominal in-control run length behavior at the same time. Bycomparing the results with Table 1, we can see open meshes have slightly longer out-of-control run lengthsthan closed meshes, but the difference is negligible. This case proved the applicability of the FEM methodsfor SPC of open meshes, a realistic case when a range sensor is mounted in a fixed position and cannot “see”the object from all perspectives. In-control Part Chipped
Table 2: Run length performance of the DFEWMA charts applied to the FEM LB spectrum for the test part examples with openmeshes. Results are obtained from 10,000 replications. Chart parameters were set at m = 100 , w min = 1 , w max = 10 , λ = 0 . , and α = 0 .
005 for the out-of-control cases and α = 0 .
05 for the in-control case. First 15 LB operator eigenvalues were used. .1.2 “Barrel-like” cylindrical parts To parameterize the out-of-control run length in a simple way, we consider cylindrical parts acquiringa more “Barrel-like” shape as an OC parameter δ > δ times the standard deviation of the noise to the radius, so the deformed radius at height h becomes 10 + 0 . δ sin( hπ/ N (0 , . I ) is added to the coordinates of the points. Table 3 comparesthe ARL and SDRL of the different methods, including the various LB spectra and the ICP based method,as a function of the OC parameter δ . The FEM spectra detect small changes into the “barrel” shape ofthe cylinder (small δ ) faster than the Li et al. (2015) spectrum. The ICP method is less efficient because itonly considers the average deviation from the CAD model, and small differences can be easily masked bythe overall natural variability of the in-control process.Linear FEM Cubic FEM Li et al. (2015) ICP objectiveARL SDRL ARL SDRL ARL SDRL ARL SDRL δ = 0 20.07 19.57 20.09 19.60 20.46 20.24 20.25 19.88 δ = 0 . δ = 0 .
005 2.03 0.16 2.00 0.00 2.03 0.19 39.76 65.89 δ = 0 . δ = 1 2.00 0.00 2.00 0.00 2.00 0.00 5.19 3.01 δ = 2 2.00 0.00 2.00 0.00 2.00 0.00 2.03 0.17 δ = 3 2.00 0.00 2.00 0.00 2.00 0.00 2.00 0.00 δ = 10 2.00 0.00 2.00 0.00 2.00 0.00 2.00 0.00 Table 3: Phase II out-of-control run length performance of the DFEWMA charts applied to the different LB spectra and ICPobjective for barrel-shaped cylindrical parts. 10,000 replications, each with 100 IC cylinders followed by a sequence of defectivecylinders until detection. The DFEWMA chart parameters are: m = 100 , w min = 1 , w max = 10, and λ = 0 . α = 0 .
005 for δ > α = 0 .
05 for δ = 0, corresponding to an in-control ARL of 200 and 20, respectively. First 15 LB operator eigenvalueswere used, and mesh sizes varied between 1995 and 2005 non-corresponding points. A more general case in practice is when the noise is correlated and non-isotropic, created, e.g., by man-ufacturing noise that changes spatially on the surface of the objects depending on how the cutting tooloperates. To analyze this case, we repeated the analysis for the barrel-shape cylinders to show the runlength performance of the FEM methods under correlated spatial noise. The defective shape is introducedin the same way as before, so the deformed radius at height h is still 10 + 0 . δsin ( hπ/ .
05 are the nominal radius, nominal height, and the standard deviation of noise, respectively. At eachpoint p i = (cid:0) p i,x , p i,y , p i,z (cid:1) , non-isotropic and spatially correlated noise e i = (cid:0) e i,x , e i,y , e i,z (cid:1) is added to thepoint coordinate. The covariance functions between different noise terms are: Cov ( e i,k , e j,l ) = σ e −| p i,k − p j,l | /r k if i (cid:54) = j, k = lσ + σ if i = j, k = l k (cid:54) = l Here i, j are point indices and k, l ∈ { x, y, z } indicate the axes. To keep a constant total level of noise,we fixed σ + σ = 0 . . Similarly to the previous example, the mesh sizes for the cylindrical partsrandomly vary between 1995 and 2005 points to model non-corresponding, different size meshes, and thefirst 15 eigenvalues are used. Table 4 shows the results. As it can be seen, again the FEM spectra havea similar and slightly better detection compared to the Li et al. (2015) spectrum when σ = 0 .
02, or inother words, σ (cid:54) = 0. When σ = 0 .
05, the FEM spectra are affected by the spatial correlation and fail to18etect small changes quickly. This is likely due to the irregular triangulations that result from the correlateddata, as irregular meshes cause stability problems for FEM. Applying a mesh pre-processing method priorto computation of the LB spectra notably improves the run length performance. σ = 0 . , r x = r y = 2 . σ = 0 . , r x = r y = 5 . σ = 0 . , r x = r y = 2 . δ Table 4: Phase II out-of-control run length performance of the DFEWMA charts applied to the different LB spectra and ICPobjective with barrel-shaped cylindrical parts under spatially correlated, non-isotropic noise. 10,000 replications, each consistingof 100 IC cylinders followed by defective cylinders until detection. DFEWMA charts were applied to both LB spectrum methodand ICP with parameters: m = 100 , w min = 1 , w max = 10 , λ = 0 .
01 and α = 0 . r z = 16 . We also compare the Phase II run length behavior of our FEM LB spectrum methods with an existing SPCmethod for 3D geometrical data due to Colosimo et al. (2014), which is based on Gaussian Processes. Itshould be pointed out that this is a method aimed at contact sensed data and hence assumes small, equallysized meshes with corresponding points from part to part distributed in a lattice pattern, and is a methodthat performs GPA registration of the points first. Their method cannot handle the harder problem ofnon-contact data, where the numbers of points per part varies and points do not correspond from part topart, and would have trouble if points did not form a lattice. Still, Table 5 shows how that the spectralFEM SPC method is very competitive in these unfavorable circumstances, and even sometimes providesbetter run length performance. The FEM methods once again show better run length performance in thiscase than the localized Laplacian spectral method used by Zhao & Del Castillo (2020).LB Spectra GP sub unif GP sub lh Linear FEM Cubic FEM Li et al. (2015)In Control ARL 100.57 96.63 99.85 99.69 100.77(SDRL) (100.71) (91.74) (94.33) (97.41) (100.94)Quadrilobe ARL 4.17 3.60 6.29 4.70 1.39 δ = 0 . δ = 0 . Table 5: Out-of-control run length performance comparisons of different LB spectra versus two Gaussian processes (GPs)models studied in Colosimo et al. (2014). Results obtained from 1,000 replications. DFEWMA chart parameters are: m =100 , w min = 1 , w max = 10 , λ = 0 .
01 and α = 0 .
01, corresponding to an in-control ARL of 100. First 15 LB operator eigenvalueswere used, and the mesh size is fixed to 1054 points. Mesh pre-processing based on the Loop method is applied for all threeLB spectrum methods, resulting in around 2000 points. Results of the GP methods were originally reported in Table 3, inColosimo et al. (2014). igure 12: Voxel objects used for the SPC run length analysis in Table 6. The first row shows the noise-free models of fourparts with an inner cylinder of different eccentricity and orientation. The second row of figures shows sample realizations ofthe corresponding parts when noise is added. The CAD model is the second part from the left, with a perfect circular cylinder(eccentricity=0) that corresponds to a radius value Rx = 8. The other three parts have elliptical cross-sections of differenteccentricities of 0.4581, 0.4841, and 0.6614, corresponding to Rx parameter values of 9, 7 and 6, respectively. Finally we consider a run length analysis for 3D voxel data. To the best of our knowledge, there has not beena Statistical Process Control scheme proposed for voxel sensor data in the literature, so only the performanceof the linear and cubic FEM spectra are compared. To simulate the volumetric datasets obtained via CTscans of a part with inner features, we consider a cube with a hollow cylinder inside it, see Figure 12. Thedimension of the cube is 20 × ×
10 voxels, and the nominal radius of the hollow cylinder is 8 voxels,shown in the second column. For the out-of-control scenarios, we vary the radius of the cylinder alongone of the axis, denoted by Rx , to make increasingly more elliptic cylinders. The in-control part has acylindrical hole with a cross section with eccentricity= 0 (parametrized with a value of Rx = 8). We chose Rx = 9, Rx = 7, Rx = 6 as three types of defectives, corresponding to ellipses of different orientations andeccentricities, namely, 0.4581, 0.4841, and 0.6614, respectively (recall circles have eccentricity = 0 and forparabolas eccentricity = 1, with ellipses in the range 0 < eccentricity < xyz coordinates and an additional value givingthe opacity or intensity of the material. In our case, we will assume a rendering algorithm is applied suchthat each voxel is classified simply as “active” or “inactive” based on its intensity compared against a giventhreshold value. If the intensity is higher, then it indicates that the voxel contains enough material in itand should be denoted as “active”. Therefore, noise is more likely to occur in the voxels near the boundaryof the object being scanned due to the natural variability of the intensity measurement. For this reason,we added noise in two different ways, making “active” boundary voxels that in reality should be “inactive”,and vice versa. The level of noise is parametrized by the maximum number of voxels allowed to have their“active/inactive” statuses switched, specified in column “Max Noise” of Table 6. For example, when “MaxNoise” is 25, then for each noisy realization, an integer between 1 and 25 is randomly selected with equalprobability and that particular number of “active” boundary voxels are randomly selected to be “inactive”.Next, this procedure is repeated to select “inactive” boundary voxels to be “active”. Simulated parts withnoise are shown in the bottom row of plots in Figure 12.The run length results of applying the DFEWMA nonparametric chart for 3D solids with different levelsof noise and types of defects are summarized in Table 6. The Dirichlet boundary condition is applied20ecause it gives smaller Laplacian matrices, which is beneficial especially in the voxel case since the numberof vertices increases dramatically as we have more voxels. Overall, the FEM spectra have excellent detectionpower, and the cubic FEM is consistently outperforming the linear FEM because it is more accurate asstated in Section 3.5.2. As expected, increasing noise level makes it harder and results in slightly longerARLs and SDRLs. Among the 3 types of defects, Rx = 9 is the hardest to detect, because it has the smallesteccentricity and hence is the closest to the CAD model. On the other hand, Rx = 6 with eccentricity 0.6614can be easily detected by the 3D FEM spectra regardless of the noise level.Max Noise Rx Linear FEM Cubic FEMARL SDRL ARL SDRL25 9 2.00 0.02 2.00 0.007 2.00 0.00 2.00 0.006 2.00 0.00 2.00 0.0050 9 2.53 0.69 2.23 0.447 2.00 2.00 2.00 0.006 2.00 0.00 2.00 0.00100 9 4.67 2.52 3.95 1.897 2.03 0.16 2.00 0.036 2.00 0.00 2.00 0.00
Table 6: Phase II out-of-control run length performance of the DFEWMA charts applied to the linear and cubic FEM LB spectrawith voxel data. The Dirichlet boundary condition is applied. 10,000 replications, each consisting of 100 IC parts followed bydefective parts until detection. DFEWMA charts were applied with parameters: m = 100 , w min = 1 , w max = 10 , λ = 0 . α = 0 . We have presented a new approach for the Statistical Process Control of 3-dimensional parts whose metrologyis acquired with either range sensors (surface data) or CT scanners (volumetric data). The new approachis based on the computation of the Laplace-Beltrami operator spectrum, an operator that codifies thegeometrical properties of an object around a point. The LB spectrum, being an intrinsic geometricalproperty of an object, permits the comparison of different parts without the need of “part localization”(registration) algorithms, which are hard nonconvex optimization problems whose ad-hoc solutions mayresult in increased noise (Zhao & Del Castillo, 2020). In contrast to the spectral method recently proposedby Zhao & Del Castillo (2020), the new method, based on estimating the LB spectrum from solving aHelmholtz boundary equation via FEM, is more accurate, can be estimated to do SPC on both mesh andvoxel data, and can be applied to the important practical case of “open” meshes with holes, thanks to theexplicit incorporation of the boundary conditions in the PDE. We have shown how the SPC run lengthperformance of a nonparametric control chart that uses the FEM LB spectrum as a “profile” feature frompart to part, is very competitive and in most cases, better, than that of existing state of the art SPCmethods for 3D data (without requiring registration/part localization). The method was also demonstratedon meshes with boundaries, typical when the scanner is mounted at a fixed position and does not havereach to all the part, and on 3D volumetric data, returned by CT scans. The cubic FEM method providesconsistently better run length performance than the linear FEM, at a higher storage and computationalcost. Fortunately, the linear FEM has almost as good run length performance as the cubic FEM approach.As can be seen from our discussion on computational and storage requirements (section 3.5.2), boththe computational and storage cost mainly depend on the number of nodal points, which directly affectsthe dimension of the A and B matrices. This can potentially be a limitation for the cubic FEM method,especially for the 3D volumetric case, where two additional nodes are inserted for each edge. The secondlimitation for the FEM methods is revealed in section 4.2, where they are less sensitive to detect part defectsunder stronger correlated noise, although the effect of this problem is reduced when a pre-processing of the21esh is applied prior to our methods. Still, this reflects the dependency of FEM methods on the propertiesof the mesh, well-known in the field of PDEs, which, in our case, are reflected in the ability of a SPC chartmechanism to detect an out of control state.One potential limitation of the spectral FEM methods in general, not directly related to our SPCproposal, is that they return two sparse matrices ( A and B ) instead of one discretized Laplacian matrix L ,which, though does not affect the eigenvalues and hence our methods, may impose difficulties if one wishesto derive other geometrical properties from the estimated Laplacians. For example, equation (2) shows howthe mean curvature at a point can be estimated by simply multiplying the LB operator times the pointcoordinates (since ∆ M p ( u, v ) = 2 H n ( u, v )). However, this cannot be done for the FEM Laplacians sincematrix B is not invertible and hence L = B − A is not possible to compute. Supplementary materials
MATLAB code that implements the FEM LB spectrum estimation methods for both mesh and voxel casesis provided, including the models used in this paper.
References
Arnold, D. N., & Awanou, G. (2011). The serendipity family of finite elements.
Foundations of ComputationalMathematics , , 337–344.Babu, M., Franciosa, P., & Ceglarek, D. (2017). Adaptive measurement and modelling methodology forin-line 3d surface metrology scanners. Procedia CIRP , , 26–31.Belkin, M., Sun, J., & Wang, Y. (2008). Discrete laplace operator on meshed surfaces. In Proceedings of thetwenty-fourth annual symposium on Computational geometry (pp. 278–287).Borg, I., & Groenen, P. J. (2005).
Modern multidimensional scaling: Theory and applications . SpringerScience & Business Media.del Castillo, E., & Zhao, X. (2020). Statistical process monitoring for manifold data. In
Wiley Stat-sRef: Statistics Reference Online (pp. 1–8). John Wiley & Sons. doi: https://doi.org/10.1002/9781118445112.stat08276 .Chen, N., Zi, X., & Zou, C. (2016). A distribution-free multivariate control chart.
Technometrics , ,448–459.Colosimo, B. M., Cicorella, P., Pacella, M., & Blaco, M. (2014). From profile to surface monitoring: Spc forcylindrical surfaces via gaussian processes. Journal of Quality Technology , , 95–113.Helffer, B., & Sundqvist, M. (2016). On nodal domains in euclidean balls. Proceedings of the AmericanMathematical Society , , 4777–4791.Huang, D., Du, S., Li, G., Zhao, C., & Deng, Y. (2018). Detection and monitoring of defects on three-dimensional curved surfaces based on high-density point cloud data. Precision Engineering , , 79–95.Ihlenburg, F. (2006). Finite element analysis of acoustic scattering, . .Kac, M. (1966). Can one hear the shape of a drum? The American Mathematical Monthly ,
73, part II ,1–23.Kreyszig, E. (1991).
Differential Geometry . Differential Geometry. Dover Publications. URL: https://books.google.com/books?id=P73DrhE9F0QC .Kruth, J. P., Bartscher, M., Carmignato, S., Schmitt, R., De Chiffre, L., & Weckenmann, A. (2011).Computed tomography for dimensional metrology.
CIRP annals , , 821–842.22e Dret, H., & Lucquin, B. (2016). Partial differential equations: modeling, analysis and numerical approx-imation volume 168. Springer.Lee, J. M. (2018).
Introduction to Riemannian manifolds . Springer.Li, X., Xu, G., & Zhang, Y. J. (2015). Localized discrete laplace–beltrami operator over triangular mesh.
Computer Aided Geometric Design , , 67–82.Loop, C. (1987). Smooth subdivision surfaces based on triangles. Master’s thesis, University of Utah,Department of Mathematics , .Marsden, J. E., & Tromba, A. (2012).
Vector calculus . (6th ed.). W.H. Freeman & Co.Montgomery, D. C. (2020).
Introduction to Statistical Quality Control . John Wiley & Sons.Niethammer, M., Reuter, M., Wolter, F.-E., Bouix, S., Peinecke, N., Koo, M.-S., & Shenton, M. E. (2007).Global medical shape analysis using the laplace-beltrami spectrum. In
International Conference on MedicalImage Computing and Computer-Assisted Intervention (pp. 850–857). Springer.Patan´e, G. (2016). Star-laplacian spectral kernels and distances for geometry processing and shape analysis.In
Computer Graphics Forum (pp. 599–624). Wiley Online Library volume 35.Reuter, M. (2006).
Laplace spectra for shape recognition . Books on Demand.Reuter, M., Biasotti, S., Giorgi, D., Patan`e, G., & Spagnuolo, M. (2009). Discrete laplace–beltrami operatorsfor shape analysis and segmentation.
Computers & Graphics , , 381–390.Reuter, M., Niethammer, M., Wolter, F.-E., Bouix, S., & Shenton, M. (2007). Global medical shape analysisusing the volumetric laplace spectrum. In (pp.417–426). IEEE.Reuter, M., Wolter, F.-E., & Peinecke, N. (2006). Laplace–beltrami spectra as ‘shape-dna’ of surfaces andsolids. Computer-Aided Design , , 342–366.Strang, G., & Fix, G. (1973). An Analysis of the Finite Element Method . Prentice Hall, Englewood Cliffs,NJ.Wardetzky, M., Mathur, S., K¨alberer, F., & Grinspun, E. (2007). Discrete laplace operators: no free lunch.In
Symposium on Geometry processing (pp. 33–37). Aire-la-Ville, Switzerland.Wells, L. J., Megahed, F. M., Niziolek, C. B., Camelio, J. A., & Woodall, W. H. (2013). Statistical processmonitoring approach for high-density point clouds.
Journal of Intelligent Manufacturing , , 1267–1279.Zang, Y., & Qiu, P. (2018). Phase ii monitoring of free-form surfaces: An application to 3d printing. Journalof Quality Technology , , 379–390.Zhao, X., & Del Castillo, E. (2020). An intrinsic geometrical approach for statistical process control ofsurface and manifold data. Technometrics (to appear) , (pp. 1–39).
Appendix A . Some details on the variational solution of Helmholtzequation
To get the weak or variational form of Helmholtz equation (3), we first multiply it times a test function φ ∈ C φ ∆ f = λφf (23)23ntegrating both sides of the equation over the surface area and applying Green’s first identity (a directconsequence of Gauss’ divergence theorem, see Marsden & Tromba (2012), p. 475) yields: (cid:90) M φ ∆ f dV = (cid:90) Γ φ ( ∇ f · n ) ds − (cid:90) M ∇ φ · ∇ f dV = (cid:90) M λφf dV (24)where dV is either the surface element on a 2-dimensional manifold (surface) or a volume element in a3-manifold M , Γ is the boundary of M , ds is either a length or area element on the boundary Γ and ∇ f · n = ∂f∂ n . Both the Dirichlet ( f, φ ≡
0) and Neumann ( ∂f∂n ≡
0) boundary conditions satisfy (cid:90) Γ φ ( ∇ f · n ) ds = 0which simplifies (24) to the so-called Garlekin, weak, or variational form of the PDE, equation (4).The underlying variational problem solved by the weak form equation (4) is not that frequently discussedin applied FEM references. Strang & Fix (1973) indicate how the equation of the weak form corresponds tothe Euler-Lagrange equation of a minimum energy functional, which, for our Helmholtz problem (∆ − λ ) u ≡L u = 0 is of the form: I ( u ) = (cid:104)L u, u (cid:105) = (cid:90) M L u u dV A perturbation around the minimizer f of I ( u ) is next introduced of the form f + (cid:15)φ where (cid:15) is arbitrarilysmall. It can then be shown that the variation of the functional I ( u ) is minimized when (cid:104)L f, φ (cid:105) = 0 , butnote how this expression yields precisely the weak formulation, equating (24) to zero. Appendix B. Operation and parameters of the DFEWMA SPCchart
All the run length analysis in the paper were obtained with the Chen et al. (2016) distribution-free multi-variate exponentially-weighted moving average (“DFEWMA”) chart. This is a chart that operates in theon-line or “Phase II” stage typical of SPC methods. We refer readers to Zhao & Del Castillo (2020) formethods for “Phase I”, or the parameter learning stage with in-control data. The Phase I method shownthere can be directly applied with the more accurate and flexible FEM LB spectra presented in this paper,and was not discussed in the present paper for conciseness.The DFEWMA chart is a distribution-free multivariate control chart. Suppose m acceptable parts areavailable from Phase I, and we want to test the n th manufactured part in Phase II assuming all previous n − p estimated LB eigenvalues of each part i be X i ∈ R p ,then all existing parts can be represented by X , · · · , X m , X m + , · · · , X m + n . The chart calculates thefollowing statistic T jn ( w, λ ) = (cid:80) ni = n − w +1 (1 − λ ) n − i R jni − E (cid:2)(cid:80) ni = n − w +1 (1 − λ ) n − i R jni (cid:3)(cid:113) Var (cid:0)(cid:80) ni = n − w +1 (1 − λ ) n − i R jni (cid:1) where R jni is the rank of the j th eigenvalue from the i th part, X i,j , among the j th eigenvalues of all partsranging from X ,j to X m + n,j . Here w is the window size and λ is the weight in this EWMA-type of chart.In our run length analyses, we chose λ = 0 .
01, and w = min { max { n, w min } , w max } with w min = 1, w max = 10for quicker detection of small changes. The use of a window implies that the quickest possible detection isgreater than one part or sample. Both the expectation and variance terms in T jn ( w, λ ) can be analyticallyderived, see more details in Zhao & Del Castillo (2020). The DFEWMA chart monitors the sum of squares T n ( w, λ ) = (cid:80) pj =1 T jn ( w, λ ) given that differences in all the first pp