An Intrinsic Geometrical Approach for Statistical Process Control of Surface and Manifold Data

TL;DR

This paper introduces an intrinsic geometrical method using Laplace-Beltrami spectrum for real-time statistical process control of manufactured parts, effectively detecting small shape and size changes without registration issues.

Contribution

The paper proposes a novel intrinsic geometrical approach based on LB spectrum for SPC, improving sensitivity and computational efficiency over registration-based methods.

Findings

More sensitive to small shape and size changes

Effective with large, unequal mesh datasets

Applicable to various 3D data types

Abstract

We present a new method for statistical process control (SPC) of a discrete part manufacturing system based on intrinsic geometrical properties of the parts, estimated from three-dimensional sensor data. An intrinsic method has the computational advantage of avoiding the difficult part registration problem, necessary in previous SPC approaches of three-dimensional geometrical data, but inadequate if noncontact sensors are used. The approach estimates the spectrum of the Laplace-Beltrami (LB) operator of the scanned parts and uses a multivariate nonparametric control chart for online process control. Our proposal brings SPC closer to computer vision and computer graphics methods aimed to detect large differences in shape (but not in size). However, the SPC problem differs in that small changes in either shape or size of the parts need to be detected, keeping a controllable false alarm rate and without completely filtering noise. An online or "Phase II" method and a scheme for starting up in the absence of prior data ("Phase I") are presented. Comparison with earlier approaches that require registration shows the LB spectrum method to be more sensitive to rapidly detect small changes in shape and size, including the practical case when the sequence of part datasets is in the form of large, unequal size meshes. A post-alarm diagnostic method to investigate the location of defects on the surface of a part is also presented. While we focus in this article on surface (triangulation) data, the methods can also be applied to point cloud and voxel metrology data.

Tables

Table 1. Table 1: In-control run length performance of the DFEWMA charts applied to the LB spectrum and ICP objective. Results are obtained from 10,000 replications. Chart parameters were set at m 0 = 100 , w min = 1 , w max = 10 , λ = 0.01 , formulae-sequence subscript 𝑚 0 100 formulae-sequence subscript 𝑤 min 1 formulae-sequence subscript 𝑤 max 10 𝜆 0.01 m_{0}=100,w_{\text{\tiny min}}=1,w_{\text{\tiny max}}=10,\lambda=0.01, and α = 0.05 𝛼 0.05 \alpha=0.05 which corresponds to a geometric in-control ARL of 20. The cylinder and the part in Figure 8 are both equal to their CAD model plus isotropic noise N ( 0 , 0.05 2 ​ 𝐈 3 ) 0 superscript 0.05 2 subscript 𝐈 3 (0,0.05^{2}{\bf I}_{3}) . First 15 LB operator eigenvalues were used.

	LB spectrum		ICP objective
	ARL	SDRL	ARL	SDRL
Geometric Distribution	20.0000	19.4936	20.0000	19.4936
In-control cylinder	20.4628	20.2390	20.2492	19.8825
In-control part	20.1654	19.8863	20.1335	19.2041

Table 2. Table 2: Phase II out-of-control run length performance of the DFEWMA charts applied to the LB spectrum and ICP objective for barrel-shaped cylindrical parts, 10,000 replications, each with 100 IC cylinders followed by a sequence of defective cylinders until detection. Chart parameters are: m 0 = 100 , w min = 1 , w max = 10 , λ = 0.01 formulae-sequence subscript 𝑚 0 100 formulae-sequence subscript 𝑤 min 1 formulae-sequence subscript 𝑤 max 10 𝜆 0.01 m_{0}=100,w_{\text{\scriptsize min}}=1,w_{\text{\scriptsize max}}=10,\lambda=0.01 and α = 0.005 𝛼 0.005 \alpha=0.005 , corresponding to an in-control ARL of 200. First 15 LB operator eigenvalues were used, and mesh sizes varied between 1995 and 2005 points.

	LB spectrum		ICP objective
	ARL	SDRL	ARL	SDRL
$\delta =0.0005$	10.7949	9.9443	83.2122	122.6474
$\delta =0.005$	2.0336	0.1851	39.7576	65.8904
$\delta =0.5$	2.0000	0.0000	31.4895	51.8724
$\delta =1$	2.0000	0.0000	5.1867	3.0066
$\delta =2$	2.0000	0.0000	2.0262	0.1695
$\delta =3$	2.0000	0.0000	2.0000	0.0000
$\delta =10$	2.0000	0.0000	2.0000	0.0000

Table 3. Table 3: Phase II out-of-control run length performance of the DFEWMA chart applied to the LB spectrum and the ICP objective for the prototype part, isotropic uncorrelated noise. Results from 10,000 replications, each with 100 IC parts followed by a sequence of defective parts until detection. DFEWMA chart parameters were: m 0 = 100 , w min = 1 , w max = 10 , λ = 0.01 , α = 0.005 formulae-sequence subscript 𝑚 0 100 formulae-sequence subscript 𝑤 min 1 formulae-sequence subscript 𝑤 max 10 formulae-sequence 𝜆 0.01 𝛼 0.005 m_{0}=100,w_{\text{\scriptsize min}}=1,w_{\text{\scriptsize max}}=10,\lambda=0.01,\alpha=0.005 (in-control ARL = 200). The first 15 eigenvalues of the LB operator were used. Mesh sizes were necessarily small (1675 to 1680 points) with pre-processing based on the Loop method resulting in close to 5000 points.

Defect type	LB spectrum		LB spectrum		ICP objective
	(original mesh)		(pre-processed mesh)
	ARL	SDRL	ARL	SDRL	ARL	SDRL
Chipped #1	158.1168	182.1803	5.0943	2.7707	2.0000	0.0000
Chipped # 2	91.3750	135.3055	4.4362	2.1155	2.0000	0.0000
Protrusion	3.6450	1.7183	2.4341	0.5115	2.0000	0.0000

Table 4. Table 4: Phase II out-of-control run length performance of the DFEWMA charts applied to the LB spectrum and ICP objective with barrel-shaped cylindrical parts under spatially correlated, non-isotropic noise. 10,000 replications, each consisting of 100 IC cylinders followed by defective cylinders until detection. DFEWMA charts were applied to both LB spectrum method and ICP with parameters: m 0 = 100 , w min = 1 , w max = 10 , λ = 0.01 formulae-sequence subscript 𝑚 0 100 formulae-sequence subscript 𝑤 min 1 formulae-sequence subscript 𝑤 max 10 𝜆 0.01 m_{0}=100,w_{\text{\scriptsize min}}=1,w_{\text{\scriptsize max}}=10,\lambda=0.01 and α = 0.005 𝛼 0.005 \alpha=0.005 , corresponding to an in-control ARL of 200. First 15 LB operator eigenvalues used, mesh size varied between 1995 and 2005 points. r z = 16.7 subscript 𝑟 𝑧 16.7 r_{z}=16.7 for all three cases.

${\sigma }_1^2$	$r_x,r_y$	$\delta$	LB spectrum		ICP objective
			ARL	SDRL	ARL	SDRL
${0.02}^2$	$2.6$	$0.0005$	11.4954	13.8908	198.0817	196.8560
		$0.005$	2.0315	0.1852	204.6340	200.7492
		$1$	2.0000	0.0000	14.5372	17.6595
${0.05}^2$	$2.6$	$0.0005$	8.4815	6.8473	202.2944	201.3962
		$0.005$	2.0087	0.0950	203.9053	199.6954
		$1$	2.0000	0.0000	26.6083	38.8395
${0.05}^2$	$5.2$	$0.0005$	7.9981	6.0534	198.0303	202.7248
		$0.005$	2.0077	0.0886	194.4391	195.0273
		$1$	2.0000	0.0000	28.6818	46.9065

Table 5. Table 5: Out-of-control run length performance comparisons of LB spectrum versus a Gaussian process (GP) model method in Colosimo et al.  , ( 2014 ) . Results obtained from 1,000 replications. DFEWMA chart parameters are: m 0 = 100 , w min = 1 , w max = 10 , λ = 0.01 formulae-sequence subscript 𝑚 0 100 formulae-sequence subscript 𝑤 min 1 formulae-sequence subscript 𝑤 max 10 𝜆 0.01 m_{0}=100,w_{\text{\scriptsize min}}=1,w_{\text{\scriptsize max}}=10,\lambda=0.01 and α = 0.01 𝛼 0.01 \alpha=0.01 , corresponding to an in-control ARL of 100. First 15 LB operator eigenvalues were used, and the mesh size is fixed to 1054 points. Mesh pre-processing based on the Loop method is applied for the LB spectrum, resulting in around 2000 points. Results of the GP methods were originally reported in Table 3, in Colosimo et al.  , ( 2014 ) .

		LB spectrum	GP${}_{\text{sub\_unif}}$	GP${}_{\text{sub\_lh}}$
In Control	ARL	99.8490	99.6890	100.7650
	(SDRL)	(94.3316)	(97.4138)	(100.9427)
Quadrilobe	ARL	6.2860	4.6980	1.3890
$\delta =0.00185$	(SDRL)	(2.6505)	(3.9684)	(0.7709)
Half frequency	ARL	3.2710	14.1130	4.5090
$\delta =0.00075$	(SDRL)	(1.4372)	(13.4398)	(4.0450)
Tapering	ARL	61.5650	80.8110	83.1830
$\delta =0.1$	(SDRL)	(76.9659)	( 77.6667)	(83.5048)

Table 6. Table 6: In-control alarm probability of the “Phase I” scheme by Capizzi & Masarotto, ( 2017 ) applied to a cylinder and the part design used before, with different number of LB eigenvalues. The nominal false alarm probability is α = 0.05 𝛼 0.05 \alpha=0.05 . Results are obtained from 10,000 replications. All parameters are at their default values as in the R package dfphase1 .

	Eigenvalues used:
	2nd-5th	2nd-15th	2nd-50th	2nd-75th	2nd-100th
In-control cylinder	0.0496	0.0501	0.0527	0.0504	0.0540
In-control prototype part	0.0478	0.0465	0.0505	0.0473	0.0550

Table 7. Table 7: Out-of-control alarm probability of the “Phase I” scheme by Capizzi & Masarotto, ( 2017 ) for the cylindrical parts with increasing out-of-control parameter δ 𝛿 \delta . Different number of LB eigenvalues were investigated. The IC nominal FAP is α = 0.05 𝛼 0.05 \alpha=0.05 . Results are obtained from 10,000 replications. All default parameters in R package dfphase1 were used.

	Eigenvalues used:
	2nd-5th	2nd-15th	2nd-50th	2nd-75th	2nd-100th
$\delta =0.0005$	0.4224	0.4860	0.2360	0.1326	0.0553
$\delta =0.005$	1	1	1	1	0.0741
$\delta =0.5$	1	1	1	1	0.0948
$\delta =1$	1	1	1	1	0.0929
$\delta =10$	1	1	1	1	0.0876

Table 8. Table 8: Out-of-control alarm probabilities of the “Phase I” SPC scheme by Capizzi & Masarotto, ( 2017 ) for the part defects in Figure 8 . Different number were investigated. The nominal in-control FAP is α = 0.05 𝛼 0.05 \alpha=0.05 . Results obtained from 10,000 replications. All default parameters in Rpackage dfphase1 were used.

	Eigenvalues used:
Defect type	2nd-5th	2nd-15th	2nd-50th	2nd-75th	2nd-100th
Chipped #1	0.0595	0.0640	0.0709	0.0975	0.0531
Chipped #2	0.0920	0.0969	0.1534	0.5443	0.0562
Protrusion	0.9001	0.9952	1	1	0.0921

Equations

$$\displaystyle G(\mbox{\boldmath$X$}_{1},\mbox{\boldmath$X$}_{2},..,\mbox{\boldmath$X$}_{N})$$

$$\min_{\mbox{\boldmath$M$},\mbox{\boldmath$\Gamma$},\mbox{\boldmath$\gamma$}}L(\mbox{\boldmath$M$},\mbox{\boldmath$\Gamma$},\mbox{\boldmath$\gamma$})={\min_{\mbox{\boldmath$M$},\mbox{\boldmath$\Gamma$},\mbox{\boldmath$\gamma$}}}\sum_{i=1}^{m_{1}}\sum_{j=1}^{m_{2}}M_{ij}\cdot C(\mbox{\boldmath$\Gamma$}\mbox{\boldmath$x$}_{q,i}+\mbox{\boldmath$\gamma$},\mbox{\boldmath$x$}_{p,j})$$

$${\mathbf{p}}_u=\frac{\partial \mathbf{p}(u,v)}{\partial u}={\left(\frac{\partial x(u,v)}{\partial u},\frac{\partial y(u,v)}{\partial u},\frac{\partial z(u,v)}{\partial u}\right)}^{\prime }$$

$${\mathbf{p}}_v=\frac{\partial \mathbf{p}(u,v)}{\partial v}={\left(\frac{\partial x(u,v)}{\partial v},\frac{\partial y(u,v)}{\partial v},\frac{\partial z(u,v)}{\partial v}\right)}^{\prime }.$$

$$\frac{d\alpha (t)}{dt}=\frac{\partial \mathbf{p}}{\partial u}\frac{du(t)}{dt}+\frac{\partial \mathbf{p}}{\partial v}\frac{dv(t)}{dt}={\mathbf{p}}_u\frac{du(t)}{dt}+{\mathbf{p}}_v\frac{dv(t)}{dt}.$$

$${\text{I}}_p\left(\frac{d\alpha (t)}{dt}\right)$$

$$s(t)={\int }_0^t\sqrt{d{s}^2}dt.$$

