{'fraction_non_alphanumeric': 0.049312124642631, 'fraction_numerical': 0.019859163037012875, 'mean_word_length': 4.028268210323203, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 53, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonparametric Regression for 3D Point Cloud Learning

Over the past two decades, we have seen an exponentially increased amount of point clouds collected with irregular shapes in various areas. Motivated by the importance of solid modeling for point clouds, we develop a novel and efficient smoothing tool based on multivariate splines over the tetrahedral partitions to extract the underlying signal and build up a 3D solid model from the point cloud. The proposed smoothing method can denoise or deblur the point cloud effectively and provide a multi-resolution reconstruction of the actual signal. In addition, it can handle sparse and irregularly distributed point clouds and recover the underlying trajectory. The proposed smoothing and interpolation method also provides a natural way of numerosity data reduction. Furthermore, we establish the theoretical guarantees of the proposed method. Specifically, we derive the convergence rate and asymptotic normality of the proposed estimator and illustrate that the convergence rate achieves the optimal nonparametric convergence rate. Through extensive simulation studies and a real data example, we demonstrate the superiority of the proposed method over traditional smoothing methods in terms of estimation accuracy and efficiency of data reduction.


Introduction
Recent advances in computer and information technology have dramatically boosted the availability of three-dimensional (3D) point clouds in many fields, such as geography, environmental science, computer graphics, engineering, economics, and medical imaging. These point clouds are usually collections of enormous measurements in space defined by a given coordinates system. Sometimes, these point clouds are rendered and inspected directly; examples can be found in Rusinkiewicz and Levoy (2000). More often, the point clouds are converted to polygon mesh or triangle mesh models, e.g., 3D reconstruction is widely used for automated driving vehicles Shao, 2021). Many efforts have been devoted to investigating the converting technique, and a survey of surface reconstruction from point clouds can be found in Berger et al. (2017).
These 3D data typically contain more information than the shape of the object as they represent some meaningful values in real life. For example, healthcare professionals can use brain activity levels contained in 3D neuroimages to access new angles, resolutions, and details to understand the tissue better. Another example is using the cloud-based 3D scans to explore the location and access the grade and type of ore in the mining industry. Unfortunately, little exploration has focused on identifying the underlying signal from 3D solid objects. There are many challenges when trying to achieve this goal.
First of all, as shown in the point cloud inputs in Figure 1.1, the shape of the point clouds is usually irregular. To extract the underlying signal over the irregular-shaped 3D solid object, many conventional methods, such as kernel smoothing (Zhu et al., 2014), tensor product smoothing (Reiss and Ogden, 2010), thin plate spline smoothing (Duchon, 1976;Yue and Speckman, 2010), and wavelet smoothing (Morris and Carroll, 2006), suffer from the problem of "leakage", which is referred as the inaccurate signal reconstruction across the complex domains; see Ramsay (2002); Wang and Ranalli (2007); Kim and Wang (2021) for more details. Several techniques have been developed in the last decades to tackle the "leakage" problem for twodimensional (2D) smoothing, such as spline smoothing (Wang and Ranalli, 2007;Wilhelm et al., 2016;Lai and Wang, 2013;Awanou et al., 2005), finite element analysis (Arnone et al., 2022;Lindgren et al., 2011), and kernel smoothing (Guo et al., 2010). For 3D images, Chung et al. (2018) applied the discrete heat kernel smoothing, and for irregular-shaped 3D objects, Huo et al. (2020) considered the finite element method.
Second, the information comprised in the 3D point clouds usually contains unwanted noise that obscures the features of the underlying signal of interest. Consequently, denoising is another fundamental problem in the modeling process. Smoothing methods are widely adopted for denoising in many 2D studies (Goldsmith et al., 2014;Kang, 2020). However, the characteristics of 3D point clouds create substantial challenges in denoising and other analysis (Wong et al., 2016). Last but not least, as mentioned before, the size of the point clouds is usually enormous, and thus both the storage and the analysis method are very challenging. For example, a 3D image produced from a standard positron emission tomography (PET) brain scan contains over half a million voxels.
To deal with these prevalent challenges in converting a 3D irregular shaped point cloud to a solid model and estimating the actual underlying function, we propose a novel and efficient smoothing method based on trivariate penalized spline over tetrahedral partitions (TPST). Another appealing aspect of the proposed TPST lies in its potentiality as a numerosity reduction method (Han et al., 2011). Assuming the underlying signal is smooth, we can reduce the massive data to a much smaller representation since we only need to store the model parameter. More specifically, once the spline basis functions are determined, we can keep the coefficients and use them to recover the underlying signal of the point clouds. Thus, the size of the point cloud is reduced by representing the signal of the point cloud via a much smallersized set of spline coefficients. For example, as demonstrated in Section 6, a point cloud of size 50,000 can be represented using a vector of spline coefficients of length 323. Let's take the biomedical image for another example, where the proposed method can reduce an image point cloud of 510,340 voxels to a vector of 4,856 coefficients with the peak signal-to-noise ratio reaching 30.76. Moreover, as a global estimation method with an explicit model expression, spline smoothing techniques are usually computationally efficient since it only requires solving a single linear system (Ma, 2012;Wang et al., 2020b;Wang and Yang, 2009;Liu and Zhao, 2021;Wang et al., 2020a). Finally, based on an appropriate partition, we can reconstruct or recover the underlying function from globally and locally missing data, effectively denoise or deblur the collected data while preserving inherent geometric features or spatial structures and providing a multi-resolution reconstruction of the actual signal. Furthermore, we investigate the theoretical properties of the proposed smoothing method.
To be more specific, we establish the convergence rate of the proposed estimator, which is governed by the fineness of the tetrahedral partition, the trivariate spline degree and smoothness, penalty parameter, and smoothness of the unknown function. We further illustrate that the rate achieves the optimal nonparametric convergence rate when the penalty parameter is zero.
We also derive the asymptotic normality for the proposed spline estimator.
The rest of the paper is structured as follows. In Section 2, we give an overview of the tetrahedral partitions and the methods for constructing the tetrahedral partitions from a 3D point cloud. Once a tetrahedral partition is obtained, we can construct the spline space over this tetrahedral partition. We introduce the penalized spline estimators in Section 3. In Section 4, we present the theoretical results. In Section 5, we discuss implementation details for the proposed method. We then demonstrate the performance of the proposed smoothers on several simulation examples and a real data example in Section 6. Finally, this paper is concluded, and future work is outlined in Section 7. We put all the theoretical details and more detailed introductions in the supplementary material.