$$\nabla f=\left(\frac{\partial f}{\partial x_1},\frac{\partial f}{\partial x_2},....,\frac{\partial f}{\partial x_n}\right)$$

$$div\mathbf{F}=\nabla \cdot F=\frac{\partial F_1}{\partial x_1}+\frac{\partial F_2}{\partial x_2}+....+\frac{\partial F_n}{\partial x_n}$$

$$\Delta f=-div\nabla f=-\sum _{i=1}^n\frac{{\partial }^{2}f}{\partial x_i^2}$$

$${\Delta }_{\mathcal{M}}f=-di{v}_{\mathcal{M}}{\nabla }_{\mathcal{M}}f$$

$${\Delta }_{\mathcal{M}}\mathbf{p}(u,v)=-di{v}_{\mathcal{M}}{\nabla }_{\mathcal{M}}\mathbf{p}(u,v)=2H\mathbf{n}(u,v)\in {\mathbb{R}}^3$$

$${\Delta }_{\mathcal{M}}f=-\frac{1}{\sqrt{det(g)}}\sum _{j=1}^k\frac{\partial }{\partial x^j}\left(\sqrt{det(g)}\sum _{i=1}^k{g}^{ij}\frac{\partial f}{\partial x^i}\right)$$

$${\Delta }_{\mathcal{M}}f=\lambda f$$

$$\underset{i\to \infty }{\lim }\frac{{\lambda }_i}{i}=\frac{4\pi }{\text{mbox}Area(\mathcal{M})}$$

$$\sum _{i=1}^{\infty }{e}^{-{\lambda }_{i}t}\le \frac{\text{mbox}Area(\mathcal{M})}{2\pi t}.$$

$${\lambda }_1\le \frac{8\pi (\mathcal{G}+1)}{\text{mbox}Area(\mathcal{M})}$$

$$L_{\mathcal{K}}^{t}f(p_i)=\frac{1}{4\pi t^2}\sum _{T\in T_{\mathcal{K}}}\frac{A(T)}{3}\sum _{p_j\in V(T)}{e}^{-||p_i-p_j|{|}^2/4t}(f(p_i)-f(p_j)),j=1,2,..,m$$

$$L_{\mathcal{K}}^t=D-W$$

$$\displaystyle{L^{t}_{\mathcal{K}}f(p_{i})=\frac{1}{4\pi t^{2}}\int_{y:d(p_{i},p_{j})\leq r}e^{-\frac{d(p_{i},p_{j})^{2}}{4t}}\big{[}f(p_{i})-f(p_{j})\big{]}dp_{j}}$$

$$\displaystyle{=\frac{1}{4\pi t^{2}}\sum_{y:d(p_{i},p_{j})\leq r}\frac{A(y)}{3}e^{-\frac{d(p_{i},p_{j})^{2}}{4t}}\big{[}f(p_{i})-f(p_{j})\big{]}}$$

$$\tilde{L}=B^{-1}L$$

$$t_{\max }=\underset{1\le j\le p}{\max }\frac{|{\sum }_{i=1}^m{X}_i(j)/m-{\sum }_{i=1}^n{X}_{m+i}(j)/n|}{S{E}_j}$$

$$S{E}_j=\frac{\sqrt{(m-1){\sigma }_{1,m,j}^2+(n-1){\sigma }_{m+1,m+n,j}^2}}{\sqrt{\frac{1}{m}+\frac{1}{n}}}$$

$$T_j=\frac{{\sum }_{i=1}^m{R}_{ji}-E[{\sum }_{i=1}^m{R}_{ji}]}{\sqrt{\text{mbox}Var({\sum }_{i=1}^m{R}_{ji})}},j=1,...,p$$

$$T_{jn}(w,\lambda )=\frac{{\sum }_{i=n-w+1}^n(1-\lambda {)}^{n-i}{R}_{jni}-E\left[{\sum }_{i=n-w+1}^n(1-\lambda {)}^{n-i}{R}_{jni}\right]}{\sqrt{Var\left({\sum }_{i=n-w+1}^n(1-\lambda {)}^{n-i}{R}_{jni}\right)}}$$

$$E\left[\sum _{i=n-w+1}^n(1-\lambda {)}^{n-i}{R}_{jni}\right]=\{\begin{array}{ll}w\frac{m_0+n+1}{2}& \lambda =0\\ \frac{1-(1-\lambda {)}^w}{\lambda }\frac{m_0+n+1}{2}& \lambda \ne 0\end{array}$$

$$\text{mbox}Var\left(\sum _{i=n-w+1}^n(1-\lambda {)}^{n-i}{R}_{jni}\right)=\{\begin{array}{ll}\frac{w(m_0+n+1)(m_0+n-w)}{12}& \lambda =0\\ \frac{1-(1-\lambda {)}^{2w}}{2\lambda -{\lambda }^2}\frac{(m_0+n+1)(m_0+n-1)}{12}& \\ -\left(\frac{1-\lambda -(1-\lambda {)}^w}{{\lambda }^2}-\frac{(1-\lambda {)}^2-(1-\lambda {)}^{2w}}{{\lambda }^2(2-\lambda )}\right)\frac{m_0+n+1}{6}& \lambda \ne 0\end{array}$$

$$Cov(e_{i,k},e_{j,l})=\{\begin{array}{ll}{\sigma }_1^2{e}^{-|p_{i,k}-p_{j,l}|/r_k}& \text{ if }i\ne j,k=l\\ {\sigma }_1^2+{\sigma }_2^2& \text{ if }i=j,k=l\\ 0& \text{ if }k\ne l\end{array}$$

$$L_{\mathcal{K}}^{t}f(p_i)$$

$$=\frac{1}{4\pi t^2}\sum _{T\in T_{\mathcal{K}}}\sum _{p_j\in V(T)}\frac{A(T)}{3}{e}^{-||p_i-p_j|{|}^2/4t}(f(p_i)-f(p_j))$$

Full text

Keywords: Laplace-Beltrami operator, Spectral method, Permutation test, Differential Geometry, non-contact sensor.

1 INTRODUCTION

Widespread use of modern sensors in engineering and industry in recent times has resulted not only in bigger datasets but also in more complex datasets. Statistical Process Control (SPC) is an area that has witnessed increased sophistication in metrology accompanied by increased complexity in the resulting data sets related to production or manufacturing processes (Colosimo *et al. *, 2014; Colosimo, 2018). We consider the case non-contact sensors (or a combination of contact and non-contact sensors) collect geometric data formed by thousands of points which actually lie on what could be described as a 2-dimensional curved space or manifold, that is, the surface of the object or part produced. Due to manufacturing and measurement errors, the observed surface will deviate from the target or nominal surface, given by the design of the part, usually available in some Computer Aided Design (CAD) file.

The purpose of this paper is to lay the foundations of a new methodology for SPC of discrete-part 3D geometrical data that can be in either of the form of a point cloud, a mesh, or, even voxel (volumetric) datasets, although we focus in this paper in the mesh case, the most common type of part data generated by non-contact scanners today. Prior approaches to this problem utilize methods where the points scanned in each part need to be registered (or “superimposed” one to one) and require that the exact same number of points in each part be scanned, tenable assumptions only when data are exclusively collected using contact sensors. Our main contributions are to provide a new SPC methodology for 3D geometrical data, based on intrinsic geometrical properties of datasets scanned from a sequence of parts, that: i) does not require registration of the points, meshes or voxel datasets scanned from each part, which is a computational expensive and non-convex problem for which no guarantee of global optimality can be given, and therefore, it is prone to error, and ii) does not require parts to have the exact same number of points.

Our approach brings SPC of discrete-part manufacturing closer to the Computer Graphics/Vision fields. The approach we propose is based on techniques popular in computer graphics to characterize 3D objects, and these methods have also been used in Machine Learning of general manifolds of much higher dimension than the 2-manifolds (surfaces) we focus in this paper. In computer graphics and computer vision applications, however, the problem to solve is the identification of large differences, usually evident to the human eye, between the shapes of objects (frequently neglecting differences in size) for classification purposes or to match a query object. In computer graphics, one works with noise free meshes drawn by an artist or CAD engineer, while in computer vision noisy measurements are obtained but usually the noise is first filtered out. Here, in contrast, we focus on detecting considerable smaller differences in either shape or size (sometimes not perceptible to the human eye) in the presence of measurement and manufacturing noise in a sequence of objects which the manufacturer is trying to produce consistently and with low variability, hence differences will tend to be rather small. Furthermore, in addition to detecting changes in the mean shape or size of an object, a noise increment could also be considered a process signature change to be detected, rather than a signal that must be filtered and ignored as in computer vision, where shape identification methods that are robust with respect to noise are typically sought.

Traditional treatment of 2 and 3-dimensional point cloud datasets in Statistics pertains to the field of Statistical Shape Analysis (SSA, Kendall, 1984; Dryden & Mardia, 2016; for applications in manufacturing see del Castillo & Colosimo, 2011). In SSA, the $m$-point cloud data are represented by a configuration matrix {\mbox{\boldmathX}}\in\mathbb{R}^{m\times n} with $n=2$ or 3. The Shape of an object is defined as the geometrical information in $X$ that remains after discounting the effect of similarity transformations usually excepting reflections (translations, rotations and dilations). To make inferences on the shape of $N$ objects, the Generalized Procrustes Algorithm, GPA, is first applied. The GPA registers or superimposes all the $N$ objects by finding scaling factors $\beta_{i}\in\mathbb{R}$, rotation matrices \mbox{\boldmath\Gamma}_{i}\in SO(n) (the special orthogonal group, which excludes reflections and has determinant one) and $n-$dimensional translation vectors \mbox{\boldmath\gamma}_{i},i=1,...,N, such that they minimize the sum of squared full procrustes distances between all pairs of configuration matrices (d_{F}(\mbox{\boldmathX}_{i},\mbox{\boldmathX}_{j})):

[TABLE]

where \mbox{\boldmath1}_{m} is a vector of $m$ ones. Constraints must be added to avoid the trivial solution where all parameters are zero (Dryden & Mardia,, 2016). Note how two objects with different sizes may still have the same shape, given that the effect of dilations (changes of scale) is usually filtered out in SSA. Neglecting differences in size is an aspect in common to shape classification in computer vision. Other shape analysis methods based on euclidean distances between the points (Lele,, 1993) require large distance matrices and will not be reviewed further here.

The problem of registering different three dimensional objects each with a large but different number of (non-corresponding) points has been known for a long time in the computer vision literature, where the Iterated Closest Point (ICP) algorithm (Besl & McKay, 1992; Zhang, 1992) is a standard. Consider the configurations of two distinct unlabeled objects \mbox{\boldmathX}_{q}\in\mathbb{R}^{m_{1}\times n} and \mbox{\boldmathX}_{p}\in\mathbb{R}^{m_{2}\times n} (with $n=3$), not necessarily having the same pose and assume $m_{1}\leq m_{2}$. Let \mbox{\boldmathM}\in\mathbb{R}^{m_{1}\times m_{2}} with $M_{ij}=1$ if \mbox{\boldmathx}_{q,i}\in\mbox{\boldmathX}_{q} is matched with point \mbox{\boldmathx}_{p,j}\in\mbox{\boldmathX}_{p}, and zero otherwise. The objects may be located and oriented differently in space, and hence the problem is not only to find the correspondences but also a rigid body transformation T(\mbox{\boldmathx})=\mbox{\boldmath\Gamma}\mbox{\boldmathx}+\mbox{\boldmath\gamma} that registers the two objects such that the following problem is solved:

[TABLE]

subject to $\sum_{j=1}^{m_{2}}M_{ij}=1,i=1,...,m_{1}$, and $M_{ij}=0$ or 1, where C(\mbox{\boldmath\Gamma}\mbox{\boldmathx}_{q,i}+\mbox{\boldmath\gamma},\mbox{\boldmathx}_{p,j}) is the cost of matching point \mbox{\boldmathx}_{q,i} to point \mbox{\boldmathx}_{p,j}. This is a hard non-linear discrete optimization problem. Existing heuristics differ by choosing different cost functions $C$. In the ICP method, C(\mbox{\boldmath\Gamma}\mbox{\boldmathx}_{q,i}+\mbox{\boldmath\gamma},\mbox{\boldmathx}_{p,j})=||\mbox{\boldmath\Gamma}\mbox{\boldmathx}_{q,i}+\mbox{\boldmath\gamma}-\mbox{\boldmathx}_{p,j}||, the euclidean distance between \mbox{\boldmathx}_{q,i} in \mbox{\boldmathX}_{q} and its closest point \mbox{\boldmathx}_{p,j} in \mbox{\boldmathX}_{p} after the transformation.

Commercial CAD and inspection software use variants of the ICP method to align the cloud points, meshes or voxel datasets of each scanned object and that of the CAD file, in order to determine regions in the manufactured part that differ from nominal. For instance, to do this in the mesh data case, a CAD model, usually in the form of NURBS curves, is sampled to form a mesh or triangulation, after which the alignment of the CAD and scanned part triangulations can be performed. Figure 1 shows an instance of a metal part CAD model and a color-coded comparison between the CAD model and the manufactured part (this is actually voxel data, but the registration problem is essentially identical).

A first approach we will consider as a benchmark for Statistical Process Control of part surface data is based on monitoring the deviations from nominal shown in Figure 1 after applying the ICP method, similarly to a “DNOM” (deviation from nominal) control chart (see e.g., Farnum, 1994). The deviations from nominal are vectors, and either their norm or their individual components could be used for SPC. The optimal value of the ICP statistic (2), L^{*}(\mbox{\boldmathM},\mbox{\boldmath\Gamma},\mbox{\boldmath\gamma}) between CAD file and each scanned part could be monitored on-line, for instance, to provide a “generic” SPC mechanism against a wide variety of unanticipated out of control states in the geometry of a part, in a similar sense to what Box & Ramírez, (1992) thought of univariate Shewhart charts, with other SPC or diagnostic mechanisms added to detect more specific, or more localized defects on the part. As far as we know, this simple strategy has not been proposed before in the literature, so we contrast it with the new intrinsic differential-geometrical methods that conform our main contribution and with some earlier SPC approaches that also require registration.

In this paper, we follow the traditional SPC paradigm with modern statistical tools. The main goal is detection of significant part to part differences with respect to historical variation because they may indicate a manufacturing process problem, and the aim is to detect “assignable causes of variation” as soon as possible, while avoiding false alarms using a statistical monitoring scheme. We consider both of what is called “Phase 1” and ”Phase 2” SPC. Figure 2 shows a “pipeline” of the methods we propose in this paper, starting from scanned data from the surface of an object, estimation of the Laplace-Beltrami operator of the surface and its spectrum, to SPC methods for Phase I and II, and finally, post-alarm diagnostics.

The paper is organized as follows. Section 2 first reviews some preliminary mathematical concepts in order to present the new differential-geometric SPC methods, including the main concept we will use, the Laplace-Beltrami (LB) operator, and its spectrum. The spectrum of the LB operator needs to be estimated from data, and Section 3 discusses methods to do so. Section 4 presents the main practical results, where a specific distribution-free multivariate chart is used to monitor a process with respect to an in-control reference data set (i.e. “Phase II” in SPC) using the spectrum of the estimated discrete LB operator of a sequence of scanned parts. We show how the spectral methods have greater sensitivity to detect out of control conditions in a discrete-part manufacturing process than either using the deviations from nominal ICP method sketched above or using instead earlier SPC approaches. A post SPC alarm diagnostic is presented in Section 5 where we investigate the use of ICP to localize the occurrence of defects on a part. To make the presentation of our SPC proposal complete, Section 6 discusses a method based on the spectrum of the LB operator for the first phase (“Phase I”) when monitoring a process in the absence of prior in-control data. The paper ends with conclusions and some directions for further research. The Appendices in the Supplementary Materials (on-line) contain proofs and derivations, further discussion about the relation between the LB operator and both the heat equation and the combinatorial Laplacian in networks, and a discussion about the applicability to SPC of other intrinsic, spectral distances popular in Computer Graphics and Manifold Learning that have received recent attention in the literature. Supplementary materials also include Matlab code that implements our methods and the datasets used.

2 PRELIMINARIES

In this section we first define the concepts we use in the sequence, in particular, those leading to the definition of the Laplace Beltrami operator (definition 7 below) which is our main object of study. For more on these definitions and concepts see, e.g,., Kreyszig, (1991) and O’Neill, (2006) and Appendix C (supplementary materials, on-line). In Differential Geometry, one begins with properties affecting the vicinity of a point on a surface and deduces properties governing the overall structure of the surface or manifold under consideration.

The preliminary concepts in this section apply to either surfaces or volumes, from which we assume a scanner takes sample measurements that constitute the datasets to be used. The surfaces or volumes are instances of a $k$-dimensional manifold ($k=2$ or 3, respectively) contained in 3-dimensional Euclidean space. Informally, a manifold $\mathcal{M}$ is a $k$-dimensional space that resembles $\mathbb{R}^{k}$ (Euclidean space) on small domains around each of its points (we will assume $\mathcal{M}$ is compact and connected). The manifold hypothesis, useful if true in machine learning and in engineering data analysis, indicates that high ($n$-) dimensional data frequently lie on or near a lower, often curved $k$-dimensional manifold, where $k<n$. In so-called Riemannian manifolds, smooth manifolds with an inner product, we can measure distances, areas, volumes, etc. If the manifold hypothesis holds in a dataset, we say the intrinsic dimension of the data is $k$, and that the data manifold is embedded in an $n$-dimensional ambient space. For manifold data, any $n$-dimensional point $x$ can be completely described by defining $k$ local (intrinsic) coordinates or “parameters” $x^{1},x^{2},...,x^{k}$. For instance, consider a parametric curve in $\mathbb{R}^{3}$, such as the helicoidal curve \mbox{\boldmathp}(t)=\left(x(t)=\cos(t),y(t)=\sin(t),z(t)=t\right), where $t\in D\subset\mathbb{R}$ is the curve coordinate or parameter. Here the intrinsic dimension of the manifold is $k=1$ (i.e., a 1-manifold) with $x^{1}=t$, while the ambient space dimension is $m=3$. In this paper, we will mostly consider data sampled from surfaces or 2-manifolds of manufactured parts, although all our methods are extendable to the case of 3-manifolds, i.e., voxel (volumetric) data.

Definition 1. Any property of a manifold $\mathcal{M}$ that can be computed without recourse to the ambient space coordinates, and instead is computed using only the intrinsic or local manifold coordinates $x^{1},...,x^{k}$ is said to be an intrinsic geometrical property, or simply, an intrinsic property, of $\mathcal{M}$.

Thus, to describe data points on a surface, such as geographical data on Earth, we need only two coordinates $x^{1},x^{2}$, so the surface of a sphere is intrinsically 2-dimensional.

Definition 2. Any geometrical property of an object that remains constant after application of a given transformation is said to be invariant with respect to that transformation.

Intrinsic geometrical properties of a manifold in Euclidean space are invariant with respect to rigid transformations (rotations and translations, but not dilations), but the opposite is not true. The SPC methods we will present are intrinsic, and therefore invariant, and it is thanks to these properties that the groupwise part registration problem can be avoided. An instance of Definition 2 are Euclidean distances between points in a configuration matrix $X$, which are invariant to rigid transformations but they are evidently not intrinsic. An instance of an intrinsic property is the geodesic distance between two points located on a manifold $\mathcal{M}$, which is therefore also invariant with respect to rigid transformations. In Appendix D (supplementary materials, on-line), we discuss intrinsic distances other than the geodesics and their potential use for SPC.

A classic result in Differential Geometry indicates that intrinsic geometrical properties of a manifold depend only on the so-called first fundamental form of $\mathcal{M}$, a quadratic form defined next for 2-manifolds (surfaces), the case we concentrate in this paper. Consider a parametric surface ${\bf p}(u,v)=(x(u,v),y(u,v),z(u,v))^{\prime}$, $(u,v)\in D\subset\mathbb{R}^{2}$ (here $u=x^{1},v=x^{2})$ which defines $\mathcal{M}$, a Riemannian 2-manifold, i.e., a smooth (so derivatives can be computed) surface. Define the surface differential vectors at ${\bf p}(u,v)$ as:

[TABLE]

and

[TABLE]

Now define a parametric curve on $\mathcal{M}$, $\alpha(t)={\bf p}(u(t),v(t))$ such that ${\bf p}(u_{0},v_{0})=\alpha(0)$ and use the chain rule:

[TABLE]

Finally, take the inner product (borrowed from $\mathbb{R}^{3}$):

[TABLE]

where: $g_{11}={\left\langle{\bf p}_{u},{\bf p}_{u}\right\rangle}\quad g_{12}={\left\langle{\bf p}_{u},{\bf p}_{v}\right\rangle}\quad g_{22}={\left\langle{\bf p}_{v},{\bf p}_{v}\right\rangle}$. We then have the following.

Definition 3. The quadratic form $\textup{I}_{p}(\alpha(t)^{\prime})=(ds)^{2}$ is called the first fundamental form of the parametric surface ${\bf p}(u,v)$ describing $\mathcal{M}$. It provides a means to measure arc lengths, areas and angles on $\mathcal{M}$. It defines a new inner product for vectors on tangent planes $\left\langle,\right\rangle_{\mathcal{M}}$ and therefore, a metric on the surface, with associated matrix (tensor) $g$:

[TABLE]

In this sense, the ambient space induces a metric, the Riemannian metric, on the manifold $\mathcal{M}$. Since $|{\bf w}|=\sqrt{\left\langle{\bf w},{\bf w}\right\rangle_{\mathcal{M}}}$, the length of a curve segment on $\mathcal{M}$ is:

[TABLE]

With the Riemannian metric $g$, we can also compute differential operators acting on a function defined on $\mathcal{M}$, which are very useful for our purposes, i.e., for estimation, and therefore, statistical monitoring, of intrinsic geometrical properties of a 3D object.

Definition 4. The gradient of a function on $\mathbb{R}^{n}$ points in the direction of steepest ascent and equals to:

[TABLE]

The gradient therefore creates a vector field in $\mathbb{R}^{n}$.

Definition 5. The divergence of a vector field \bf F$$=(F_{1},F_{2},...,F_{n})^{\prime} in $\mathbb{R}^{n}$ is given by:

[TABLE]

and measures the “quantity” of the outward flux of $F$ from the infinitesimal neighborhood around each point $p$. This is a scalar-valued function that creates a scalar field.

Definition 6. The Laplace operator of a twice differentiable function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is minus the divergence of its gradient field:

[TABLE]

and measures the difference between $f(x)$ and the average $f(x)$ in a neighborhood around $x$. Given the second derivatives, it is a measure of curvature, and can be alternatively understood as minus the trace of the Hessian of $f(x)$. The minus sign is for consistency with equation (4) below. Note how the domain of the function here is n-dimensional euclidean space. We next extend this definition to general manifolds, obtaining the main differential-geometric operator used in the sequence, the Laplace-Beltrami operator, widely used in Computer Graphics and Machine Learning (Belkin,, 2003; Kimmel,, 2004; Levy,, 2006; Patané,, 2014; Reuter *et al. *,, 2006).

Definition 7. For a function $f:\mathcal{M}\rightarrow\mathbb{R}$, the Laplace-Beltrami (LB) operator (sometimes called the second differential parameter of Beltrami, see Kreyszig, 1991) is defined as:

[TABLE]

where $div_{\mathcal{M}}$ is the divergence taken on $\mathcal{M}$. For a point defined by a parametric surface ${\bf p}(u,v)$ the following relation holds:

[TABLE]

where ${\bf n}(u,v)$ is the normal at the point ${\bf p}(u,v)$ on $\mathcal{M}$ and $H$ is the mean curvature of $\mathcal{M}$ at ${\bf p}$, which equals the average of the maximum and minimum curvatures at ${\bf p}$. This relation provides a geometric interpretation of the action of the LB operator, which can be visualized as creating a vector field of normals on $\mathcal{M}$ such that the “height” of the normal is twice the mean curvature of $\mathcal{M}$ at that point.

The LB operator extends the definition of the Laplacian to functions defined on manifolds, and is an intrinsic measure of local curvature of a function at a point. Intuitively speaking, this curvature of the function needs to consider also the curvature of the manifold itself, which, contrary to the Laplacian case for a function defined on Euclidean space, is not flat. The LB operator, “contains” the local manifold curvature. The intrinsic nature of the LB operator can be seen from defining a local coordinate system (or parametrization) on the manifold (${\bf x}(x^{1},...,x^{k}$), with $k=2$ for surfaces). Then, the LB operator applied to a function $f(x^{1},...,x^{k})\in\mathcal{C}^{2}$ is defined as:

[TABLE]

where $g^{ij}$ are the elements of \mbox{\boldmathg}^{-1}. The LB operator on $f$ is therefore a function of elements in the metric tensor $g$ only, and thus it is intrinsic and invariant with respect to rigid transformations. This is the key property we exploit for SPC: when considering a sequence of part surface datasets, the spectrum (i.e., eigenvalues, see below) of the corresponding estimated LB operators can be compared directly without any need to register the parts, since it does not matter how the parts were oriented or located when measured, the LB operator remains the same. As discussed below, the spectrum of the LB operator contains considerable additional geometric information about the manifold, and is widely used for this reason in both machine learning and computer graphics/computer vision. For concrete examples of the computation of the LB operator (4) for functions $f$ defined on 3D objects, see Appendix C (on-line supplementary materials).

2.1 The Laplace-Beltrami operator spectrum and some of its properties

We propose to work with the spectrum (collection of eigenvalues) of the estimated LB operator of a part dataset. Formally, the Laplacian eigenvalue problem is

[TABLE]

sometimes called the Helmholtz partial differential equation (see Evans, 2010, pg. 323, and Appendix C, on-line supplementary materials), with an infinite number of pairs of eigenfunctions $f$ and eigenvalues $\lambda$. The collection of eigenvalues $\{\lambda_{i}\}_{i=0}^{\infty}$ (in ascending order) is called the spectrum of the LB operator and the eigenfunctions form an orthonormal basis in $L_{2}(\mathcal{M})$, see Chavel, (1984). In the particular case $\mathcal{M}$ is a circle (1-dimensional manifold) the corresponding basis functions consist of the usual Fourier harmonics $sin(k\pi x)$ and $cos(k\pi x)$.

The spectrum of the LB operator is always discrete, non-negative, and contains considerable geometrical and topological information about a manifold that can be used for shape identification. For instance, Weyl’s law (see Chavel, 1984) indicates that, for a surface $\mathcal{M}$:

[TABLE]

(thus the area of $\mathcal{M}$ can be inferred from the asymptotic slope of the spectrum; note $\lambda_{i}$ is proportional to the index $i$). Another result shown in a classic paper by Kac, (1966) is

[TABLE]

Also, the spectrum contains topological information about $\mathcal{M}$. For instance, one result showing dependency of the spectrum on topological information is that for a surface without boundary (Yang & Yau,, 1980),

[TABLE]

where $\mathcal{G}$ is the genus of the surface (number of holes).

Scaling a surface $\mathcal{M}$ by a factor of $s$ changes the eigenvalues by $1/s^{2}$ (see Figure 4). Also, the spectrum changes in a continuous form as the shape of the manifold deforms, making it possible to monitor small shape changes through the spectrum. Furthermore, useful information can be extracted from only the lower part of the spectrum. According to Reuter *et al. *, (2009), the LB operator spectrum of a manifold has more discrimination power than simpler measures like surface area. These authors provided examples of shapes with the same surface area but different spectrum.

3 CHARACTERIZING THE GEOMETRY OF A 3D OBJECT USING THE LB OPERATOR SPECTRUM

The true spectrum of very few manifolds is known. One instance where it is known is the case of a unit sphere (Figure 4). Repeated eigenvalues are the result of perfect symmetries in the geometry of an object and are therefore rare in practice. To characterize the geometry of general 3D objects, such as discrete parts from a manufactured process, we need to use a discrete approximation of the LB operator based on a sample of points, possibly with adjacency information, forming a mesh. Representations based on voxel data are also possible. In this paper we focus on mesh data and leave voxel laplacians for further research.

3.1 Approximating the LB operator and its spectrum

The LB operator is linear, taking functions into functions. In practice, we have no expression for the surface $\mathcal{M}$ of a part, we only have a (large) sample of points (point cloud dataset) typically returned by most 3D scanners with adjacency information (mesh or triangulation dataset). We can discretize the manifold $\mathcal{M}$ based on the data, and functions defined on $\mathcal{M}$ are reduced to vectors whose elements are the function values on the sampled points. Then an approximate, or discrete LB operator is obtained in the form of a matrix acting on vectors, returning a vector. One then works with the eigenvalues of the matrix which approximate the true spectrum of the LB operator on the original manifold. There are different ways to discretize the manifold where the data lie. In computer graphics and machine learning, it is common to work with triangulations, sometimes generated automatically by non-contact sensors. If the data are in the form of a point cloud, a graph can also be constructed by connecting the nearest neighbors to each point. Finite-element methods (FEM) approximations of the LB operator for mesh data (Reuter *et al. *,, 2006) and for voxel data (Reuter *et al. *,, 2007) have also been developed, which increases the practical use of this differential geometry tool.

Motivation for some of the most popular LB operator approximations used for analyzing the shape of an object comes from the theory of heat diffusion and wave propagation in Physics (see Appendix C, supplementary materials, on-line). For instance, if the sensor data available have the form of a triangulation $\mathcal{K}$, the Mesh Laplacian approximation (Belkin *et al. *,, 2008) is given by:

[TABLE]

where $T_{\mathcal{K}}$ denotes the set of all triangles in the mesh, $A(T)$ denotes the area of triangle $T$, and $V(T)$ denotes the set of vertices in triangle $T$. The parameter $t$ comes from the heat equation, whose relationship to the LB operator is explained in Appendix C (supplementary materials, on-line). Larger $t$ values imply the approximate Laplacian is considering larger areas of interest around a given point.

We will use a recent modification of the Mesh Laplacian (7). We first show how the Mesh Laplacian has a simple expression that connects it to the underlying graph Laplacian, as the next simple result, which is new as far as we know, indicates. This Lemma also applies to the Localized Mesh Laplacian we use. The result is relevant in practice given the interpretability and applications of the graph Laplacian eigenvectors (Appendix C, supplementary materials, on-line).

Lemma 1. For a Riemannian manifold $\mathcal{M}$ sampled at $m$ points, the discretized approximation of its LB operator (7) results in an $m\times m$ Laplacian matrix acting on vectors f$$\in\mathbb{R}^{m} which can be written as:

[TABLE]

with $W_{ij}=\frac{1}{12\pi t^{2}}\sum_{T:p_{j}\in V(T)}A(T)e^{-||p_{i}-p_{j}||^{2}/4t}$ for $i=1,...,m$ and $j=1,...,m$ and diagonal matrix $D$ with entries $D_{jj}=\sum_{i}W_{ij}$.

For a proof, see Appendix A (supplementary materials, on-line). The estimated spectrum of the LB operator is then given by the eigenvalues of matrix $L_{\mathcal{K}}^{t}$. Li *et al. *, (2015) noted how the Mesh Laplacian approximation (7) is a global approximation, since computing it at a given point requires the integral over the whole surface $\mathcal{M}$. Instead, these authors proposed a modification of the Mesh Laplacian that uses geodesic distances between points $p_{i}$ and $p_{j}$ (as opposed to euclidean distances) and that considers in the last sum only points within a certain radius $r$ of each point $p_{i}$, resulting in the alternative discrete Laplacian approximation:

[TABLE]

where $A(y)$ denotes the area of the one ring neighborhood of point $y$. By writing the equation above in vector form, it is easy to obtain the Laplacian matrix in the form of equation (8), $L^{t}_{\mathcal{K}}=D-W$, with $W_{ij}=\frac{1}{12\pi t^{2}}A(p_{j})e^{-d(p_{i},p_{j})^{2}/4t}$ for $i=1,...,m$ and $j=1,...,m$ and diagonal matrix $D$ with entries $D_{ii}=\sum_{j}W_{ij}$. We call this LB approximation the Localized Mesh Laplacian. It has the merit of resulting in sparser $L$ matrices, reducing storage and computational requirements, and was used in all the examples and figures shown below.

**Convergence. **Both Belkin *et al. *, (2008, 2009) and Li *et al. *, (2015) show how as the triangulation $\mathcal{K}$ gets finer, their Laplacians $L_{\mathcal{K}}^{t}$ (7 and 9) converge pointwise to the continuous LB operator $\Delta_{\mathcal{M}}$ defined on a smooth manifold. Dey *et al. *, (2010) further proved the convergence and stability of the spectra of the Mesh Laplacian in Belkin *et al. *, (2008). The following is a Corollary of these results, important for our purposes as we use the spectrum of (9), and its proof is shown in Appendix A (on-line supplementary materials).

Lemma 2. For a smooth manifold $\mathcal{M}\subset\mathbb{R}^{3}$ with spectrum $\{\lambda_{i}\}$ and for the spectrum of the Localized Mesh laplacian (9), ${\{}\lambda^{L}_{i}{\}}$, we have that for fixed $i$, $|\lambda_{i}-\lambda^{L}_{i}|\rightarrow 0$ as $\epsilon\rightarrow 0$, $t\rightarrow 0$, and $\epsilon/t^{4}\rightarrow 0$.

The parameter $\epsilon$ is a function of the density of the mesh (see Appendix A, supplementary materials, on-line), going to zero with denser meshes. Figure 3 illustrates how as the mesh gets denser, it approximates the underlying surface $\mathcal{M}$ better, and this implies the convergence of the spectrum of the Localized Mesh Laplacian (9).

Real eigenvalues of discrete Laplacian approximations. Patané, (2017) pointed out that many of the discretized LB operators proposed in the literature can be represented in a unified way as the multiplication of two matrices

[TABLE]

where $B$ and $L$ are symmetric, positive semi-definite ($B$ is positive definite) matrices and called the mass matrix and the stiffness matrix, respectively. This means $\tilde{L}$ is symmetrizable (Liu *et al. *,, 2012) and guaranteed to have real eigenvalues, as it can be easily seen that $(\lambda,\phi)$ is an eigenvalue-eigenvector pair of (10) if and only if $(\lambda,B^{1/2}\phi)$ is an eigenvalue-eigenvector pair of the symmetric matrix $B^{-1/2}LB^{-1/2}$. We note that although a symmetrizable LB approximation has therefore the desirable property of having real spectrum, not all other desirable properties of their continuous LB operator counterparts can be achieved with a discrete approximation, so there is a “no free lunch” situation, explaining why there are so many discrete approximations proposed to the LB operator, see Wardetzky *et al. *, (2007). For instance, the eigenfunctions of the continuous LB operator form an orthogonal basis in $L_{2}(\mathcal{M})$, however, the eigenvectors of the discrete LB operator approximations $\tilde{L}=B^{-1}L$ do not form an orthogonal basis if using the euclidean inner product, a result of the underlying meshes used for their computation not being uniform on $\mathcal{M}$ (Rustamov,, 2007). The non-uniform mesh in turn results in a non-symmetric discrete Laplacian. The eigenvectors of a discrete approximation $\tilde{L}=B^{-1}L$ do form an orthogonal basis but with respect to the $B$-inner product $\langle\cdot,\cdot\rangle_{B}$, i.e., $\langle\phi_{i},\phi_{j}\rangle_{B}=\phi_{i}^{\prime}B\phi_{j}=0$ for $i\neq j$. For this reason, Rustamov, (2007) indicates that only when the mesh is uniform one can expect a discrete Laplacian to be “faithful” to the continuous LB operator. This indicates that mesh pre-processing techniques may be valuable, and we explore this possibility in some of the examples in Section 4.3 below. The next result, which is new as far as we know, assures both the Mesh and Localized Mesh Laplacians have a real spectrum (for a proof, see Appendix A, supplementary materials, on-line):

Proposition 1. The Mesh Laplacian (7) and the Localized Mesh Laplacian (9) can be written as (10) and therefore their eigenvalues are all real.

In what follows, we use the Localized Mesh Laplacian due to its sparseness advantage over (7). To demonstrate numerically the behavior of this discrete Laplacian, consider Figure 4, which shows its first 10 eigenvalues for a mesh with 3000 points on the unit sphere. Both the true spectrum of the sphere and the spectrum of the discrete approximation are invariant with respect to rigid transformations.

The spectra of widely different objects are considerably different, as Figure 5 indicates. As expected, if the noise is very high and dominates, the spectrum cannot be estimated well (Figure 6). For moderate levels of noise (which includes both manufacturing errors and measurement errors), we demonstrate in Section 4 how the spectrum of the LB operator can still be used for process monitoring purposes.

3.2 Computational complexity and stability of the discrete LB operator spectrum

Our motivation to use intrinsic geometrical methods in SPC is to avoid the registration problem, which is a combinatorial hard, non-convex computational problem. In our method, the computational cost of obtaining the eigenvalues of $L_{\mathcal{K}}^{t}$ is an $O(m)$ to $O(m^{3})$ operation depending on the sparseness of the $m\times m$ matrix.

However, extra speed can be gained from some other special properties of the LB operator besides the sparsity of its estimator. First, as it will be shown when evaluating the run-length performance of our SPC method, only the lower part of the spectrum is needed, so there is no need to compute all the eigenvalues. Second, as discussed after Lemma 2, symmetrizable LB approximations share the same spectrum with some symmetric matrices, which one can work with instead of the original approximation to take advantage of the computational benefits of symmetry, exploited by the popular Arnoldi algorithm, used in Matlab’s function eigs. The Arnoldi algorithm finds the first $k$ eigenvalue-eigenvector pairs of a sparse symmetric matrix, and has a typical complexity of $O(mk)$. Switching between the two eigenvector systems is also easy since the mass matrix $B$ is diagonal as discussed in the proof of Proposition 1. We show below how $k=15$ provides good run length performance for our LB-based SPC method, and hence, the applicability of our methods to very large datasets is possible with a desktop computer.

It is relevant to point out that there are also methods to compute spectral quantities such as heat kernels and diffusion distances (discussed in Appendix D, supplementary materials, on-line) that do not require to compute the LB spectrum first (Patané,, 2014).

A numerical question of importance is if the computation of the spectrum of the LB operator is stable or not. Patané, (2017) gives an analysis of the stability of this computation for each single eigenvalue $\lambda_{i}$, showing how the computation of eigenvalues with multiplicity one is stable numerically. Furthermore, our numerical experiments (e.g. Figure 6) show that the lower spectrum is stable with respect to moderate levels of noise, typically encountered in manufacturing. The computation of the spectra is not stable when there are higher multiplicities (Patané,, 2017), but as mentioned before, repeated eigenvalues occur due to exact symmetries which will be rare from real scanned objects, and our numerical methods did not show evidence of this kind of potential problem.

4 USING THE ESTIMATED LB SPECTRUM AS A TOOL FOR SPC

Our proposal is to use the lower part of the Localized Mesh Laplacian spectrum (9), i.e., the ordered spectrum cropped up to certain maximum index, obtained from scanned parts, and consider each spectrum a profile from which we derive a general statistical process monitoring technique supplemented by additional post-alarm tools to aid the localization of the defects on a part.

4.1 Permutation tests based on the LB spectrum

Empirical evidence presented by Sun *et al. *, (2009b) indicates that commercial 3D scanner noise is not Gaussian. But even for normally distributed, isotropic errors added to a surface, the estimated LB spectra are not multivariate normal. To illustrate, Figure 7 shows the Shapiro-Wilks marginal tests of normality for each of the first 500 eigenvalues of 200 “acceptable” parts simulated based on the prototype part depicted in Figure 8 below to which we added isotropic N$(0,\sigma^{2}\bf{I}_{3})$ measurement noise. It is therefore necessary to develop distribution-free process monitoring methods for the LB spectrum.

To obtain an initial assessment of the detection capabilities of the spectrum of the LB operator, we use 2-sample permutation tests to compare the mean LB spectra between two groups of parts. Represent the spectra of the estimated LB operator (sorted from smallest to largest, up to some given index $p$) of each part $i$ by ${\bf X}_{i}\in\mathbb{R}^{p}$ (not to be confused with configuration matrices $X$ used in section 1), and let the two samples be \{{\bf X}_{1},...,{\bf X}_{m}\}\sim F({\mbox{\boldmath\mu}}_{0}), when the process is in control or acceptable, and \{{\bf X}_{m+1},...,{\bf X}_{m+n}\}\sim F({\mbox{\boldmath\mu}}_{1}) for parts that have an out of control condition. Given the non-normality of the spectrum data, in our preliminary tests we use a nonparametric, distribution-free permutation test for H_{0}:{\mbox{\boldmath\mu}}_{0}={\mbox{\boldmath\mu}}_{1} using two different types of statistics: the maximum t-statistic (Reuter *et al. *,, 2007) and the component-wise Wilcoxon rank-sum statistics (Chen *et al. *,, 2016). The maximum t-statistic is defined as

[TABLE]

Here $X_{i}(j)$ is the $j$th element of $X_{i}$ (in our case, the $j$th eigenvalue of part $i$), and $SE_{j}$ is the pooled standard error estimate of the $j$th eigenvalue

[TABLE]

where $\sigma_{i,k,j}$ is the standard deviation of the $j$th element in sample $\{{\bf X}_{i},{\bf X}_{i+1},...,{\bf X}_{k}\}$.

For specifying the component-wise Wilcoxon rank-sum statistic, denote the rank of the jth eigenvalue from the ith part, with respect of the pool of $m+n$ parts, as $R_{ji}$ ($i=1,...,m+n$ and $j=1,...,p$) and define:

[TABLE]

where $\sum_{i=1}^{m}E[R_{ij}]=m(m+n+1)/2$ and Var$(\sum_{i=1}^{m}R_{ij})=mn(m+n+1)/12$ (see Appendix B, on-line supplementary materials, for a proof). The test statistic is formed by combining the $T_{j}$’s using $T_{0}=\sum_{j=1}^{p}T_{j}^{2}$, with a large value of the test statistic leading to rejection of $H_{0}$.

Both statistics $t_{\max}$ and $T_{0}$ are intuitive, as there is no explicit mapping available from the eigenvalues of the LB operator to the manifold $\mathcal{M}$, and differences in all the first $p$ eigenvalues should be considered jointly as evidence of a difference between the two objects.

Figure 8 shows the results of permutation tests between two samples of a prototypical part of realistic size, with the number of points per part varying between 26750 and 26850 points (it is important to point out that real datasets from a non-contact sensor will not have identical number of points/part). The prototype part is typical of an additive manufacturing process, and the 3 types of defects we consider below, two types of “chipped” corner parts and a part with a “protrusion” in one of the top edges are also typically encountered in an AM process. The first sample (with sample size 10) is a group of acceptable parts, while the second sample (with sample size 5) is a group of parts of the type shown on top of each column, where the first three parts have different types of defects (chipped corners or a protrusion) while the last one is acceptable in the sense of being equal the CAD design of the part plus isotropic white noise. In this exercise, we added isotropic $N(0,0.05^{2}{\bf{I}_{3}})$ noise to all points of all parts, defective and non-defective. The second row of plots in the figure are results from the permutation test using the maximum t-statistic as in (11), while the last row is using the Wilcoxon rank-sum statistic $T_{0}=\sum_{j}T_{j}^{2}$ with $T_{j}$ as in (12), both using $p=500$. In each of the small figures showing the test results, the blue bars are the empirical pdf’s of the test statistic from all permutations, while the red line indicates the observed value of the test statistic, $t_{\max}$ or $T_{0}$ . Both the maximum-t and rank-sum tests are one-sided. The red numbers under the red lines are the estimated p-values for the corresponding tests, defined as the number of permutations with more extreme (larger) values than the observed test statistic, divided by the total number of permutations (3003 for sample sizes 10 and 5, note this is an exact permutation test).

We repeated the experiments in Figure 8 using eigenvalues normalized by surface area, given that from (6) they relate inversely to the surface area of $\mathcal{M}$. For small, localized defects on the surface of a part, the detection capabilities of the LB spectrum using either normalized or unnormalized are very similar, and therefore the normalized spectrum results are not shown here. Since we wish to detect changes in shape and in size, and not only in shape, we suggest using the unnormalized spectrum, which will be used in the following sections.

Figure 9 shows the distributions of the estimated p-values (as defined above) when the permutation test procedure is repeated 1000 times. In these figures, the last column of plots shows the case when both groups of parts consist of acceptable parts, i.e., the null hypothesis is true. As it is easy to show, in such case the theoretical p-value should follow a standard uniform distribution, and this is approximately the case in the depicted histograms of p-values. In the other cases, when we are testing defective parts against acceptable parts, it is desired to have p-values very close to zero, as it is indeed the case.

These comparisons indicate the potential for using the LB spectrum and permutation tests for shape difference detection. These are not however, an SPC scheme, since for on-line, “Phase II” monitoring we require a sequential test, as further discussed next, and a way to initialize the scheme, or “Phase I” SPC, which we detail in Section 6 .

4.2 On-line SPC scheme (“Phase II”)

Given the non-normality of the LB spectrum, for on-line statistical control we recommend to use a multivariate permutation-based control chart. We have used, with some modifications as discussed below, Chen *et al. *, (2016) distribution-free multivariate exponentially-weighted moving average (“DFEWMA”) chart, and applied it to the lower part of the estimated LB spectra. Following their notation, suppose we have $m_{0}$ parts from Phase I and $n$ parts from Phase II, labeled as $\{-m_{0}+1,-m_{0}+2,...,0,1,...,n\}$. Consider the set of the $j$th component of the vectors one wishes to monitor (eigenvalues of the LB spectra, in our case), taken from observation $k$ to the most recent observation $n$, $\mathcal{X}_{k,j}^{n}=\{X_{jk},...,X_{jn}\}$. We wish to test the equality of the location of the samples $\mathcal{X}_{-m_{0}+1,j}^{n-w}$ and $\mathcal{X}_{n-w+1,j}^{n}$, that is, the in-control observations compared to the most recent observations in a “window” of $w$ observations.

The idea of the DFEWMA chart is to compute the “exponentially weighted” rank statistic $(1-\lambda)^{n-i}R_{jni}$ of the last $w$ observations among all IC observations thus far, where $R_{jni}$ is the rank of $X_{ji}$ among $\mathcal{X}_{-m_{0},j}^{n}$. If these ranks are extreme (large or small) this is evidence the process has changed from its IC state. The exponential weights give more weight to the more recent observations within the last $w$ and can be useful to detect smaller process changes. We therefore use the statistic:

[TABLE]

For given $m_{0}$ in-control observations and a false alarm probability $\alpha$ determining the geometrically distributed in-control run lengths (hence the nominal IC average run length is $1/\alpha$), the remaining chart design parameters are thus the weight $\lambda$ and the window size $w$. In Appendix B (supplementary materials, on-line) we derive the following moment expressions, which consider the covariances of the weighted ranks in the sum:

[TABLE]

and

[TABLE]

We note these expressions are corrected from those in Chen *et al. *, (2016) (see Appendix B, on-line supplementary materials, for derivations). A sequential test statistic based on $T_{jn}(w,\lambda)$ statistics is the sum of squares $T_{n}(w,\lambda)=\sum_{j=1}^{p}T_{jn}^{2}(w,\lambda)$, used by Chen *et al. *, (2016) (maximum statistics could also be used instead). Here we report results based on $T_{n}(w,\lambda)$ and the moment expressions above.

To illustrate the use of the resulting permutation chart, we simulated again parts from the CAD part model shown in Figure 8. Mesh sizes varied randomly (in the range 26750 to 26850 vertices). Simulations were comprised of a “Phase I” of 50 in-control parts (nominal plus noise), followed by an on-line “Phase II” where defectives (parts with a protrusion in one of its “teeth”, see Figures 8-9, third column) were introduced starting at part no. 21. Figure 10 shows simulations of the modified DFEWMA charts under different parameters $\lambda$ and $\alpha$. The chart has a variable control limit that depends on the observed realization of the underlying process, thanks to which it has a geometric in-control run length distribution (the main results in Chen *et al. *, (2016) hold as well under the different test statistic (13) we use). In every case, detection of the out of control (OC) state occurred within the first 3 observations.

**DFEWMA chart parameter selection.-**In common to properties of other EWMA-type of charts, smaller weight $\lambda$ leads to quicker detection of small changes, as confirmed in our simulations. Therefore, $\lambda=0.01$ was used in the following numerical results. As for the window size $w$, we use $w=\min\{n,w_{\text{\scriptsize max}}\}$ with $w_{\text{\scriptsize max}}=10$ as the largest window size (as opposed to 20, used in Chen *et al. *, 2016) where $n\in\mathbb{N}$ is the number of currently available Phase II parts.

Chen *et al. *, (2016) recommend that when the number of Phase II parts is smaller than 5, some Phase I parts should be included in the window to keep a smallest window size of 5 and prevent the number of all possible permutations from being too small, resulting in an inaccurate empirical distribution and therefore an inaccurate critical value. However, they also mentioned that bounding $w$ this way may reduce the detection power when a location shift occurs at the beginning of Phase II, which was found to be the case in our run-length simulations. When the smallest window size is set to 5, the DFEWMA chart needs at least 3 out-of-control parts to signal, which means the majority of the parts in the window needs to be out-of-control even if the variables are able to evidently reflect the change right away. To avoid such “masking effect” and truly present the detection ability of our method, while keeping the nice properties of the DFEWMA, we decided not to include Phase I parts in the window and allow the smallest window size, $w_{\text{\scriptsize min}}$, to be 1. To keep an adequate empirical distribution, we set the number of Phase I parts, $m_{0}$, to 100, so that for the most extreme case, when there is only 1 part in Phase II, the number of possible permutations is ${101\choose 1}$, providing enough permutations to form the empirical distribution. From our simulation results, this modification works well and reduces the smallest possible ARL from 3 to 2.

4.3 Run length behavior

To gain a more complete sense of the effectivity of the SPC chart, we conducted a run length analysis based on simulation of cylindrical parts of increasingly more deformed shape, with parts acquiring a more “barrel-like” shape as an OC parameter $\delta>0$ is increased, to permit computation of out of control run lengths parametrized in a simple way (see Figure 11). This is one of the typical out of control signals Colosimo *et al. *, (2014) discussed in the fabrication of cylindrical parts in a lathe process. We also conducted a run length analysis simulating realizations of the prototype part (and its defect types) shown in Figure 8. Given that a run length analysis implies computation of thousands of LB spectra, to avoid long simulation times, we used smaller mesh sizes, with $1995\sim 2005$ points for the cylindrical parts and $1675\sim 1680$ points for the parts in Figure 8. We also permuted the already simulated parts instead of simulating new parts for new replications to further reduce the computational cost of the simulation while keeping the variability of the run lengths in our analysis.

We applied the DFEWMA charts –with the corrected moments as described in Section 4.2– to the top $p=15$ eigenvalues of the LB spectrum and also applied it to the optimal function value returned by the ICP algorithm (2), which is a measure of the difference between two configuration matrices after discounting similarity transformations. Table 1 shows the average run lengths (ARL) and the standard deviation of the run lengths (SDRL) in the in-control processes, where we use both a perfect cylinder with radius 10 and height 50 and the acceptable part in the last column of Figure 8, both with isotropic N$(0,0.05^{2}\bf{I}_{3})$ noise added. As the table shows, all run lengths have nearly geometrical behaviors, which is consistent with the theorem given in Chen *et al. *, (2016).

For the out-of-control analysis of the cylindrical parts, we added a first harmonic with amplitude $\delta$ times the standard deviation of the noise to the radius, so the deformed radius at height $h$ becomes $10+0.05\delta\sin(h\pi/50)$, adding also isotropic $N(0,0.05^{2}{\bf I}_{3})$ noise to the points. Table 2 compares the ARL and SDRL of both methods (LB spectrum and ICP) as a function of the OC parameter $\delta$. The LB spectrum is very sensitive to changes in the shape of an object, taking only two parts to signal the deformation of the cylinder for $\delta\geq 0.005$. On the other hand, the ICP objective is less efficient than the LB spectrum, especially when $\delta$ is very small (note also the large SDRL values). This is because the small increase in the ICP statistic caused by a slight local deformation can be masked by the overall natural variability of the in-control process, plus the inherent variability in the ICP algorithm itself (recall that (2) is a hard non-convex combinatorial problem), making it difficult for the chart to distinguish the change to an OC condition until more parts are available. However, this would not be a problem for the lower part of the LB spectrum because it is reflecting the overall shape of the parts. Chen *et al. *, (2016) recommend their chart when the ratio $m_{0}/p$ is small, although in this case we observed good run length performance when $m_{0}/p=100/15$.

We also conducted the OC run length analysis for the prototype part and defects shown in Figure 8, where we consider three types of defects corresponding to the first three columns in the figure. The ARLs and SDRLs for both methods (LB spectrum and ICP) for the three defective parts are summarized in Table 3. In this case, the ICP method works consistently better than the LB spectrum to detect all three types of defective parts. This is because the three defects are very localized and evident to the eye, making the increase in the ICP objective function quite significant. As these three defect types are local and the mesh size used is very small, they do not change the overall shape of the part enough for the LB spectrum to quickly detect the changes in the process, particularly with the chipped part #1, which is the smallest and more localized type of defect, with a faster detection for the protrusion defect part, which is the largest of the 3 defects relative to the mesh. To demonstrate that with larger meshes the LB spectrum would detect these types of defects quicker, we applied a pre-processing step to the meshes, utilizing the Loop subdivision method (Loop,, 1987), resulting in meshes of around 5000 points instead, and this notably improves the run length performance, event though no new information is added by the pre-processing apart from the original meshes.

For still larger meshes, the performance of the LB-spectrum will further improve. Consider Figure 8 where the average mesh size was 26,800 points, and the p-values of the Wilcoxon rank-sum tests (bottom row) are significant for all three cases. This indicates that the DFEWMA chart statistic will signal faster with larger mesh sizes. We also note that chipped defect #1 has a slightly larger p-value in Figure 8, and in Figure 9 the distribution of its p-values has a thicker and longer right tale than for the chipped #2 and protrusion defects. This shows the chipped #1 defect is harder to detected by nature and explains why it results in the longest run lengths.

Effect of spatially correlated, non-isotropic noise.- The run-length analysis for the barrel-shape cylinders was repeated under spatially-correlated, non-isotropic noise, as manufacturing noise may be spatially correlated on the surface of the objects depending on how the cutting tool operates on the surface. For this case, the defects are the same as before, so the deformed radius at height $h$ is still $10+0.05\delta sin(h\pi/50)$, where $10$, $50$, $0.05$ are the nominal radius, nominal height, and the standard deviation of noise, respectively. At each point $\boldsymbol{p}_{i}=\begin{pmatrix}p_{i,x},p_{i,y},p_{i,z}\end{pmatrix}$, non-isotropic and spatially correlated noise $\boldsymbol{e}_{i}=\begin{pmatrix}e_{i,x},e_{i,y},e_{i,z}\end{pmatrix}$ is added to the point coordinate. The covariance functions between different noise terms are:

[TABLE]

Here $i,j$ are point indices and $k,l\in\{x,y,z\}$ indicate the axes. To keep the same level of noise, $\sigma_{1}^{2}+\sigma_{2}^{2}=0.05^{2}$. Same as in the previous example, the mesh sizes for the cylindrical parts randomly vary between 1995 and 2005 points and the first 15 eigenvalues are used. Table 4 shows the results. As it can be seen, the effect of the spatial correlation on the run length properties is negligible in the LB-spectrum method for the same levels of noise, but badly affects the chart based on the ICP objective.

In summary, the ICP objective method is an effective method to detect evident local defects in small meshes with non-correlated noise, but as the size of the meshes grows or non-isotropic spatially correlated noise increases the registration it requires deems the method either infeasible or ineffective. On the other hand, the LB spectrum remains very sensitive to global shape changes regardless of the covariance structure of noise. For small local changes, the LB spectrum needs denser meshes to signal quickly, a condition automatically satisfied when data are obtained from non-contact sensors.

Comparisons versus other SPC methods for 3D data. We finally compare the Phase II behavior of our LB spectrum method with an existing SPC method for 3D geometrical data due to Colosimo *et al. *, (2014) which is based on Gaussian Processes. It should be pointed out that this is a method aimed at contact sensed data and hence assumes small, equally sized meshes with corresponding points from part to part distributed in a lattice pattern, and is a method that performs GPA registration of the points first. Their method cannot handle the harder problem of non-contact data, where the numbers of points per part varies and points do not correspond from part to part, and would have trouble if points did not form a lattice. Still, Table (5) shows how that our method is very competitive in these unfavorable circumstances, and, even in some cases it actually provides better run length performance.

5 POST ALARM DIAGNOSTICS

We have proposed a multivariate permutation SPC chart on the lower spectrum of the LB operator as a tool to detect general out of control (OC) states in the process that are not precisely defined, similarly to the role standard Shewhart charts have, as described by Box & Ramírez, (1992). Once an alarm is triggered, an investigation of the specific assignable cause that is normally carried out should include the localization of the defect on the part surface or manifold $\mathcal{M}$, a task we now describe.

Suppose we have a part that has triggered an alarm in the SPC charts described above and we have also have available a CAD model for the part being produced. In order to localize the defect on each part, we apply the ICP algorithm to register the CAD model and the part that triggered the alarm (we use the ICP algorithm implemented by Bergström & Edlund, 2014). Upon completion, the ICP algorithm provides for each point on the defective part the index of the closest point on the noise-free CAD model, so that deviations from target can be computed as the Euclidean distance between them (this is the minimum distance \min_{j=1,...,m_{2}}C(\mbox{\boldmath\Gamma}\mbox{\boldmathx}_{q,i}+\mbox{\boldmath\gamma},\mbox{\boldmathx}_{p,j})=\min_{j=1,...,m_{2}}||\mbox{\boldmath\Gamma}\mbox{\boldmathx}_{q,i}+\mbox{\boldmath\gamma}-\mbox{\boldmathx}_{p,j}|| in problem (2), with {\mbox{\boldmathX}}_{q} and {\mbox{\boldmathX}}_{p} being the OC part and the CAD model respectively).

Figure 12 shows three different locally defective parts, each with different number of points and with isotropic errors $N\sim(0,0.05^{2})$ added to all three coordinates. We color each point on the OC part proportionally to these deviations, with lighter colors corresponding to larger deviations. As it can be seen, the location of each of the 3 defects on a part, the two parts with chipped corners and the part with a protrusion in one of its “teeth”, is very accurately identified. We suggest to conduct this ICP localization diagnostic after each SPC alarm.

When the change in shape is rather small in the sense that point-wise deviations do not increase sharply, the ICP localization diagnostic will not work as well as in the case where the change of shape is very evident. To illustrate, Figure 11 shows four out-of-control cylindrical parts with increasing $\delta$ values. The number of points varies from part to part and the same isotropic noise we have been using is added. Each part is color coded by the deviation from CAD target and lighter colors means larger deviations. The global deformation of the cylinders is strongest along the “waist” of the cylinder, and is only detectable by the ICP registration when $\delta$ is large enough to be quite evident to the eye. This is consistent with our findings in the OC run length analysis, where the ICP objective is more effective with relatively larger $\delta$ values. A similar “defect localization” could be performed with a registration method based on other distance functions between points, such as spectral distances, see Appendix D (supplementary materials, on-line).

An additional diagnostic worth mentioning is the nonparametric estimator of the change point suggested by Chen *et al. *, (2016), which works from the test statistics used by their DFEWMA chart and evidently can be used in the present situation as well, if desired, to determine the first part in the sequence that needs to be investigated with the ICP diagnostic.

6 A PERMUTATION-BASED SPC SCHEME FOR “PHASE I”

Although the chart by Chen *et al. *, (2016) is self-starting, these authors recommend to perform a “Phase I” of $m_{0}$ parts to avoid masking effects in their chart if a problem occurs very early after startup, with recommended values of $m_{0}$ of at least 50. If parts are expensive, it is important to have an additional scheme that can detect an out of control process within this period with high probability. We therefore suggest using the distribution-free multivariate chart proposed by Capizzi & Masarotto, (2017), available in their R package dfphase1, applied to the lower spectrum of the LB operator. Similar to the DFEWMA chart, these authors used rank related statistics, so their method is more sensitive to small rather than large changes in the process. Following Capizzi & Masarotto, (2017), the performance of the chart is evaluated based on the false alarm probability (FAP) for an in-control process and the alarm probabilities for out-of-control scenarios. Given the “Phase I” dataset consisting of the cropped spectra, we first discard the first eigenvalue of each part because it is theoretically zero (non-zero first eigenvalues from the estimated LB operator are pure numerical error). We compare the performance when varying the number of eigenvalues used. The same prototype part and cylinder as in Section 4.3 is used.

Table 6 summarizes the in-control false alarm probability (FAP), where similarly to Section 4.3, we sampled and permuted 50 IC parts from pre-simulated 40,000 IC parts instead of simulating new parts for new replications to reduce the computational effort. As the table shows, all cases have a FAP close to the nominal $\alpha=0.05$, indicating that as long as the FAP is concerned, the method works well regardless the number of eigenvalues selected.

The out-of-control alarm probabilities for the cylindrical parts with increasing “barrel” shape parametrized by parameter $\delta$ are shown in Table 7. Each replication consists of 25 simulated IC parts followed by 25 simulated OC parts under the same isotropic noise as before (N$(0,0.05^{2}{\bf I}_{3})$). The table indicates how the power to detect changes in the shape of the cylinders is concentrated in the lower part of the spectrum. Using up to the 100th eigenvalue is counterproductive: the detection capability goes down as the higher eigenvalues are associated with geometrical noise. Using the first 15 eigenvalues, as was also recommended for Phase II, provides good detection power in Phase I as well.

The OC alarm probabilities for the defective parts displayed in Figure 8 are shown in Table 8. Each new set of “Phase 1” data consists of 25 simulated IC parts and 25 simulated OC parts with the same isotropic noise level as before. As it can be seen, the protrusion defect is the easiest to detect, followed by chipped #2 defect, with the chipped #1 type of defect being the hardest to detect. This is consistent with our OC results from “Phase II” (Table 3). Similarly to the cylindrical-barrel defect parts, using the top 15 eigenvalues provides good detection capabilities, unless the change to detect is small, which is the case of the chipped parts. In such case, using up tp the 75th eigenvalue provides better detection, but again, adding eigenvalues up to the 100th is counterproductive due to their modeling of noise.

7 CONCLUSIONS AND FURTHER WORK

We have presented a fundamentally new approach for the Statistical Process Control of discrete-part manufacturing processes based on intrinsic geometrical properties of the sequence of parts that, contrary to existing methods, does not require registration of the parts and does not require equal number of points per part. Our proposal brings SPC closer to Computer Graphics/Vision methods. The SPC problem, however, is inherently different than the shape similarity problem from these other fields, given that contrary to them, a method to be useful for SPC must be able to distinguish small but significant shape and size differences in a sequence of very similar parts measured with noise, avoiding false alarms but considering increments in noise a potential additional source of an out-of-control condition. In contrast, computer graphics/vision methods for shape similarity assessment typically aim to detect large shape differences in a manner that is robust with respect to any measurement noise, which (if existent) is filtered out, and usually neglect differences in size.

The main differential-geometric tool we use is the unnormalized spectrum of the discrete Laplace-Beltrami operator, cropped to consider only its lower part. We discussed two different discrete LB operator approximations which are symmetrizable (ensuring a real spectrum and providing computational advantages) and pointwise convergent (providing theoretical guarantees), and adopted the localized mesh Laplacian of Li *et al. *, (2015) due to its sparseness. Other discrete approximations of the LB operator are also symmetrizable, based on Finite Element Methods (FEM, Reuter *et al. *, 2006), and we leave their study and comparison with the Localized Mesh Laplacian used here for future work. The LB-spectrum chart method is intrinsic and hence avoids registration of the parts, which is a hard-to-solve combinatorial problem.

Given the non-normality of the discrete LB spectrum, we proposed to use (with some modifications) a multivariate, nonparametric permutation-based control chart due to Chen *et al. *, (2016) for on-line or “Phase II” SPC and a similar permutation and rank-based approach for the startup or ”Phase I” due to Capizzi & Masarotto, (2017). Run length analyses and detection probability assessments, respectively, indicate the practical feasibility of the methods, even with relative large meshes (with tens of thousands of points) on a modest desktop computer assuming enough storage. The on-line (Phase II) method is especially sensitive to detect small changes in the shape or size of the surfaces, while providing an easy to tune in-control Average Run Length. The Phase I method applied to only the first several eigenvalues has excellent detection performance while controlling the false alarm probability. We compared our Phase II method with a nonparametric univariate chart based on registration of the parts using the Iterated Closest Point (ICP) algorithm, considering its objective function as a monitoring statistic, but found the LB-spectrum chart to be much more sensitive to detect process changes. An ICP-based method was presented to determine the localization of the defect on the part surface as a post-alarm diagnostic only to be used after the generic, or overall SPC mechanism provided by the LB spectrum chart, triggers an alarm. Phase II run length performance comparisons were made also for isotropic and non-isotropic noise, and further comparisons with a registration-based SPC method shows very competitive behavior for our method, even under conditions that clearly favor registration methods (small, equal size lattice meshes).

We focused in this paper on surface data, where the intrinsic dimension of the manifold is 2 and the topology of the object has genus 0. Our methods carry over to the case of voxel data, where the intrinsic dimension of the manifold is 3 and objects with holes and internal features can be modeled. These are of particular interest in SPC of additive manufacturing data obtained via computed tomography scans of a part, in order to determine the inner features of printed parts. The intrinsically 3-manifold data can be represented with a tetrahedralization, and FEM methods exist for approximating the LB operator from such data structure (Reuter *et al. *,, 2007), thus we will consider voxel extensions of our methods in future work.

Further work is also needed to develop charts that detect not only changes in the mean geometry of a part, but also changes in the overall variance of noise. Our methods assumed no systematic local bias due to optical aberration in the scanner, and extensions to deal with such bias, if significant, are of interest if calibration is not efficient.

Acknowledgements.- We thank four anonymous referees, an Associate Editor, and the Editor (Prof. Roshan Joseph) for several comments and suggestions which have resulted in an improved manuscript. We also thank Professors Bianca Colosimo (Politecnico di Milano, Italy) and Massimo Pacella (Università del Salento, Italy) for their kindness in making available the data and code from their paper (Colosimo *et al. *,, 2014).

SUPPLEMENTARY MATERIALS

Appendices

: A) Proofs; B) Exact moments of the DFEWMA chart statistics; C) Further discussion on the heat equation and the LB operator of manifolds and graphs; D) Other intrinsic geometrical statistics and their use for SPC.