Trivariate splines over tetrahedral partitions and basic properties
In this section, we provide a basic framework for tetrahedral partitions and trivariate splines on those partitions. Please refer to Appendix A for more detailed introductions.

Tetrahedral partitions
For any i = 1, . . . , n, let point p i = (x i , y i , z i ) ∈ R 3 be the location of the ith observation in a point cloud of sample size n. We assume there exists a tight polyhedral domain, Ω, that contains all the data locations p i 's of a point cloud. Furthermore, we consider the domains which may have one or multiple holes not containing any observations, for example, the 3D domain shown in Figure 2.1 (b). Tetrahedral partitions are very popular in approximating the domain Ω. In the following, we use T to denote a tetrahedron, that is, a convex hull of four noncoplanar points in R 3 . Then, a collection = {T 1 , . . . , T N } of N tetrahedra is called a tetrahedral partition of Ω = ∪ N h=1 T h , provided that any pair of tetrahedra in intersect at most at a common vertex, along a common edge, or along a common triangular face. An Figure 2.2 (a), where v ı , ı = 1, . . . , 5, are the vertices of the tetrahedra. In contrast, as illustrated in Figure 2.2 (b) for tetrahedra T 3 = v 2 , v 1 , v 3 , v 6 and T 4 = v 2 , v 1 , v 6 , v 4 , = {T 1 , T 3 , T 4 } does not form a tetrahedral partition of Ω, because of the intersections between the pairs of tetrahedra (T 1 , T 3 ) and (T 1 , T 4 ).
Furthermore, let |T | be the length of the longest edge of T , and T be the radius of the largest ball that can be inscribed in T , then the ratio β T := |T |/ T is the shape parameter of T ; see an illustration example in Figure 2.3 (a). In general, the shape parameter β T describes the shape of T . If T is a regular tetrahedron whose six edges are all of the same lengths, then β T = 2 √ 6. For any other tetrahedron, β T > 2 √ 6, and the larger T is, the flatter the tetrahedron T becomes. Based on this shape parameter, we can define that a tetrahedral partition is β-quasi-uniform if there is a positive value β such that satisfies where | | := max{|T |, T ∈ } is referred as the size of , i.e., the length of the longest edge of . Let N be the number of the tetrahedra in the polygonal domain Ω. From (2.1), we have for a β-quasi-uniform partition, N ≤ (4π| | 3 ) −1 3V Ω β 3 , where V Ω denotes the volume of Ω.
To illustrate the definition of quasi-uniform partition, we show an example of non-β-quasiuniform tetrahedral partition in Figure 2.3 (b), where the edges v 4 , c 1 and c 2 , c 3 are perpendicular to the edge v 2 , v 3 , and c 1 , c 2 and c 3 , c 4 are perpendicular to v 3 , v 4 .

Point clouds to tetrahedral partitions
We now describe how to construct a tetrahedral partition given a 3D point cloud. We start with the surface reconstruction process of the point cloud, which converts it into more practical triangular meshes to find the tight domain, Ω. As Ω has piecewise linear boundary faces, we (a) Shape parameters (b) Non-β-quasi-uniform tetrahedral partition Figure 2.3: (a) shows an illustration of the shape parameters of a tetrahedron; (b) shows an example of a non-β-quasi-uniform tetrahedral partition.
have a normal direction perpendicular to each planar face of Ω, which can be used to determine whether a point is inside or outside of Ω. Then, we use this normal direction information of all boundary faces and the observed data locations to partition Ω into a collection of tetrahedra.
Many mesh generation algorithms for 3D domains available in various software packages and toolboxes, such as the MATLAB functions delaunay and distmesh/distmeshnd (Persson and Strang, 2004), the C++ library CGAL (Jamin et al., 2015;Project, 2020), as well as the TetGen (Si, 2015) and iso2mesh (Fang and Boas, 2009) where (b 1 , b 2 , b 3 , b 4 ) are called the barycentric coordinates of p relative to the tetrahedron T , and they are nonnegative if p is inside or on the faces of T . Accordingly, for some nonnegative integers i, j, k, l with i + j + k + l = d, define trivariate Bernstein basis polynomial of degree d relative to T as For any positive integer d, let P d be the space of all trivariate polynomials with total degrees less than or equal to d. Note that the dimension of P d is d+3 3 . According to Theorem 15.8 in Lai and Schumaker (2007) and Lemma A.2 in the Appendix A, the set of Bernstein basis polynomials in (2.2) forms a basis for the space of polynomials P d . Thus, any polynomial φ(p) ∈ P d on with domain T can be written uniquely as the B-form, For any nonnegative integer r, we use C r (Ω) to denote the collection of all r-th continuously differentiable functions over Ω. For tetrahedral partition = {T 1 , . . . , T N }, let S r d ( ) = {s ∈ C r (Ω) : s| T ∈ P d (T ), T ∈ } be a spline space of degree d and smoothness r over , where s| T is the polynomial piece of spline s restricted on tetrahedron T . According to (2.3), for any For noise-free data, the approximation order of trivariate spline spaces was studied in Lai and Schumaker (2007) when d ≥ 8r + 1 and Lai (1989) when d ≥ 6r + 3. Particularly, Lai (1989) proved that when d ≥ 6r + 3, the space S r d ( ) can attain the full approximation power; see Lemma 1 below. For any index α = (α 1 , α 2 , α 3 ) of order |α| = α 1 + α 2 + α 3 , we denote the derivatives D α s = ∂ |α| s/∂ α 1 x ∂ α 2 y ∂ α 3 z . For any function f over the closure of domain Ω, denote the L q (Ω) norm (1 ≤ q < ∞) and supremum norm as f L q (Ω) = { Ω |f (v)| q dv} 1/q and f ∞,Ω = sup v∈Ω |f (v)|, respectively. Next, for any k ≥ 0, denote Lemma 1. (Theorem 3.5.2 in Lai, 1989) For all f ∈ W d+1,q (Ω) with 1 ≤ q ≤ ∞, r ≥ 0, and for all 0 ≤ |α| ≤ m, where K > 0 is a constant independent of f and | | but is dependent on the geometry of .
Remark 1. Although trivariate splines with d ≥ 6r + 3 is required to achieve the full approximation order, we can still use spline space S r d ( ) with d < 6r + 3 to approximate the underlying signal from a given point cloud dataset. Indeed, S r d ( ) is usually not an empty set as long as d ≥ r, e.g., S r d ( ) can be used to reproduce polynomials of degree d. We use r = 1 and various values of d in our numerical experiments.

Penalized spline estimators
With all the preparations in the previous introduction, we apply the trivariate spline over the tetrahedral partition to recognize the underlying signal from 3D point clouds in this section. In the following, we let W i be the response variable observed on the ith location or design point p i = (x i , y i , z i ) ∈ R 3 , i = 1, . . . , n. Then, we regard any point cloud as a set To extract the underlying signal, we consider the following nonparametric regression model where m(·) is some smooth but unknown 3D function, and ε i is the random error term with mean zero.
In nonparametric smoothing, the roughness penalty approach is widely used when smoothing noisy data (Green and Silverman, 1994;Wood, 2003;Lai and Wang, 2013). Including a roughness penalty and choosing a proper tuning parameter can avoid overfitting problems and balance the bias and variance of the estimator of the function. To estimate the underlying function m in (3.1), we formulate the roughness penalty approach as the following penalized least squares problem: where the roughness penalty and ρ n is the roughness penalty parameter. We propose to seek for the minimizer of (3.2) in S r d ( ) and denote it as m ρn , which is the Trivariate Penalized Spline over Tetrahedral partition (TPST) estimator of m. The tuning parameter ρ n controls the smoothness of the fitted spline function. A larger ρ n leads to a less fluctuating s. If ρ n goes to infinity, our estimator shrinks to linear functions where the roughness penalty E(s) = 0. On the other hand, when ρ n = 0, the estimator becomes the standard unpenalized least squares spline estimator. A proper penalty parameter ρ n balances the goodness of fit for the data and the volatility of estimated functions.
The penalty E(s) is a commonly used penalty, the second-order derivative to zero for surface smoothing; see Green and Silverman (1994) for the 2D case. For a spline function s ∈ S r d ( ), combining (2.4), the roughness penalty in (3.3) can be written as follows: where P T is the corresponding penalty matrix. Let B d = (B d (p 1 ), . . . , B d (p n )) , then we have E(B d γ) = γ Pγ, where P is the block diagonal penalty matrix. See Section 5.2 for more details in derivations and calculations of penalty matrices.
Since s ∈ C r ( ), the coefficients of s satisfy some smoothness conditions across each interior faces of . One can obtain a smoothness constraint matrix H such that Hγ = 0 by repeatedly applying the equation (A.8) in the Appendix A over all shared triangular faces, and more details are available in Section A.2 in the supplementary material. A MATLAB implementation is discussed in Awanou et al. (2005). Thus, the objective function can be written as where H is the matrix for all smoothness conditions across shared edges or faces of tetrahedra, which depends on the smoothness parameter r and the structure of the tetrahedral partition.
One can use QR decomposition to get rid of the constraint in (3.4). Specifically, where Q is an orthogonal matrix, R is an upper triangular matrix, R 1 is a full rank matrix with the same rank as H, and R 2 is a matrix of zeros. Note that for any vector γ satisfying Hγ = 0, there exists some θ such that γ = Q 2 θ. Also, for any θ, H(Q 2 θ) = 0 holds. Thus, (3.4) is equivalent to a penalized regression problem without constraint: which leads to a closed form of the solution. To be specific, let W = (W 1 , . . . , W n ) , then the minimizer of (3.5) is given by Consequently, the spline coefficient in (3.4) can be estimated by γ ρn = Q 2 θ ρn , which yields the TPST estimator m ρn (p) = B d (p) γ ρn . Several methods can be used for choosing the penalty parameter ρ n , such as the block cross-validation (block CV) and generalized cross-validation (GCV). Detailed discussion is given in Section 5.3.

Theoretical results
In the previous section, we discuss how to construct TPST and capture the underlying signal using TPST from the point clouds. In the following, we investigate the theoretical support of TPST.
For random variables X n , X n = O P (b n ) if lim c→∞ lim sup n P (|X n | ≥ cb n ) = 0, n ≥ 1.
Similarly, for any constant c > 0, X n = o P (b n ) if lim n→∞ P (|X n | ≥ cb n ) = 0. And a n b n if there exist two positive constants c 1 , c 2 such that c 1 |a n | ≤ |b n | ≤ c 2 |a n |, for all n ≥ 1.