Matlab code

for the computation of the Localized Mesh Laplacian and for the modified DFEWMA control chart in the examples.

Data set:

Prototypical part CAD model and CAD model for cylinder (mesh data), both in-control noise-free and noise-free defect versions included.

Appendix A Appendix. Proofs.

A.1 Proof of Lemma 1. Mesh Laplacian derivation

We show how the discrete Mesh Laplacian can be written as the matrix $L_{k}^{t}=D-W$. Define $\mathcal{K}$=given triangulation, $T_{\mathcal{K}}$=set of all triangles in $\mathcal{K}$, $V_{\mathcal{K}}$=set of all vertices in $\mathcal{K}$, $A(T),V(T)$=area and set of vertices of triangle $T$, respectively. The discrete Mesh Laplacian can be written as follows:

[TABLE]

For any function $f$, we first discretize it based on the points in the mesh as a vector $f=\begin{pmatrix}f(p_{1})&f(p_{2})&\cdots&f(p_{m})\end{pmatrix}^{T}$. Then let $W_{ij}=\frac{1}{12\pi t^{2}}\sum_{T:p_{j}\in V(T)}A(T)e^{-||p_{i}-p_{j}||^{2}/4t}$ and $D_{ii}=\sum_{j}W_{ij}$, we have

[TABLE]