Convergence rate
Before we state our main results, we make the following assumptions, which are standard in nonparametric literature (Lai and Wang, 2013;Yu et al., 2020;Huang, 2003).
(A3) The density function of the observations is bounded below and above.
(A4) The number of the tetrahedra N and the sample size n satisfy that N = Cn γ for some constant C > 0 and γ < 1.
To obtain the asymptotic analysis of spline estimators, we account on an important property that the data-driven norm uniformly approach to its expectation uniformly over the entire spline space. To see this, for any function f over the closure of domain Ω, let E n (f ) = . Define the empirical inner product and norm as f 1 , f 2 n,Ω = E n (f 1 f 2 ) and f 1 2 n,Ω = f 1 , f 1 n,Ω , respectively, for measurable functions f 1 and f 2 on Ω. The theoretical L 2 inner product and the induced norm are given by f 1 , f 2 L 2 (Ω) = E(f 1 f 2 ) and f 1 2 L 2 (Ω) = f 1 , f 1 L 2 (Ω) . We illustrate the uniform convergence rate of empirical inner product to the theoretical one in the following Lemma 2.
Lemma 2. Let {B ξ } ξ∈M be the basis for S r d ( ) constructed in Lai and Schumaker (2007), where M stands for the index set of spline bases. Let g 1 = ξ∈M c ξ B ξ , g 2 = ζ∈Mc ζ B ζ be any spline functions in S r d ( ). Under Assumptions (A3) and (A4), we have over Ω = ∪ T ∈ T equipped with the inner product n f 1 , f 2 n,Ω + ρ n f 1 , f 2 E , where Next, we introduce a measure of the complexity of the spline space S r d ( ), and another measure which bounds the size of the derivatives: g n,Ω , g n,Ω = 0 , V n = sup g∈S r d ( ) g E g n,Ω , g n,Ω = 0 . (4.1) These two measures will play an important role in developing the asymptotic results. We use the following Lemma 3 to demonstrate the upper bounds of V n and V n .
We then define a linear operator P ρn : B → S r d ( ) such that P ρn W = m ρn . Note that in general, P ρn is not a linear projection. Thus, we have P ρn W = P ρn m + P ρn ε, where P ρn m and P ρn ε are the penalized spline estimators based on {m(p i )} n i=1 and {ε i } n i=1 , respectively. Under some conditions (Huang, 2003), P 0 is a bounded operator on S r d ( ) in the maximum norm. Denote s ρn,m = P ρn m and s ρn,ε = P ρn ε. Consequently, for the penalized spline estimator m ρn in (3.2), where s ρn,m (p) − m(p) and s ρn,ε (p) are referred to as the bias and noise terms, respectively.
According to the error decomposition in (4.2), to derive the convergence rate of m ρn to m, it is sufficient to evaluate the size of the bias and noise terms. The following propositions provide the upper bound of the bias size and noise size.
Based on Propositions 1-3, we illustrate the convergence rates of the TPST estimator in the following Theorem 1, in terms of both the L 2 (Ω) and supremum norms.
Remark 3. Assumption (A3) is standard in nonparametric literature. Even though this assumption may not be satisfied with missing data, we usually can find a proper choice of the tetrahedral partition by adjusting the size of some of the tetrahedra to obtain a decent penalized spline fitting.

Asymptotic normality
To derive the asymptotic normality of the TPST estimator, we further assume the following conditions.
(A4') The number of the tetrahedra N and the sample size n satisfy that N = Cn γ for some constant C > 0 and 1/( + 2) < γ < 1.

Remark 4. A sufficient condition for a negligible bias term is provided by Assumptions (A4')
and (A5). Compared with Assumption (A4) in Section 4.1, (A4') further assumes that the number of tetrahedra needs to be greater than a lower bound which depends on the degree of the function. A similar assumption for the univariate case has been discussed in Li and Ruppert (2008). Meanwhile, Assumption (A5) requires smaller ρ n , which reduces the bias through under smoothing.
where P is the collection of all p i , i = 1, . . . , n.
The above asymptotic distribution result can be used to construct asymptotic confidence intervals. For example, if we estimate m(p) using piecewise constant splines, Lemma B.5 in the Appendix B gives the size of the pointwise variance where V T is the volume of tetrahedron T . Therefore, an asymptotic 100(1 − α)% pointwise confidence envelop for m(p) over Ω is where f stands for the density function of p i , i = 1, . . . , n.

Implementation details
This section provides some implementation details on how to construct the penalty matrix P and select the penalty parameter ρ n in (3.4). To facilitate discussion, we first introduce the directional derivatives for basis functions, followed by the construction details of penalty matrix, the selection criteria of penalty parameter, and tetrahedral partition selection.

Directional derivatives for basis functions
To generalize the smoothness restriction over the joint triangular face for two adjacent tetrahedra, we need to introduce the definitions of directional derivative first. For a general multivariate smooth function φ, the directional derivative at point p with respect to direction u is defined as Accordingly, for vector u := (u x , u y , u z ) ∈ R 3 and trivariate function φ, the directional deriva- Direct calculation gives the directional derivative of the Bernstein basis polynomial B d ijkl . Based on the conclusion in (A.7) in the Appendix A, the mth order directional derivatives for all the Bernstein basis functions with degree d can be written as some linear combination of Bernstein basis functions with degree d − m. Specifically, for directions u 1 , . . . , u m , there Therefore, we can obtain the explicit form of the mth order directional derivatives for the Bernstein basis functions with degree d, given the explicit form of matrix C where (a 1 , a 2 , a 3 , a 4 ) is the barycentric coordinate of direction u. Then, the , where e I is unit vector with Ith element being one, I 1 , I 2 , I 3 and I 4 are indexes of basis functions

Details of constructing penalty matrix
Next, we introduce the details of constructing penalty P matrix. Following the definition of penalty matrix P T in the Section 3, we can further implement as Applying the results of directional derivatives in Section 5.1, the second order derivative of B d (p) at directions u 1 and u 2 are For notation simplicity, we denote C α d the derivative coefficient matrix for D α B d (p) with respect to direction(s) of Cartesian coordinate system, i.e., By the Lemma 15.29 in Lai and Schumaker (2007), we have Recall that E(s) = T ∈ E(s T ) and E(s T ) = γ T P T γ T . Therefore, E(s) = γ Pγ, where γ = (γ 1 , . . . , γ N ) , and P = diag(P T , T ∈ ) is a block diagonal matrix.

Penalty parameter selection
To balance the bias and variance of the proposed estimator and achieve a good estimation and prediction performance, it is crucial to choose a suitable value of the penalty parameter ρ n .
Since the in-sample fitting errors can not gauge the prediction accuracy of the fitted function, we select a criterion function that attempts to measure the out-of-sample performance of the fitted model. Minimizing the generalized cross-validation (GCV) criterion is one computationally efficient approach to selecting smoothing parameters that also has good theoretical properties.
Note that 3D object data are often generated with spatial dependence. When performing cross-validation (CV), these dependence structures are usually ignored, leading to underestimating the predictive error (Roberts et al., 2017). To tackle this problem, one can adopt the block CV strategy (Roberts et al., 2017;Valavi et al., 2019). To be specific, all the sample points are first divided into 3D blocks with similar volumes. Then, these blocks are randomly allocated to the CV folds. In this paper, we adopt the tetrahedral partition to divide the 3D domain into small 3D blocks. Each tetrahedron is considered as one single 3D block. Figure   5.1 shows an example of block CV using a tetrahedral partition. In Figure 5.1, we divide the domain into 504 tetrahedra and randomly assign these tetrahedra into five folds with colors indicating different folds. Figure 5.1: An illustration of block CV using a tetrahedral partition. Tetrahedra with the same color belong to the same fold.

Tetrahedral partition selection
As illustrated in Figure 2.4, to form a proper tetrahedral partition, we start with surface reconstruction. A wide range of techniques has been developed to reconstruct the surface from point clouds. Our paper focuses on the methods that generalize a well-sampled point cloud to arbitrary shapes and produce a watertight surface mesh, such as a triangular mesh. Since we are interested in recovering the actual signal within the 3D point clouds, we skip the details of the surface reconstruction here. The triangulation surface representations can be coarse or fine and affect the construction of the tetrahedral partition. A very coarse tetrahedral partition could give a poor approximation to the object, while a very fine tetrahedral partition could introduce a more serious computation burden. Note that Assumption (A4') requires that the number of the tetrahedra, N , is larger than some minimum depending upon the degree of the spline, so theoretically, we can determine the fineness of a tetrahedral partition by setting N ≈ c 1 n γ log(n) + c 2 , in which · denotes the integer part, c 1 , c 2 are tuning constants, and 1/( + 2) < γ < 1. In practice, we can implement the proposed method in a coarse-to-fine resolution manner with similar criteria discussed in Section 5.3. Our extensive Monte Carlo simulation studies suggest that once a tetrahedral partition is fine enough, further partition refinement usually has little effect on the fitting accuracy. To avoid model over-fitting, we can stop the partition refinement when the model performance on the test set sits flat or even worsens.

Numerical studies
In this section, we conduct various simulation studies to assess the performance of the proposed TPST method. We use studies in Sections 6.1 and 6.2 to illustrate the capabilities of TPST in handling a variety of complex data structures, and we use Section 6.3 to show the superiority of TPST in data reduction as well as signal enhancement, compared to other smoothing methods. More specifically, depending on the nature of the point clouds, we conduct experiments with unstructured (random design) and structured (fixed design) point clouds. The case of random design in Section 6.1 mimics the scenario that the point clouds are collections of 3D points distributed randomly in space. In contrast, the case of structured design is studied in Section 6.2. In this structured design setting, the point cloud is reconstructed from grid data; in other words, the locations of the points are deterministic. In Section 6.3, we explore the practical performance of TPST on a structured point cloud by using PET scan data as a reference.

Unstructured complete point clouds
In this example, we consider an unstructured (random design) point cloud in which the observations are randomly generated over the entire domain. We set the number of observations in each point cloud as n = 20, 000 and 50, 000. To mimic some complicated scenarios in practice, we generate the point clouds from the following two domains: (i) a cuboid with a hole inside (Ω 1 ); (ii) a 3D horseshoe (Ω 2 ). These domains are illustrated in Figure 6.1. To extract the underlying signal from point clouds, we consider the model in (3.1), where the random noises, ε i 's, are assumed to be independent and identically distributed and follow a normal distribution N (0, σ 2 ). For each domain, we consider two types of underlying functions with different degrees of variation: m 1 and m 2 for Ω 1 (Figure 6.1 (c) and (e)), and m 3 and m 4 for Ω 2 (Figure 6.1 (i) and (k)). The noises level, σ, is chosen according to the peak signal-to-noise ratios (PSNR) defined as In this study, we set PSNR = 5 and 10, representing scenarios of high and moderate noise levels, respectively. We use this study to investigate the effect of the size of the point cloud, degree of spline polynomial, and tetrahedral partition. We also compare the proposed method with the traditional tensor product spline method (Stone, 1994) and thin plate spline smoothing (Wood, 2003). We implement the tensor product spline and the thin plate spline smoothing using the function gam in R package mgcv. For a fair comparison, we set the dimension of the basis used to represent the underlying function for the thin plate spline and tensor product spline to be very similar to the dimension we used for the TPST method. As a result, the dimension after the numerosity reduction is comparable among different methods.
To evaluate the estimation and prediction performance of each method, we calculate the out-of-sample mean integrated squared error (MISE). To evaluate the effect of tetrahedral partitions on the TPST, we consider a relatively coarse mesh 11 for Ω 1 and 21 for Ω 2 , and a fine mesh 12 for Ω 1 and 22 for Ω 2 . An illustration of these partitions is given in the Figure 6.1 (a), (b), (g), and (h), and the number of vertices and the number of tetrahedra are summarized in Table 6.1. From Figure 6.2, we can see that for each domain, the TPST estimator has a similar performance based on the two partitions. As discussed in Section 5.4, for TPST, when the number of the tetrahedra is sufficiently large to capture the pattern and features of the underlying function, more delicate tetrahedral partitions will not benefit the estimation or the data reduction. For example, for domain Ω 2 , when we keep d = 3, the remaining dimensions after the numerosity reduction are 323 and 340 for the two different tetrahedral partitions, respectively.
Furthermore, based on Figure 6.2, one can see that the difference between estimates with d = 3 and d = 4 are relatively small. When we use a larger d, the estimators are usually less    (120) 456 (180) 504 (207) 1236 (426) biased but with larger variance and more computationally intensive. In this example, d = 3 is preferred as they are more efficient in data reduction. In general, the choice of d depends on the smoothness of the underlying function, the strength of signals, and computing resources.
As discussed previously, the spline methods are computationally efficient since they provide a global estimator. In this simulation example, the TPST method takes less than twenty seconds to fit the model for most of the simulation samples on a single Intel E5-2640 v3 core.
The speed is comparable to tensor product and thin plate spline estimators.

Structured point clouds with missing data
One of the critical aspects of analyzing point clouds data is handling the missing values.
We use this example to investigate further the proposed method on a structured point cloud at different resolutions with different missing schemes. To be more specific, we mimic three types of missing schemes, including (i) complete data with no missing, (ii) missing at random, and (iii) missing in a contiguous block as well as at random. Figures 6.3 illustrate different types of missing data we deal with in this example. Furthermore, we explore various missing rates under different missing mechanisms. For missing at random, we consider the missing rates ranging uniformly from 0 to 0.5, where 0 represents no missing voxels, and 0.5 means half of the voxels are missing. The contiguous block shown in Figure 6.3 contains 12% of the data.

Biomedical imaging analysis
In our increasingly aging societies, Alzheimer's disease (AD) has become the most frequent cause of dementia. Much progress has been made in assisting the early diagnosis of AD with neuroimaging techniques. One widely used neuroimaging technique is PET imaging, which can also be considered as an example of structured point clouds. However, the traditional PET scanning technique limits the overall resolution of the brain image, and there is a lack of effective and efficient image reconstruction methods. In this example, we apply the proposed TPST method to denoise and enhance the resolution of an actual PET image while improving the storage efficiency. We consider the PET image from the longitudinal Alzheimers Disease Neuroimaging Initiative (ADNI) study (http://adni.loni.usc.edu). The illustration of the transverse, coronal, and sagittal planes is shown in Figure 6.5 (a).
As discussed in Section 6.1, the partition only has a very small effect on the performance of the TPST estimator as long as the number of the tetrahedra is sufficient to capture the underlying pattern. Therefore, when implementing the TPST method for brain imaging analysis, we consider the tetrahedral partition shown in Figure 6.5 (b) with 317 tetrahedra and 117 vertices. To capture the high-frequency oscillations in the human brain, we set d ≥ 5, which is higher than the degree used in the previous two simulation examples.
To compare the estimation accuracy, we consider four different measures: (i) root mean squared error, iii) peak signal-to-noise ratio, PSNR = 20 log 10 {max(W i )/RMSE}, and (iv) mean relative absolute error, MRAE = n −1 n i=1 | W i − W i |/W i . The estimation and numerosity reduction results of different methods are shown in Table 6.2. From this table, one sees that the TPST, Thin Plate smoothing, and Tensor Product spline provide similar RMSE, MAE and PSNR when the target percentage of the data reduction is similar. However, the MRAEs for the Thin Plate smoothing and Tensor Product spline are significantly larger than those for TPST. For brain image, although it is collected on rectangular parallelepiped grid points, the accurate signals are only contained in the voxels within the human brain. As a result, the Thin Plate and Tensor Product methods are highly affected by the noises outside the actual brain. On the other hand, the TPST method provides a more accurate estimation for deblurring and enhancing the underlying signal for the original image. When d = 5, the TPST method achieves a similar dimension from data reduction: from a size of 510,340 voxels to a vector of spline coefficients of the length around 2000. If a higher degree d is affordable, we can achieve an even more accurate estimate but with a larger size of basis expansion (or input dimension).
In practice, the selection of d is the optimal balance between the percentage of data reduction and visual quality.  Figure 6.5: An illustration of an actual PET image and the tetrahedral partitions of the brain domain.

Conclusions and discussion
Challenges in handling and analyzing the irregularly-shaped 3D point cloud data motivated the advanced statistical methods to uncover the underlying trajectory over the point clouds.
Unlike conventional smoothing methods, the proposed TPST methods are proved to perform well for complex objects. Moreover, the proposed approaches are able to handle the irregularly and regularly missing data problem, effectively denoise or deblur the data while preserving inherent geometric features or spatial structures, and providing multi-resolution reconstruction.
The experimental results demonstrate the effectiveness of the proposed approaches compared to existing smoothing techniques.
Modern geoinformation technologies have provided a variety of feasible means for generating point clouds with billions of points. The enormous size of point clouds data always leads to serious time and memory consumption problems in data processing, storage, visualization, and transmission. Two approaches are currently considered as the most promising ones: parallel computing and data reduction. As a numerosity reduction method, the proposed method is an efficient data reduction technique which replaces the original big cloud point by a much smaller form of data representation. On the hand other, parallel computing is very appealing to address such "big" data issues. Unlike tensor product spline or thin-plate spline, one unique feature of the proposed TPST is its great scalability in computing. Specifically, the spline basis function is generated restricted to each tetrahedron without any overlap, and the smoothness is achieved only by the constraints on the spline coefficients. Based on this feature, we plan to develop a parallel computing algorithm based on domain decomposition in future works.

Supplemental Materials for "Nonparametric Regression for 3D Point Cloud Learning"
A. Detailed introduction of trivariate splines on tetrahedra In this section, we give the detailed discussion of the trivariate splines on tetrahedra, along with the introduction for some theoretical properties. In Section A.1, we first introduce the barycentric coordinates associated with a tetrahedron and show that a trivariate polynomial can be written in a convenient form using the barycentric coordinates. In Section A.2, we describe the directional derivatives of a polynomial, and smoothness conditions for polynomials on adjoining triangular face. The introduction of trivariate splines on a tetrahedra partition is given in Section 2.3 in main part.

A.1. Barycentric coordinates and Bernstein basis polynomials
Given a tetrahedron T = v 1 , v 2 , v 3 , v 4 , any fixed point p := (x, y, z) ∈ R 3 has a unique representation in terms of The barycentric coordinates (b 1 , b 2 , b 3 , b 4 ) also have an interesting geometric interpretation. As shown in Figure A.1, for any p ∈ T , it divides the tetrahedron T to four sub-tetrahedra, T 1 , T 2 , T 3 and T 4 . Notice that det(M) = 6V T , then b ı = V Tı /V T , where T ı replace the vertex v ı in V T with p, ı = 1, . . . , 4.
For a nondegenerate tetrahedron T and a point p ∈ T with barycentric coordinates (b 1 , b 2 , b 3 , b 4 ), for nonnegative integers i, j, k, l with i + j + k + l = d, define trivariate Bernstein basis polynomial of degree d relative to T as For any positive integer d and tetrahedron T , let P d (T ) be the space of all trivariate polynomials defined on T with degrees less than or equal to d. Note that the dimension of P d (T ) is d+3 3 . According to Theorem 15.8 in Lai and Schumaker (2007) and Lemma A.2 in Appendix Section A.4.1, the set of Bernstein basis polynomials 2. for all p ∈ T , 0 ≤ B d,T ijkl (p) ≤ 1; 3. B d,T ijkl has a unique maximum at the point d −1 (iv 1 + jv 2 + kv 3 + lv 4 ). Remark A.1. Barycentric coordinates are invariant to linear transformations of Cartesian coordinates, that is, they do not depend on the orientation or location of the tetrahedra. Consequently, the trivariate splines based on Bernstein basis polynomials, which are constructed with barycentric coordinates, are also invariant to linear transformations. where the coefficients γ T = {γ T ;ijkl } i+j+k+l=d are called B-coefficients of φ. For the purpose of computer implementation, in this paper, we employ the lexicographical order for ordering of the coefficients γ T . To be specific, γ T ;ijkl orders ahead of γ T ;i j k l either (i) i > i , or (ii) i = i and j > j , or (iii) i = i , j = j and k > k , or (iv) i = i , j = j , k = k and l > l . Consequently, we can express Accordingly, in (A.2), the index of the element γ T ;ijkl in the vector γ T is: Note that using a different ordering method will not affect the evaluation results for the trivariate polynomial functions. It is convenient to derive conditions of continuous connection for polynomials defined on adjacent tetrahedra in using barycentric coordinates and Bernstein basis polynomials. We use the following Example A.1 for illustration. Assume two sets of Bernstein polynomial basis {B d,T ijkl (p)} i+j+k+l=d and { B d, T ijkl ( p)} i+j+k+l=d are defined on T and T using the barycentric coordinates, respectively. Consider two degree-d polynomials φ(p) and φ( p) defined on T and T , respectively, with B-forms For point p F ∈ F , the barycentric coordinates with respect to T and T are (0 Therefore, φ and φ are continuous on F if and only if for j, k, l ≥ 0 and j + k + l = d.

A.2. Directional derivatives and smoothness
Following the definitions of directional derivatives in Section 5.1, we illustrate the theoretical properties of directional derivatives as follows.
Lemma A.1. (Lemma 15.12 in Lai and Schumaker, 2007) Consider direction u with directional coordinates a = (a 1 , a 2 , a 3 , a 4 ). Then for any p ∈ T and i + j + k + l = d.
Consequently, one can obtain the directional derivative for any trivariate polynomial φ.
Theorem A.1. (Theorems 15.13 and 15.14 in Lai and Schumaker, 2007) Consider direction u with directional coordinates a = (a 1 , a 2 , a 3 , a 4 ). Then for any trivariate polynomial φ with B-form in (A.1), the directional derivative is ijkl (a) = a 1 γ i+1,j,k,l + a 2 γ i,j+1,k,l + a 3 γ i,j,k+1,l + a 4 γ i,j,k,l+1 . In general, given u 1 , . . . , u m with associated directional coordinates a (ı) = (a where the coefficients are defined recursively as follows: In preparation for the discussion of trivariate spline and spline spaces, we need the following conditions for a smooth join between two polynomial on adjoining tetrahedra, like Example A.1 illustrated in Figure A Note that the barycentric coordinate of v 5 with respect to T is (2, −1, 0, 0).
If the trivariate polynomial is continuous over the whole region, then applying (A.8) for m = r = 0 generates which matches the conclusion in (A.4). Based on the de Casteljau Algorithm (Theorem A.3 in Appendix Section A.4.2), we propose the following computationally efficient algorithm to calculate the derivatives of φ(p).
Algorithm 1: Algorithm for the derivatives of φ(p).

A.3. Example of the constraint matrix
We illustrate the construction of the constraint matrix H using the Example A.1 presented in Figure A.3.
Continue our discussions in Section A.2 and consider the case of d = 2 piecewise polynomial. Recall that there are in total 10 Bernstein basis polynomials, with coefficients {γ ijkl } and { γ ijkl } in (A.9), and the barycentric coordinate of v 5 with respect to T is (2, −1, 0, 0).
As shown in (A.10), if the trivariate polynomial is continuous over the whole region, then γ 0ijk = γ 0ikj for i + j + k = d. In this case, we can write the constraint matrix H as Furthermore, if the trivariate polynomial has continuous first derivatives over the whole region, then (A.8) holds for both m = r = 0 and m = 1. Thus, in addition, we also need for any non-negative integers i, j, k such that i + j + k = 1, Therefore, for r = 1 where the trivariate spline has continuous first derivatives over the whole region, we will obtain the following H Here, the first six rows match the form in (A.11), and the last three rows correspond to (A.12).  = (b 1 , b 2 , b 3 , b 4 ) are the barycentric coordinates of p. Then the barycentric coordinates of p + tu are (b 1 + ta 1 , b 2 + ta 2 , b 3 + ta 3 , b 4 + ta 4 ). Thus, for i + j + k + l = d,
To show Theorem A.2, we need the following results.
Proof. Define operators E 1 γ ijkl = γ i+1,j,k,l , E 2 γ ijkl = γ i,j+1,k,l , E 3 γ ijkl = γ i,j,k+1,l and E 4 γ ijkl = γ i,j,k,l+1 . Thus, by (A.13), for i . Then for any 1 ≤ n ≤ d, the n-th order directional derivative of φ with respect to the direction u = v 4 − v 2 is given by Proof. The directional coordinates of u = v 4 − v 2 is a = (−1, 0, 0, 1). By (A.14), Proof of Theorem A.2. The equivalence between Statements 1 and 2 is obvious by definition. Thus, we just show the equivalence between Statements 2 and 3. i) We start to consider when r = 0. It is equivalent to consider directions u along F , where the directional coordinates are (0, 1 with respect to the tetrahedron T and T , respectively. Thus, φ and φ join continuously along F if and only if That is, γ 0jkl = γ 0jlk , with j + k + l = d, which matches our conclusion in (A.8).
By Theorem A.1, for each 0 ≤ n ≤ r, Since for points p F ∈ F , B d−n 0jkl (p F ) = B d−n 0jlk (p F ), it follows that (A.15) holds if and only if for j + k + l = d − n, n = 0, . . . , r, γ (n) 0jkl ( a) = γ In another direction, following the proof of Lemma A.3, for j + l + k = d − n,

A.5. Example of directional derivatives for basis functions and penalty matrix
In the following, we provide a simple example to illustrate how to calculate the first and second order derivatives and therefore to construct the penalty matrix P T in (5.1) in Section 5.2. Since the penalty function is calculated through the second-order derivatives in terms of the three-dimensional Cartesian coordinate system x = (1, 0, 0), y = (0, 1, 0) and z = (0, 0, 1), in the following example, we only consider directional derivatives for x, y, z.

B.1. Preliminaries
In this section, we start with the stability condition for trivariate spline bases over tetrahedralizations, and establish the uniform rate of the approximation of empirical inner product to the theoretical inner product based on the stability property.
Lemma B.1. Let {B ξ } ξ∈M be the basis for S r d ( ) constructed in Lai and Schumaker (2007), where M stands for the index set of spline bases. Then there exist positive constants C 1 , C 2 depending on degree d and partition quasi-uniform parameter β such that Proof. The proof of Lemma B.1 follows directly from the Theorem 17.18 in Lai and Schumaker (2007).
We illustrate the connection between the theoretical norm and L 2 norm in Lemma B.2.
Proof. It is straightforward to obtain upper and lower bounds of ξ∈M c ξ B ξ 2 Ω by Assumption (A3).
Therefore, Proposition 1 is established.
Next we give the proof of Proposition 2.
Let Γ ρn be the symmetric positive definite matrix We firstly show the bounds of Γ ρn .
In the other side, as shown in the proof of Proposition 1 and by Lemma B.1, for some positive constant C ρ .
(ii) By (i), Γ ρn is an invertible symmetric matrix, and its condition number c d = ρ max (Γ ρn )/ρ min (Γ ρn ) satisfies 1 < c d ≤ c −1 ρ C ρ . According to the definition of banded matrix in DeVore and Lorentz (1993), a matrix A = (a ij ) is said banded with bandwidth b if a ij = 0, |i − j| ≥ b, and if b is the smallest integer with this property. Based on the construction of trivariate splines, Γ ρn is a banded matrix with bandwidth b = d+3 3 . By Theorem 13.4.3 in DeVore and Lorentz (1993), Lemma B.4. Under Assumptions (A2) and (A4) Proof. Note that s 0,ε (p) = ξ∈M c 0,ξ B ξ (p) for some coefficients c 0,ξ , so the order of s ρn,ε (p) is related to that of c 0,ξ . In fact where c 0 = ( c 0,ξ ) ξ∈M and Γ 0 is the symmetric positive definite matrix Next we show that with probability 1 To prove (B.15), we decompose the noise variable ε i into a truncated part and a tail part It is straightforward to verify that the mean of the truncated part is uniformly bounded by D −2 n , so the boundedness of trivariate spline basis and of the function σ 2 entail that Next we show that tail part vanishes almost surely. Recall that E|ε 2+η n | ≤ υ η in Assumption (A2), so (B.17) By the Borel-Cantelli Lemma, we have P {ω|∃N (ω) , |ε n (ω)| ≤ D n for n > N (ω)} = 1.
Proof of Proposition 3. Note that the penalized spline s ρn,ε of ε is characterized by the orthogonality relations n ε − s ρn,ε , u n,Ω = ρ n s ρn,ε , u E , for all u ∈ S r d ( ).
Note that B(p) 2 = ξ∈M B 2 ξ (p) is bounded above by a constant C 2 < ∞ and below by C 1 > 0 for any point p ∈ Ω. Indeed, if C 1 = 0 for a point p, then B ξ (p) = 0 for all ξ ∈ M. That is, we have 0 = ξ∈M B ξ (p). It follows that these basis functions are linearly dependent which is a contradiction. On the other hand, for any fixed point p ∈ Ω, there is at most C 2 nonzero terms in the above summation, where C 2 is dependent on the smallest angle of the tetrahedralization . Note that B ξ is uniformly bounded for all ξ ∈ M, say bounded by 1. Thus, we know C 2 < ∞. We summarize the above discussion to get n −1 C 1 c 2 σ α max (ρ n ) −2 α min (0) ≤ Var {s ρn,ε (p)|P} ≤ n −1 C 2 C 2 σ α min (ρ n ) −2 α max (0) with probability approaching one. We now spend some effort to estimate the largest and smallest eigenvalues of Γ ρn . It is easy to see that for any vector a = [a ξ , ξ ∈ M] , a Γ ρn a = a 1 n Let s a = ξ∈M a ξ B ξ ∈ S r d ( ) be the spline associated with vector a = (a ξ , ξ ∈ M) T . By (B.4), we have s a 2 n,Ω = 1 + O P log n n/N s a 2 L 2 (Ω) ≤ K 2 1 + O P log n n/N | | 3 a 2 .
Here we have used the stability conditions in Lemma B.1. Furthermore, using Markov's inequality, we have ρ n n E ξ∈M a ξ B ξ ≤ ρ n n C | | 4 s a 2 L 2 (Ω) ≤ ρ n n C | | K 2 a 2 .
Thus, the largest eigenvalue α max (ρ n ) of the matrix Γ ρn in (B.13) is less than or equal to K 2 1 + O P log n n/N | | 3 + K 2 ρ n n C | | .
This establishes the result in the case for ρ n = 0. Similar for the case ρ n > 0. We have therefore completed the proof.
The above variance result can be more precise when spline space S −1 0 ( ) is considered. The next lemma provides the pointwise variance of s ρn,ε when using spline space S −1 0 ( ).
Lemma B.5. Consider piecewise constant spline space S −1 0 ( ) and suppose that the tetrahedra are of equal size. Under Assumptions (A1)-(A4), if the density function f (p) of p is continuous and positive on Ω, then Var {s ρn,ε (p)} = 1 n where V T is the volume of the tetrahedron T as defined before.
Hence, finding the asymptotic variance of s ρn,ε (p) is equivalent to finding the asymptotic variance ofs ρn,ε (p). Next we calculate the pointwise variance ofs ρn,ε (p). Note that Var {s ρn,ε (p)} = B(p) E c ρn,εc ρn,ε B(p) For any p ∈ Ω, let T ξ(p) be the tetrahedron that contains p. It is easy to see that f (p) dp.
For any continuous function g, let ω(g, ) = sup p,p ∈Ω, p−p ≤ |g(p) − g(p )| be the moduli of continuity of g on Ω. Then for any p ∈ T ξ(p) , we can write {f ( L 2 (Ω) B 2 ξ (p) The conclusion follows.