Define $D=\mbox{diag}[D_{ii}]$ and $W=[W_{ij}]$, then $L^{t}_{\mathcal{K}}f=(D-W)f$ and $L^{t}_{K}=D-W$. This is in agreement with the Graph Laplacian, defined as $L=D-W$ where the weights $W$ are given by the adjacency matrix of the graph. $\;\;\blacksquare$

A.2 Proof of Lemma 2, convergence of the spectrum of the Localized Mesh Laplacian.

Here we prove the convergence of the spectrum of the Localized Mesh Laplacian (9) as the mesh $\mathcal{K}$ gets finer and closely follows the manifold $\mathcal{M}$, a result that in itself is new but follows directly from previous results in Belkin *et al. *, (2008), Li *et al. *, (2015) and Dey *et al. *, (2010). Following these authors, denote by $\rho(\mathcal{M})$ the reach of the manifold, defined as the infimum of the closest distance from any point $p\in\mathcal{M}$ and the medial axis of $\mathcal{M}$. Define an $({\epsilon,\eta})$ approximation mesh to a manifold $\mathcal{M}$ as a triangulation $\mathcal{K}$ where: 1) for any point $p_{i}\in\mathcal{M}$, there is vertex in $\mathcal{K}$ that is at most $\epsilon\rho(\mathcal{M})$ away, and 2) for a face $t\in\mathcal{K}$ and vertex $p\in t$, the angle between the normals at $p$ and $t$ is at most $\eta$, where $\eta=O(\epsilon)$. Condition one assures the mesh is fine enough and condition 2 ensures the distortion between $\mathcal{K}$ and $\mathcal{M}$ is small.

Consider an $(\epsilon,\eta)$ mesh approximation $\mathcal{K}$ to a smooth surface $\mathcal{M}\subset\mathbb{R}^{3}$. In Dey *et al. *, (2010), an intermediate operator $\mathcal{L}_{t}^{\mathcal{M}}$ defined as

[TABLE]

is used to connect the Mesh Laplacian $L^{t}_{\mathcal{K}}$ in (7) and the manifold LB operator $\Delta_{\mathcal{M}}$. The authors proved the pointwise convergence between $\mathcal{L}_{t}^{\mathcal{M}}f$ and $\Delta_{\mathcal{M}}f$, and between $L^{t}_{\mathcal{K}}f$ and $\mathcal{L}_{t}^{\mathcal{M}}f$ in sequence, and finally arrived at the convergence of the spectrum of $L^{t}_{\mathcal{K}}$ to the true LB spectrum as $\epsilon\rightarrow 0,t\rightarrow 0$, and $\epsilon/t^{4}\rightarrow 0$.

As Dey *et al. *, (2010) indicate, to also show convergence between the corresponding spectra of $L_{\mathcal{K}}^{t}$ and $\Delta_{\mathcal{M}}$, a stronger (than pointwise) convergence is needed between the intermediate operator $\mathcal{L}_{t}^{\mathcal{M}}$ and the Mesh Laplacian $L^{t}_{\mathcal{K}}$. These authors show convergence between these operators in a Sobolev 2-norm, defined for a function $f\in H_{2}$ (Sobolev space) by $||f||_{H_{2}}=\sqrt{||f^{(1)}+f^{(2)}||^{2}}$. Given the small differences between the Mesh Laplacian (7) and the Localized Mesh Laplacian (9), we follow their proof using the latter Laplacian, which requires geodesic distances on $\mathcal{M}$.

First, we modify the intermediate operator $\mathcal{L}_{t}^{\mathcal{M}}$ to use the geodesic distances instead of the Euclidean distances in the exponential term. Theorem 1 in Li *et al. *, (2015) proved the pointwise convergence between the new intermediate operator and the true LB operator as $t\to 0$. To show convergence in $H_{2}$ between the Localized Mesh Laplacian and the new intermediate operator and complete the proof, we need to substitute the geodesic distance for the Euclidean distance in the heat kernel $G_{t}(x,y)\triangleq e^{-\frac{[\text{dist}(x,y)]^{2}}{4t}}$, whose structure is explicitly used only in Lemma 4.3 in Dey *et al. *, (2010), which asserts that the derivatives needed in the $H_{2}$ norm have the form:

[TABLE]

Dey *et al. *, (2010) proved this Lemma for $\text{dist}(x,y)=||x-y||$. We need to prove it for $\text{dist}(x,y)=g(x,y)$, the geodesic distance between $x$ and $y$, which can be done in 2 steps, as follows.

First we observe that

[TABLE]

where the term $\frac{[\text{dist}(x,y)]^{i-2j}}{(2t)^{i-j}}G_{t}(x,y)$ has exactly the same format (and also the same power) regardless the distance used. This equivalence remains true for higher order derivatives, because this structure (and also the power in $O(i^{i})$) comes from differentiating an exponential term (namely $G_{t}(x,y)$) repeatedly, regardless of which distance is used.

Secondly, we consider the difference between $\frac{\partial||x-y||}{\partial x}$ and $\frac{\partial g(x,y)}{\partial x}$. Since small neighborhoods in Riemannian manifolds can be seen as Euclidean spaces, $\frac{\partial||x-y||}{\partial x}$ and $\frac{\partial g(x,y)}{\partial x}$ at least have the same magnitude, and the exact ratio can be omitted in the $O(i^{i})$ term. Putting these facts together shows how Lemma 4.3 also holds for the geodesic distance as well. Once this is done, all remaining steps in the proof of Dey *et al. *, (2010) can be exactly followed to obtain the convergence of the spectrum of the Localized Mesh Laplacian to the true LB spectrum under the same limit conditions, namely, $\epsilon\rightarrow 0,t\rightarrow 0$, and $\epsilon/t^{4}\rightarrow 0$. $\;\;\blacksquare$

A.3 Proof of Proposition 1, on the real spectrum of the LB operator

Here we prove the following:

Proposition 1. The Mesh Laplacian (7) and the Localized Mesh Laplacian (9) can be written as (10) and therefore their eigenvalues are all real.

Proof. From Lemma 1, both the Mesh Laplacian and the Localized Mesh Laplacian can be written as $L^{t}_{\mathcal{K}}=D-W$, where $W_{ij}=\frac{1}{12\pi t^{2}}A(p_{i},p_{j})e^{-dist(p_{i},p_{j})^{2}/4t}$ and $D_{jj}=\sum_{i}W_{ij}$ is diagonal, for $i=1,...,m$ and $j=1,...,m$. The only differences between these two discrete Laplacians are the definitions of area $A(p_{i},p_{j})$ and distance $dist(p_{i},p_{j})$, which are summarized in the following table:

[TABLE]

Where recall $A(p_{j})$ is the area of the one ring neighborhood of point $p_{j}$. Now define a diagonal matrix $B$ where $B_{ii}=A(p_{i})$, which is obviously positive definite. Pre-multiplying a matrix by $B$ is the same as multiplying the $i$th row of the matrix times the $i$th diagonal element in $B$, $i=1,...,m$, so matrix $BW$ has elements:

[TABLE]

It is easy to see that $A(p_{i})A(p_{i},p_{j})=A(p_{j})A(p_{j},p_{i})$ holds with either of the two area definitions above, and indices $i$ and $j$ are interchangeable in both types of distances. Therefore, $BW$ is symmetric for both types of discrete Laplacians. In addition, $BD$ is a diagonal matrix and automatically symmetric. So $BL^{t}_{\mathcal{K}}=BD-BW{\triangleq}L$ is symmetric, and $L^{t}_{\mathcal{K}}=B^{-1}L$ is the multiplication of two symmetric, positive semi-definite ($B$ is positive definite) matrices. Therefore, both discrete Laplacians, (7) and (9) are symmetrizable and thus always have real eigenvalues. $\blacksquare$

Appendix B Appendix. Exact moments for DFEWMA chart statistic

Here we give the derivation of the correct moments in the statistic (13). Since this is a weighted sum, the covariances between the terms added need to be considered222We point out that the expression for $T_{jn}(w,\lambda)$ in Chen *et al. *, (2016) has an extra $w$ in the numerator (the mean), and the denominator (the variance) is only correct for $\lambda=0$..

First consider $R_{jni}$ as defined in Section 4.2. Under the null hypothesis, all parts come from the same distribution. Thus, $R_{jni}$ can be any number from 1 to $m_{0}+n$ with equal probability, and therefore, it is easy to show that

[TABLE]

Now consider $\sum_{i=n-w+1}^{n}(1-\lambda)^{n-i}R_{jni}$. Its expectation is:

[TABLE]

When $\lambda=0$,

[TABLE]

When $\lambda\neq 0$,

[TABLE]

where $\lim_{\lambda\rightarrow 0}\;L_{1}(\lambda,w)=w$, i.e.

[TABLE]

Next we look at the variance:

[TABLE]

In particular, let $N=m_{0}+n$ and consider:

[TABLE]

[TABLE]

[TABLE]

When $\lambda=0$, the variance reduces to:

[TABLE]

When $\lambda\neq 0$,

[TABLE]

where $\lim_{\lambda\rightarrow 0}L_{2}(\lambda,w)=w$ and $\lim_{\lambda\rightarrow 0}L_{3}(\lambda,w)=\frac{w(w-1)}{2}$, i.e.

[TABLE]

Summarizing, the expectation equals to:

[TABLE]

and the variance equals to:

[TABLE]

Appendix C Appendix: further discussion on the LB operator of manifolds and graphs.

C.1 The heat equation and the LB operator

The LB operator emerges from the spatial part of the solution to the heat and wave partial differential equations on a manifold $\mathcal{M}$ (Evans,, 2010). The intuition is that to model both heat and wave phenomena on a manifold it is necessary to consider the geometry or curvature of $\mathcal{M}$, and the LB operator encodes it. A heat diffusion process over $\mathcal{M}$, for our purposes, the surface of a 3D object, is governed by the heat equation (or diffusion equation):

[TABLE]

where $\Delta_{\mathcal{M}}$ is the LB operator on $\mathcal{M}$ and $u(x,t)$ is the temperature at $x$ at time $t$, assumed enough times differentiable in each argument. The heat equation intuitively says that the rate of change in time of the temperature is proportional to the “curvature” of $u$ in $\mathcal{M}$.

Separation of variables $u(x,t)=f(x)g(t)$ in the heat equation leads to $f(x)g^{\prime}(t)=\Delta f(x)g(t)$, which, dividing both sides by $f(x)g(t)$ we get $\frac{g^{\prime}(t)}{g(t)}=\frac{\Delta f(x)}{f(x)}=\lambda$ from which we obtain the so-called Helmholtz differential equation (5):

[TABLE]

Similarly, the Helmholtz equation can also originate from the spatial part of the solution to the wave equation $\Delta u(x,t)=\frac{\partial^{2}u(x,t)}{\partial t^{2}}$.

Given an initial heat distribution $u(x,0)=f:\mathcal{M}\rightarrow\mathbb{R}$, the complete or fundamental space-time solution to the heat equation is given by the integral equation (Evans,, 2010, pg. 46):

[TABLE]

where the bilinear function $k_{t}(x,y):\mathbb{R}^{+}\times\mathcal{M}\times\mathcal{M}\rightarrow\mathbb{R}$ is called the heat kernel which can be thought of as the amount of heat transmitted from $x$ to $y$ in time $t$ if initially there is a unit heat source at $x$ and where $dy$ stands for area (of $\mathcal{M}$). $H_{t}$ is called the heat operator. If $\mathcal{M}$ has boundaries, the solution to the heat and wave equations requires additional conditions such as the so-called Dirichlet boundary condition $u(x,t)=0$ for all $x\in\partial\mathcal{M}$(so the boundary $\partial\mathcal{M}$ acts as an absolute refrigerator for all $t>0$ in the case of the heat equation).

In case of a flat 2D space (i.e., when $\mathcal{M}=\mathbb{R}^{2}$), the heat kernel is given by the Gaussian function (Belkin & Niyogi,, 2008):

[TABLE]

(in the more general manifold case, when $\mathcal{M}\subset\mathbb{R}^{n}$, the heat kernel can be obtained only after obtaining the spectrum of the LB operator, using expression (18) shown in Appendix D). Note how the Mesh Laplacian (7) is based on the heat kernel.

If $\mathcal{M}=\mathbb{R}^{2}$, the eigenfunctions $f(u,v)_{i}$ satisfying the Helmholtz differential equation ($\Delta_{\mathcal{M}}f=\lambda f$) can be interpreted as the modes of vibration of a 2D drum membrane with resonant or natural frequencies $\lambda_{i}$. The question of whether one can determine the shape of the drum from its spectrum entertained mathematical physicists for many years until it was shown by Gordon *et al. *, (1992) that there are pairs of 2D drums that have the same spectrum, yet their shapes are different. In general, for $\mathcal{M}\subset\mathbb{R}^{3}$, this result implies that the geometric information contained in the spectrum is not enough to uniquely identify the shape of the surface $\mathcal{M}$. However, these pairs of figures are always concave polygons and correspond segment by segment, so they are rare.

C.2 Relation between the mesh Laplacians and the combinatorial graph Laplacian

The discrete Mesh Laplacian in equation (7) and the Localized Mesh Laplacian (9), have the same form as the combinatorial Graph Laplacian (Chung & Lu,, 2006), $L=D-A_{\mathcal{K}}$, where $A_{\mathcal{K}}$ is the adjacency matrix of the mesh $\mathcal{K}$ and $D$ is a diagonal matrix with entries equal to the degree of each vertex. The graph Laplacian is an operator applied to a function defined on the vertices of the graph. It corresponds to the discrete Laplacian if the edge weights are $W_{ij}=1$, so the discrete LB operator (7) can be thought of as a weighted version of the combinatorial Laplacian where the edge weights are given by the integrated heat kernel (with either euclidean or geodesic distances over the neighborhood of a point, depending on the approximation used) over one third of the area of the one-ring neighborhood of either of its two end points. The graph Laplacian arises when modeling a diffusion process on a network or mesh. Let $\psi_{i}$ be the amount of some substance at vertex $i$, and $c(\psi_{j}-\psi_{i})$ be the rate at which the substance flows from vertex $j$ to vertex $i$, where $c$ is a constant. Then

[TABLE]

where $A_{ij}$ is the $(i,j)$th element of $A_{\mathcal{K}}$. It is easy to show (Newman,, 2010) that

[TABLE]

where $\psi$ is a vector with all $\psi_{i}$’s, from which the diffusion equation on a network results, compare to the heat equation (17).

The spectrum of the combinatorial graph Laplacian $L$ has some properties with counterparts in the approximated spectrum of $L_{\mathcal{K}}^{t}$ as the latter is also based on a mesh or network: first, $\lambda_{1}=0$, and the algebraic multiplicity of this eigenvalue gives the number of connected components in the graph. The first eigenvector is constant, and the second eigenvalue $\lambda_{2}$ is greater than zero only if the graph is connected (so $\lambda_{1}$ is not repeated). The signs in the second eigenvector, called Fiedler’s vector, can be used to cluster the vertices in two sets, a notion related to the nodal sets of the eigenvectors of the LB operator (see Chung & Lu, 2006).

C.3 Computation of the true LB operator for a parametric surface

The following material can be found in Differential Geometry books (e.g. Kreyszig, 1991). Here we illustrate with an example the first fundamental form and the LB operator for a parametric surface ($k=2$). In this case, (see Definition 3) the first fundamental form is given by the quadratic form:

[TABLE]

Consider a Torus in $\mathbb{R}^{3}$ described parametrically by:

[TABLE]

where the surface coordinates are $x^{1}=u$ and $x^{2}=v$ (see Figure 13) and where $R$ is the outer radius and $r$ is the inner radius. Then, $g_{11}=(R+r\cos v)^{2}$, $g_{12}=0$ and $g_{22}=r^{2}$, and the first fundamental form is:

[TABLE]

Therefore,

[TABLE]

From equation (4),

[TABLE]

Consider the function $f=u^{2}$ defined on the Torus $\mathcal{M}$. The LB operator applied on $f$ equals

[TABLE]

which is only a function of $v$ (Figure 13, middle). If $f=u\cdot v$ instead, its Laplacian is:

[TABLE]

which varies as a function of both angles (Figure 13, right).

Appendix D Appendix: Other intrinsic geometrical statistics and their use for SPC.

In addition to the spectrum of the LB operator, there are several other descriptors of shape used for object recognition in Computer Vision. Here we review some of these alternatives and report on their applicability to the SPC problem. We focus only on intrinsic geometrical properties, which can be computed for each object without reference to the ambient space, and hence in principle could be used for SPC schemes that do not require registration. We also consider these statistics as further diagnostics in case an alarm has been triggered, similar to the ICP diagnostic presented in Section 5.

D.1 Heat kernel and diffusion distances

For a long time, researchers in computer vision have tried to find shape descriptors that use the properties of a cloud of points from the point of view of each point in the cloud. In a highly cited paper, Belongie *et al. *, (2002) introduce the concept of a shape context, a local description of the shape in the vicinity of a given point, and use it for registration of objects. The main idea is to measure the frequency and location of other points in the neighborhood of each point of a shape and use differences between these measures as costs to be minimized in a classical weighted matching problem (thus points with similar local information tend to be matched), solvable via Linear Programming. For a point $i$ in an object, Belongie *et al. *, (2002) propose to compute a 2-dimensional histogram where the number of points nearby are counted. If $r$ is the Euclidean distance between two points of the shape, the 2-dimensional histogram $h_{i}(\cdot,\cdot)$ extends along $log\;r$ and $\theta$, measuring the distance and direction where the nearby points are located. The 2-dimensional histogram formed by the frequencies $h_{i}(\cdot,\cdot)$ is called the “context” of point $i$ by these authors and was computed for recognition of 2-dimensional shapes via registration. Extending this concept to 3D objects is harder as higher dimensional “histograms” need to be computed at each point to locally describe the points in a neighborhood of a given point.

In contrast, the lower part of the spectrum of the LB operator provides a global description of the shape, with more details included in higher eigenvalue-eigenvector pairs. Interestingly, the geometric interpretation in expression (3) of the LB operator indicates a similarity with the point histograms of Belongie *et al. *, (2002) in the sense that it generates a local characterization at each point on the manifold (the normal whose length relates to the mean curvature). Using the spectrum of the LB operator results in a very effective summary of this local information, providing a global description of the geometry (and even providing some topological information, as discussed in Section 2) of the part. As discussed, the LB spectrum results from the solution of the Helmholtz equation, which relates to the spatial part of the heat equation. It seems possible to also look at the temporal solution of the heat equation to try to find additional local intrinsic descriptors of the geometry of an object, a matter for further work.

A result due to Huber (see e.g., Buser, 1992) proves how two Riemannian surfaces have the same sequence of eigenvalues of their LB operator if and only if they have the same sequence of lengths of their closed geodesics, thus, in principle, the geodesics, evidently intrinsic, could be used for detecting part to part differences in a SPC scheme. For our prototypical part, Figure 14 illustrates some geodesics between a point and some others in the part under consideration, as well as the geodesic distances between the point and all others on the mesh defining the part. Geodesics, being intrinsic to the surface, are simpler to compute than the spectrum, but in practice, the question, if one wishes to apply Huber’s theorem, is that geodesics are greatly affected by noise (Lee & Wasserman,, 2010) and it is not clear how many geodesic distances are needed to compare between 2 objects, a manufactured/scanned one, and a noise-free CAD model, in order to find local defects on the surface of the manufactured part. In this section, we report our experiments with geodesic distances and other more robust intrinsic distance functions based on the heat kernel defined on the surface of a part.

Given the spectrum $\{\lambda_{i}\}$ and associated eigenfunctions $\{\phi_{i}()\}$ of the LB operator on a general Riemannian manifold $\mathcal{M}$ (not necessarily flat euclidean space), the heat kernel is given by the decomposition:

[TABLE]

which, as mentioned before, satisfies the heat equation. For small values of $t$, $k_{t}(x,y)$ reflects local properties of $\mathcal{M}$ around point $x$; for larger values of $t$ it reflects the global structure of $\mathcal{M}$ from the point of view of $x$. It is, in a sense, also conveying analogous information as Belongie *et al. *, (2002) “shape contexts”. There is a close relationship between the heat kernel and geodesic distances $g(x,y)$ for $x,y\in\mathcal{M}$:

[TABLE]

(see Sun *et al. *, 2009a), which can be shown directly by taking the limit to $k_{t}(x,y)\approx\frac{1}{4\pi t}e^{-g(x,y)^{2}/(4t)}$ which is an approximation of (18) for small $t$.

The heat kernel can be also be interpreted as the transition density function of a Brownian motion on $\mathcal{M}$ (Sun *et al. *,, 2009a), or alternatively, the estimated heat kernel, using a discrete LB approximation based on a mesh, can be related to the transition probabilities of a random walk process defined on the mesh, a property used in computer graphics by some authors (see, e.g., Sinha & Ramani, 2014) for shape similarity assessment.

Sun *et al. *, (2009a) propose to use the “autodiffusion” $k_{t}(x,x)$ at different points on a mesh as a “Heat Kernel Signature” used to identify shapes of objects. A useful notion related to the heat kernel is the diffusion distance between $x,y\in\mathcal{M}$ defined by the interplay of heat from $x$ to $y$ and from $y$ to $x$:

[TABLE]

The diffusion distances are intrinsic and more robust to changes in the mesh as they are an “average” distance over all possible paths between two points on the mesh, and are not only the shortest path as the geodesic distances are, a property that has received recent attention in the Statistics and Manifold Learning literature (e.g., see Lee & Wasserman, 2010). This property is shared by heat kernels. This might be a disadvantage for SPC: changes in the mesh may reflect true changes in the shape of the parts that need to be detected. Figure 15 illustrates the evolution, for different values of $t$, of the heat kernels and diffusion distances from a point on the prototype part shown in earlier figures, given the LB spectrum of its mesh. These are the quantities that were computed in the histograms of the last two columns of Figure 16, whose quantile-quantile (Q-Q) plots appear in Figure 17 (below). While the heat kernel can be easily understood as the diffusion in time of heat from a point, the diffusion distances have no easy interpretation from their time evolution plot shown in the figure.

Consider the three parts displayed on the leftmost column of Figure 16 (measurements in mm.), where (as before) we have added isotropic $N(0,0.05^{2}{\bf I}_{3})$ noise to the points of a CAD model defined as a triangulation on a (rather small) point cloud with 1680 points. The euclidean distances between all points for the in-control and defective parts are displayed in the second column; these histograms are invariant, but not intrinsic. The geodesic distances between all points, computed on the triangulation (using the “Fast marching algorithm” on the triangulation, see Kimmel, 2004) is displayed in the third column; this is both invariant and intrinsic.

Figure 17 shows Q-Q plots for an in-control part and two of the types of local defects studied before, the chipped part #1 shown in Figure 16, and the part with a protrusion on the edge at the top. Here the “in control” distributions of the heat kernels and diffusion distances are plotted against the distributions of the defective parts, so the third row shows in-control vs. in-control distributions. We can see deviations from linearity in the profiles of both heat kernels and diffusion distances, but it is not clear what the role of the $t$ parameter in the “detection” capabilities of a profile monitoring scheme would be. It seems plausible that smaller values of $t$ allow to detect local changes better, while larger values of $t$ would detect global changes better, but considerable more research is needed to assess this statement in order to provide useful tools for SPC. There are also deviations even in the third row, where an in-control part is plotted against another in-control part. Kolmogorov-Smirnov tests performed on all the cases in the figure are recognized as significantly different regardless of the type of the parts (in control or defective). This occurs due to the number of pairs of heat kernel or diffusion distances is so large (proportional to $n^{2}$, with $n$ being the number of points), that any slight difference is detected as significant.

D.2 Other spectral distances

Both the heat kernel and the diffusion distances, and the other spectral distances we describe next, can be used as an alternative criterion for registering two parts, using the “shape context” idea of Belongie *et al. *, (2002) to find matches of points between two surfaces. The matches between points in a part that triggered an alarm and a CAD model could be then used as an additional post-alarm SPC diagnostic to determine where the produced part differs from the CAD model, analogously to the ICP post-alarm diagnostic presented earlier (this would require more computational effort than using ICP, as it requires both finding the spectrum and solving the combinatorial registration problem, but could provide better information specially when the defect is small and not evident to the eye). It is unclear, however, how to use these spectral distance diagnostic tools without recourse to registration or how to set in practice the $t$ parameter, present in both heat kernel and diffusion distances.

There are, however, spectral distances that eliminate the $t$ parameter. For a point $x$ on a surface $\mathcal{M}$, Rustamov, (2007) defines the “Global Point Signature” (GPS) as the map of $x$ into the infinite dimensional vector

[TABLE]

In practice, the GPS embedding is finite dimensional using the lower part of the discrete LB spectrum and the associated eigenvectors. The idea of embedding a point from a manifold $\mathcal{M}$ into a different dimensional space is a theme common to both Computer Graphics (e..g, “Heat Kernel Signatures” Sun *et al. *, 2009a), where a higher dimensional embedding is sought to define a similarity metric in that space between points originally in $\mathbb{R}^{3}$, and Manifold Learning (see, e.g., Belkin & Niyogi, 2002 “Laplacian eigenmaps”), where the goal is usually to find a lower dimensional embedding for understanding the data manifold structure. In the case of GPS, a distance function between two points on a 2-manifold can be defined by taking the inner product of their embeddings:

[TABLE]

which can be obtained also as the integral of the heat kernel $k_{t}(x,y)$ with respect to $t$. A distance (metric) between $x,y\in\mathcal{M}$ is then given by $\sqrt{G(x,y)}$.

A distance (metric) related to the GPS distance is the so-called “Commute Time Distance”, defined as:

[TABLE]

which is related to the GPS distance since $d_{c}(x,y)^{2}=G(x,x)+G(y,y)-2G(x,y)$. Lipman *et al. *, (2010) indicate that both the GPS and the commute time distances have the disadvantage of diverging when $x=y$ since $\sum_{i=1}^{\infty}1/\lambda_{i}\approx\sum_{i=1}^{\infty}1/i$ diverges, and proposed instead a modification, called “biharmonic distance”:

[TABLE]

which does not diverge when $x=y$ since $\sum_{i=1}^{\infty}1/\lambda_{i}^{2}\approx\sum_{i=1}^{\infty}1/i^{2}=\pi^{2}/6<\infty$. In practice, all these distances need to be computed using the lower part of the spectrum of the discrete LB operator, up to a given eigenvalue-eigenvector pair index. The applicability of biharmonic distances to the SPC problem stems from the absence of a $t$ parameter, necessary in other spectral distances reviewed, and in principle they could be used for SPC by comparing their distribution with that of an in-control part using Q-Q charts as in Figure 17, but further work is necessary to examine this matter.