Spline Smoothing of 3D Geometric Data

Xinyi Lia, Shan Yub, Yueying Wangc, Guannan Wangd and Li Wangc †† _Address for
correspondence_ : Li Wang, Department of Statistics and the Statistical
Laboratory, Iowa State University, Ames, IA, USA. Email<EMAIL_ADDRESS>
aClemson University, bUniversity of Virginia, cIowa State University and
dCollege of William & Mary

> Abstract: Over the past two decades, we have seen an increased demand for 3D
> visualization and simulation software in medicine, architectural design,
> engineering, and many other areas, which have boosted the investigation of
> geometric data analysis and raised the demand for further advancement in
> statistical analytic approaches. In this paper, we propose a class of spline
> smoothers appropriate for approximating geometric data over 3D complex
> domains, which can be represented in terms of a linear combination of spline
> basis functions with some smoothness constraints. We start with introducing
> the tetrahedral partitions, Barycentric coordinates, Bernstein basis
> polynomials, and trivariate spline on tetrahedra. Then, we propose a
> penalized spline smoothing method for identifying the underlying signal in a
> complex 3D domain from potential noisy observations. Simulation studies are
> conducted to compare the proposed method with traditional smoothing methods
> on 3D complex domains.
>
> Key words and phrases: Complex domain; Nonparametric smoothing; penalized
> splines; Trivariate splines; Tetrahedra partitions.

1\. Introduction

Recent advances in computer and information technology have dramatically
boosted the availability of three-dimensional (3D) geometric data in many
fields, as illustrated in Figure 1.1, including manufacturing ((a)–(c)),
medical science ((d)–(e)), arts ((f)–(g)) and animation ((h)–(i)). For
example, by utilizing 3D technology within medical imaging, we are able to
take and combine those two-dimensional (2D) cross-section slices into a
concise 3D visual of the area being scanned. With 3D medical imaging,
healthcare professionals can now access new angles, resolutions, and details
that offer an all-around better understanding of the body part in question. As
a result, statistical learning for 3D geometry is emerging as a critical
research area in many fields.

|  |  |  |
---|---|---|---|---
(a) | (b) | (c) | (d) | (e)
|  |  |  |
(f) | (g) | (h) | (i) |
Figure 1.1: Examples of 3D geometric data and corresponding tetrahedral meshes
generated with iso2mesh toobox and distmeshnd function using MATLAB.

In reality, 3D geometric data usually contains unwanted noise that obscures
the features of the volumes of interest. Consequently, denoising is critical
in the 3D geometric analysis; for example, in many imaging studies, smoothing
methods are widely adopted in image processing tasks (Goldsmith et al., 2014;
Kang, 2020). As shown in Figure 1.1, the objects (or regions) of interest are
usually irregularly shaped. Many conventional denoising 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” across the complex domains, which refers to the poor
estimation over difficult regions caused by smoothing inappropriately across
boundary features (Ramsay, 2002; Wang and Ranalli, 2007; Kim and Wang, 2021).
Several techniques have been developed in the last decades to tackle the
“leakage” problem when recovering the true signals. Many of the work focus on
2D manifold domains, such as the spline smoothing (Wang and Ranalli, 2007;
Wilhelm et al., 2016), finite element analysis (Lindgren et al., 2011), and
Kernel smoothing (Guo et al., 2010). For 3D cases, Chung et al. (2018) applied
the discrete heat kernel smoothing to irregular image domains, and Huo et al.
(2020) considered the finite element method for irregular 3D subjects.

To deal with this prevalent challenge in 3D geometric data analysis, we
propose to use trivariate penalized spline smoothing methods based on
tetrahedral partitions. To handle the complex 3D domains, it is crucial to
construct an ideal representation of the 3D geometry. In computer graphics and
geometric-aided design, tetrahedral partitions are the standard representation
of 3D objects. For example, many 3D scanners use the laser to capture an
object’s exact size and shape into the computer as a digital 3D triangle mesh
representation. Although traditional medical planning systems often use either
voxel-based or surface-based models, tetrahedral meshes are more flexible than
both data structures because complex volumetric data sets can be represented
with few primitives. After an appropriate partition of the domain, we can
construct the trivariate spline basis function on the tetrahedral partition to
approximate data distributed over 3D complex domains.

Moreover, the characteristics of 3D geometry data create substantial
challenges in analysis (Wong et al., 2016). First, we need to take into
account the large scale and high resolution of the collected data. For
example, a standard PET scan from one subject contains ($79\times 95\times
68$) 510,340 voxels. Therefore, even for a single image, the number of
observations is very high. One of the significant advantages of the spline
smoothing method lies in its computational simplicity. As a global estimation
with an explicit model expression, the fitting of splines techniques is
achieved by solving a single linear system.

Besides the computation issues raised by the large volume of data, the
challenges also arise from the complicated structure of the signal. On one
hand, the underlying data may vary mildly or sharply according to the
location. We propose the Adaptive Trivariate Penalized Spline on Tetrahedral
Partitions (ATPST) by adopting adaptive weights in smoothing spline to
overcome this challenge. On the other hand, missing data is common and
inevitable in 3D geometry data acquisition. For example, in 3D seismology, the
incomplete and sparse data in the collection is often caused by physical and
economic constraints (Chai et al., 2020). In our simulation study, we mimic
the cases in real life and conduct experiments on missing data with multiple
missing schemes. The simulation results reveal the superiority of the proposed
estimators, regardless of the setting.

Other than denoising or deblurring, the reconstruction of 3D geometry data in
a different (usually higher) resolution is another important task in analysis.
In terms of statistics, this is equivalent to providing an accurate prediction
based on estimation. The proposed spline smoothing methods illustrate supreme
prediction power in simulation studies.

The rest of the paper is structured as follows. In Section 2, we give an
overview of the tetrahedral partitions of a 3D domain. Once a tetrahedral
partition has been constructed, we can construct the spline space over the
tetrahedral partition. In Section 3, we introduce the penalized spline
estimators. We then demonstrate the smoothers’ performance on simulated
problems in Section 4. Finally, this paper is concluded, and future work is
outlined in Section 5.

2\. Tetrahedral partitions and trivariate splines on tetrahedr

In this section, we provide a basic framework for tetrahedral partitions and
trivariate splines on those partitions. More detailed introductions can be
referred to online supplementary materials (Appendices A and B).

Let $Y_{i}$ be the response variable observed on the $i$th location or design
point $\boldsymbol{v}_{i}=(x_{i},y_{i},z_{i})$, $i=1,\ldots,n$. Suppose all
$\boldsymbol{v}_{i}$’s range over a bounded domain
$\Omega\subseteq\mathbb{R}^{3}$ of arbitrary shape, for example, the 3D
horseshoe domain shown in Figure 1.1 (a). For any given dataset with design
points located on $\Omega$, one can find a polygonal domain (a domain with a
piecewise linear boundary) to include all the design points. With this
principle, in the following, we assume that $\Omega$ is a polygonal domain
itself, which may have one or multiple holes not containing any design points.

2.1. Tetrahedral partitions

In the following, we use $T$ to denote a tetrahedron, which is a convex hull
of four noncoplanar points in $\mathbb{R}^{3}$. A collection
$\triangle=\\{T_{1},\ldots,T_{N}\\}$ of $N$ tetrahedra is called a tetrahedral
partition of $\Omega=\cup_{h=1}^{N}T_{h}$ provided that any pair of tetrahedra
in $\triangle$ intersect at most at a common vertex, along a common edge, or
along a common triangular face. In this paper, we restrict our attention to
tetrahedralizations of $\Omega$ because any polygonal domain of arbitrary
shape can be partitioned into finitely many tetrahedra. Without loss of
generality, we assume that all the location points $\boldsymbol{v}_{i}$’s are
inside tetrahedra of $\triangle$, that is, they are not on edges, vertices, or
faces of tetrahedra in $\triangle$. Otherwise, we can simply count them twice
or multiple times if any observation is located on a face, an edge, or at a
vertex of $\triangle$ (Lai and Schumaker, 2007).

2.2. Tetrahedral partitions in practice

We use Figure 2.1 to illustrate the flowchart of obtaining a tetrahedral
partition in practice, given a 3D object (taking the cow model in Figure 1.1
(h) as an example) with arbitrary complex shape. This mesh generation
procedure also works for other typical 3D data objects, like those in Figure
1.1.

Figure 2.1: A flowchart demonstrating how a tetrahedral partition can be
obtained from a 3D object with arbitrary complex shape.

For practical purposes, there are plenty of tetrahedral partitions that
satisfy our definition in Section 2.1, wherein the Delaunay tetrahedral
partition is most preferred. The Delaunay tetrahedral partition has been well-
studied with superior properties for 2D objects (Chew, 1989; Lee and
Schachter, 1980), for example, a 2D Delaunay triangulation is “angle optimal”,
i.e., it maximizes the minimum angle in a triangulation. Extending the 2D
Delaunay triangulation to higher dimensionality, a 3D Delaunay tetrahedral
partition possesses similar property that it maximizes the minimal solid angle
in the tetrahedra.

In the past few decades, Delaunay-based mesh generation algorithms for 3D
domains have also been implemented in various software packages and toolboxes,
such as the functions delaunayTriangulaion, delaunay and distmesh/distmeshnd
(Persson and Strang, 2004) in MATLAB, QualMesh (Cheng et al., 2005) and the
CGAL library (Jamin et al., 2015; The CGAL Project, 2020) for C++, TetGen(Si,
2015) and iso2mesh (Fang and Boas, 2009) built on CGAL as well as TetGen for
MATLAB/Octave. Note that even though all the algorithms can be used to
generate triangular meshes, the implementation details and usages differ. For
example, the functions distmesh and distmeshnd take distance functions as
representations of boundary geometry, while delaunayTriangulation generates
Delaunay 2D/3D triangulation based on given point sets. In practice, when the
mathematical representations of boundary geometry or the vertex sets for
Delaunay triangulations are not immediately available, one can employ TetGen
instead, whose inputs are triangular surface mesh on domain boundaries, as
shown in Figure 2.1.

2.3. Trivariate splines on a tetrahedra partition

Given a tetrahedron
$T=\langle\boldsymbol{v}_{1},\boldsymbol{v}_{2},\boldsymbol{v}_{3},\boldsymbol{v}_{4}\rangle$,
any fixed point $\boldsymbol{v}:=(x,y,z)\in\mathbb{R}^{3}$ has a unique
representation in terms of
$\langle\boldsymbol{v}_{1},\boldsymbol{v}_{2},\boldsymbol{v}_{3},\boldsymbol{v}_{4}\rangle$,

$\boldsymbol{v}=b_{1}\boldsymbol{v}_{1}+b_{2}\boldsymbol{v}_{2}+b_{3}\boldsymbol{v}_{3}+b_{4}\boldsymbol{v}_{4},~{}~{}\text{
with }b_{1}+b_{2}+b_{3}+b_{4}=1,$

where $(b_{1},b_{2},b_{3},b_{4})$ are called the barycentric coordinates of
$\boldsymbol{v}$ relative to the tetrahedron $T$, and they are nonnegative if
$\boldsymbol{v}$ 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

$B_{ijkl}^{d,T}(\boldsymbol{v}):=\frac{d!}{i!j!k!l!}b_{1}^{i}b_{2}^{j}b_{3}^{k}b_{4}^{l},\text{
with }i+j+k+l=d.$

For any positive integer $d$ and tetrahedron $T$, let $\mathcal{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 $\mathcal{P}_{d}(T)$ is
$\binom{d+3}{3}$. According to Theorem 15.8 in Lai and Schumaker (2007) and
Lemma D.1 in Section D.1 in the supplement, the set of Bernstein basis
polynomials

$\boldsymbol{B}^{d}_{T}(\boldsymbol{v}):=\\{B_{ijkl}^{d,T}(\boldsymbol{v}):i,j,k,l\geq
0,i+j+k+l=d\\}$

forms a basis for the space of polynomials $\mathcal{P}_{d}(T)$. Thus, given
Bernstein basis functions $\\{B_{ijkl}^{d,T}(\boldsymbol{v})\\}_{i+j+k+l=d}$,
any polynomial $p(\boldsymbol{v})\in\mathcal{P}_{d}(T)$ can be written
uniquely as B-form,

$p(\boldsymbol{v})|_{T}=\sum_{i+j+k+l=d}\gamma_{T;ijkl}B_{ijkl}^{d,T}(\boldsymbol{v})=\boldsymbol{B}^{d}_{T}(\boldsymbol{v})^{\top}\boldsymbol{\gamma}_{T},$

where the coefficients
$\boldsymbol{\gamma}_{T}=\\{\gamma_{T;ijkl}\\}_{i+j+k+l=d}$ are called
B-coefficients of $p$.

Consider the whole domain $\Omega$ and let $\mathcal{C}^{r}(\Omega)$ be the
collection of all $r$-th continuously differentiable functions over $\Omega$
for some nonnegative integer $r$. For tetrahedral partition
$\triangle=\\{T_{1},\ldots,T_{N}\\}$, let
$\mathcal{S}_{d}^{r}(\triangle)=\\{s\in\mathcal{C}^{r}(\Omega):s|_{T}\in\mathcal{P}_{d}(T),T\in\triangle\\}$
be a spline space of degree $d$ and smoothness $r$ over $\triangle$, where
$s|_{T}$ is the polynomial piece of spline $s$ restricted on tetrahedron $T$.

######  Remark 1.

There are a variety of choices for $d$ and $r$. According to Lai (1989), space
$\mathcal{S}_{d}^{r}(\triangle)$ with tetrahedral partition $\triangle$ and
$d\geq 6r+3$, reach the full approximation power.

According to (B.1) in Section B.1 in the supplement, for any
$s\in\mathcal{S}_{d}^{r}(\triangle)$, there exists
$\boldsymbol{\gamma}_{T_{J}}=\\{\gamma_{T_{J};ijkl}\\}$ such that

$s|_{T_{J}}(\boldsymbol{v})=\sum_{i+j+k+l=d}\gamma_{T_{J};ijkl}B_{ijkl}^{d,T_{J}}(\boldsymbol{v})=\boldsymbol{B}^{d}_{T_{J}}(\boldsymbol{v})^{\top}\boldsymbol{\gamma}_{T_{J}},~{}~{}J=1,\ldots,N.$

Consequently, let
$\boldsymbol{B}_{d}=((\boldsymbol{B}^{d}_{T_{1}})^{\top},\ldots,(\boldsymbol{B}^{d}_{T_{N}})^{\top})^{\top}$
and
$\boldsymbol{\gamma}=(\boldsymbol{\gamma}_{T_{1}}^{\top},\ldots,\boldsymbol{\gamma}_{T_{N}}^{\top})^{\top}$,
one has
$s(\boldsymbol{v})=\boldsymbol{B}_{d}(\boldsymbol{v})^{\top}\boldsymbol{\gamma}$.
To satisfy the smoothness conditions in $s$, certain smoothness constraints
are imposed on the B-coefficients $\boldsymbol{\gamma}$. One can obtain a
smoothness constraint matrix $\mathbf{H}$ such that
$\mathbf{H}\boldsymbol{\gamma}=\boldsymbol{0}$ by repeatedly applying the
equation (B.8) over all shared triangular faces, and more details are
available in Section B.2 in the supplement.

3\. Penalized spline estimators

With all the preparations in the previous introduction, we apply the
trivariate spline over tetrahedral partition to estimate functions on a
topologically and/or geometrically complex 3D domain $\Omega$ in this section.
For observations $\\{(\boldsymbol{v}_{i},Y_{i})\\}_{i=1}^{n}$, we consider the
following nonparametric regression model

$Y_{i}=m(\boldsymbol{v}_{i})+\epsilon_{i},$

where $m$ is some smooth but unknown function, and $\epsilon_{i}$ is random
error 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 balances the bias and variance of the
estimator of $m$. We formulate the roughness penalty approach as the following
penalized least squares problem:

$\min_{m}\sum_{i=1}^{n}\left\\{Y_{i}-m(\boldsymbol{v}_{i})\right\\}^{2}+\lambda\mathcal{E}(m),$
(3.1)

where $\mathcal{E}(m)$ is a penalty function with details given in Section
3.1, and $\lambda$ is the roughness penalty parameter. A larger $\lambda$
leads to a less fluctuating $m$. If $\lambda$ goes to infinity, our estimator
shrinks to linear functions where the roughness penalty $\mathcal{E}(m)=0$. A
proper penalty parameter $\lambda$ balances the goodness of fit for the data
and the volatility of estimated functions. We propose to seek for the
minimizer of (3.1) in $\mathcal{S}_{d}^{r}(\triangle)$ and denote it as
$\widehat{m}_{\lambda}$, which is the Trivariate Penalized Spline over
Tetrahedral partition (TPST) estimator of $m$.

3.1. Tetrahedral partitions

In this section, we describe the details of implementing the TPST method. For
the penalty function in (3.1), we consider

$\mathcal{E}(m)=\sum_{i+j+k=2}{2\choose i}{2-i\choose
j}\int_{\Omega}\\{D_{x}^{i}D_{y}^{j}D_{z}^{k}m(\boldsymbol{v})\\}^{2}\mathrm{d}\boldsymbol{v}$
(3.2)

as the roughness penalty. The penalty $\mathcal{E}(m)$ is a commonly used
penalty, the second-order derivative to zero for surface smoothing, see Green
and Silverman (1994) for the two-dimensional case.

For a spline function $s\in\mathcal{S}_{d}^{r}(\triangle)$, the roughness
penalty can be written as the following:

$\displaystyle\mathcal{E}(s)=\sum_{T\in\triangle}\mathcal{E}(s_{T})$
$\displaystyle=\sum_{T\in\triangle}\sum_{i+j+k=2}{2\choose i}{2-i\choose
j}\int_{T}\\{D_{x}^{i}D_{y}^{j}D_{z}^{k}s_{T}(\boldsymbol{v})\\}^{2}\mathrm{d}\boldsymbol{v},$

where
$s_{T}(\boldsymbol{v})=\sum_{i+j+k+l=d}\gamma_{T;ijkl}B_{ijkl}^{T,d}(\boldsymbol{v})$
is a function restricted on $T$. Specifically, we have

$\displaystyle\mathcal{E}(s_{T})$
$\displaystyle=\sum_{g,g^{\prime}\in\\{x,y,z\\}}\int_{T}\bigg{\\{}\sum_{i+j+k+l=d}\gamma_{T;ijkl}D_{g}D_{g^{\prime}}B_{ijkl}^{T,d}(\boldsymbol{v})\bigg{\\}}^{2}\mathrm{d}\boldsymbol{v}$
$\displaystyle=\sum_{g,g^{\prime}\in\\{x,y,z\\}}\boldsymbol{\gamma}^{\top}_{T}\mathbf{P}_{T}^{gg^{\prime}}\boldsymbol{\gamma}_{T}=\boldsymbol{\gamma}^{\top}_{T}\mathbf{P}_{T}\boldsymbol{\gamma}_{T},$

where
$\boldsymbol{\gamma}_{T}=\left\\{\gamma_{T;ijkl}\right\\}_{i+j+k+l=d}^{\top}$,
and $\mathbf{P}_{T}^{gg^{\prime}}$, for $g,g^{\prime}=x,y,z$ are the
$\binom{d+3}{3}\times\binom{d+3}{3}$ matrices with entries being
$\int_{T}\left\\{D^{2}_{gg^{\prime}}B_{ijkl}^{T,d}(\boldsymbol{v})\right\\}\left\\{D^{2}_{gg^{\prime}}B_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}^{T,d}(\boldsymbol{v})\right\\}\mathrm{d}\boldsymbol{v}$,
respectively. See Appendix C for more details in calculations of penalty
matrices $\mathbf{P}_{T}^{gg^{\prime}}$, for $g,g^{\prime}\in\\{x,y,z\\}$ and
$T\in\triangle$.

######  Remark 2.

For a spline function $s\in\mathcal{S}^{r}_{d}(\triangle)$ with B-coefficients
$\boldsymbol{\gamma}$, its roughness penalty $\mathcal{E}(s)$ can be written
as some quadratic function of $\boldsymbol{\gamma}$, that is,
$\mathcal{E}(s)=\boldsymbol{\gamma}^{\top}\mathbf{P}\boldsymbol{\gamma}$,
where $\mathbf{P}=\textrm{diag}(\mathbf{P}_{T},T\in\triangle)$ is a block
diagonal matrix.

Accordingly, we have
$\mathcal{E}(\mathbf{B}_{d}^{\top}\boldsymbol{\gamma})=\boldsymbol{\gamma}^{\top}\mathbf{P}\boldsymbol{\gamma}$,
where $\mathbf{P}$ is the block diagonal penalty matrix. Thus, the objective
function can be written as

$\min_{\boldsymbol{\gamma}}\sum_{i=1}^{n}\left\\{Y_{i}-\mathbf{B}_{d}(\boldsymbol{v}_{i})^{\top}\boldsymbol{\gamma}\right\\}^{2}+\lambda\boldsymbol{\gamma}^{\top}\mathbf{P}\boldsymbol{\gamma},~{}\textrm{subject
to~{}}\mathbf{H}\boldsymbol{\gamma}=\mathbf{0}.$ (3.3)

One can use QR decomposition to get rid of the constraint. Specifically,
$\mathbf{H}^{\top}=\mathbf{Q}\mathbf{R}=(\mathbf{Q}_{1}~{}\mathbf{Q}_{2})\binom{\mathbf{R}_{1}}{\mathbf{R}_{2}}$,
where $\mathbf{Q}$ is an orthogonal matrix, $\mathbf{R}$ is an upper
triangular matrix, $\mathbf{R}_{1}$ is a full rank matrix with the same rank
as $\mathbf{H}$, and $\mathbf{R}_{2}$ is a matrix of zeros. Note that for any
vector $\boldsymbol{\gamma}$ satisfying
$\mathbf{H}\boldsymbol{\gamma}=\mathbf{0}$, there exists some
$\boldsymbol{\theta}$ such that
$\boldsymbol{\gamma}=\mathbf{Q}_{2}\boldsymbol{\theta}$. Also, for any
$\boldsymbol{\theta}$,
$\mathbf{H}(\mathbf{Q}_{2}\boldsymbol{\theta})=\mathbf{0}$ holds. Thus, (3.3)
is equivalent to a penalized regression problem without constraint:

$\sum_{i=1}^{n}\left\\{Y_{i}-\mathbf{B}_{d}(\boldsymbol{v}_{i})^{\top}\mathbf{Q}_{2}\boldsymbol{\theta}\right\\}^{2}+\lambda\boldsymbol{\theta}^{\top}\mathbf{Q}_{2}^{\top}\mathbf{P}\mathbf{Q}_{2}\boldsymbol{\theta}.$
(3.4)

Let $\mathbf{Y}=(Y_{1},\ldots,Y_{n})^{\top}$,
$\mathbf{B}_{d}=(\mathbf{B}_{d}(\boldsymbol{v}_{1}),\ldots,\mathbf{B}_{d}(\boldsymbol{v}_{n}))^{\top}$.
Then, the solution of (3.4) is given by
$\widehat{\boldsymbol{\theta}}_{\lambda}=(\mathbf{Q}^{\top}_{2}\mathbf{B}_{d}^{\top}\mathbf{B}_{d}\mathbf{Q}_{2}+n^{-1}\lambda\mathbf{Q}_{2}^{\top}\mathbf{P}\mathbf{Q}_{2})^{-1}\mathbf{Q}^{\top}_{2}\mathbf{B}_{d}^{\top}\mathbf{Y}$
and
$\widehat{\boldsymbol{\gamma}}_{\lambda}=\mathbf{Q}_{2}\widehat{\boldsymbol{\theta}}_{\lambda}$,
which yields the TPST estimator
$\widehat{m}_{\lambda}(\boldsymbol{v})=\mathbf{B}_{d}(\boldsymbol{v})^{\top}\widehat{\boldsymbol{\gamma}}_{\lambda}$.
Several methods can be used for choosing the penalty parameter, such as the
block cross-validation and generalized cross-validation. Detailed discussion
about these methods is given in Section C.4 in the supplement.

3.2. The ATPST estimator

To deal with the complicated variation on the complex 3D domain, we propose a
total variation based estimator, ATPST, by adjusting the smoothness penalty
part in the TPST estimator.

In (3.2), the roughness penalty function imposes unique penalty weights across
$\Omega$, which is not ideal for functions with various degrees of volatility
among different regions. To provide a precise estimation for functions with
various volatility, we consider an adaptive penalty function where the penalty
weights are assigned according to the local volatility. To be specific, we
update the adaptive energy function $\mathcal{E}$ as

$\mathcal{E}^{\omega}(s)=\sum_{T\in\triangle}\omega_{T}\sum_{i+j+k=2}{2\choose
i}{2-i\choose
j}\int_{T}\\{D_{x}^{i}D_{y}^{j}D_{z}^{k}s_{T}(\boldsymbol{v})\\}^{2}\mathrm{d}\boldsymbol{v},$

where $\omega_{T}$ are the weights measuring local volatility. Similar to
Remark 2 in Section 3.1, for $s\in\mathcal{S}^{r}_{d}(\triangle)$ with
B-coefficients $\boldsymbol{\gamma}$,
$\mathcal{E}^{\omega}(s)=\sum_{T\in\triangle}\omega_{T}\boldsymbol{\gamma}_{T}^{\top}\mathbf{P}_{T}\boldsymbol{\gamma}_{T}$.
Then, one can obtain the ATPST by minimizing

$\sum_{i=1}^{n}\left\\{Y_{i}-\mathbf{B}_{d}(\boldsymbol{v})^{\top}\boldsymbol{\gamma}\right\\}^{2}+\lambda\sum_{T\in\triangle}\omega_{T}\boldsymbol{\gamma}_{T}^{\top}\mathbf{P}_{T}\boldsymbol{\gamma}_{T},~{}\textrm{subject
to~{}}\mathbf{H}\boldsymbol{\gamma}=\mathbf{0}.$

We define $\omega_{T}=(C-C_{m}^{-1}V_{T}^{-1}\|m\|_{T,TV})^{\tau}$, where
$C_{m}=\max\\{V_{T}^{-1}\|m\|_{T,TV},T\in\triangle\\}$, $C$ and $\tau$ are
some positive constants, which can be chosen through cross validation, and
$\|m\|_{T,TV}$ is the total variation of the function $m$ within tetrahedron
$T$. Note that $\omega_{T}$ is a monotonically decreasing function of
$\|m\|_{T,TV}$, which generates high gradient for regions with large
variation. Specifically,

$\|m\|_{T,TV}=\sup\left\\{\int_{T}m(\boldsymbol{v})\textrm{div}f(\boldsymbol{v})\mathrm{d}\boldsymbol{v}:f\in
C_{c}^{1}(\Omega,\mathbb{R}^{3}),\|f\|_{\infty,\Omega}\leq 1\right\\},$

where $C_{c}^{1}(\Omega,\mathbb{R}^{3})$ is the set of continuously
differentiable vector functions of compact support contained in $\Omega$,
$\textrm{div}f(\boldsymbol{v})$ is the divergence of function $f$, and
$\|\cdot\|_{\infty,\Omega}$ is the essential supremum norm. If $m$ is a
differentiable function, one has
$\|m\|_{T,TV}=\int_{T}[\\{D_{x}m(\boldsymbol{v})\\}^{2}+\\{D_{y}m(\boldsymbol{v})\\}^{2}+\\{D_{z}m(\boldsymbol{v})\\}^{2}]^{1/2}\mathrm{d}\boldsymbol{v}$.

In practice, we suggest using the TPST estimator as an initial estimator for
the total variation at tetrahedron $T$. Figure 3.1 shows the plot of the
estimated $\omega_{T}$ based on a typical data sample in our simulation study,
from which one sees that larger weights are assigned to the less fluctuant
region. By the property of the first-order partial derivative of spline
functions, $\|m\|_{T,TV}$ can be approximated by
$\int_{T}[\\{\widehat{\boldsymbol{\gamma}}_{T}^{\top}\mathbf{C}^{(1)}_{d}(x)\mathbf{B}_{d-1}(\boldsymbol{v})\\}^{2}+\\{\widehat{\boldsymbol{\gamma}}_{T}^{\top}\mathbf{C}^{(1)}_{d}(y)\mathbf{B}_{d-1}(\boldsymbol{v})\\}^{2}+\\{\widehat{\boldsymbol{\gamma}}_{T}^{\top}\mathbf{C}^{(1)}_{d}(z)\mathbf{B}_{d-1}(\boldsymbol{v})\\}^{2}]^{1/2}\mathrm{d}\boldsymbol{v}$,
where $\widehat{\boldsymbol{\gamma}}_{T}$ is the spline coefficient of the
TPST estimator at tetrahedron $T$.

Figure 3.1: Estimated $\omega_{T}$ under a random design.

4\. Simulation studies

In this section, we conduct various simulation studies to assess the
performance of the proposed methods, TPST and ATPST. For the ATPST, we
consider the weight functions defined in Section 3.2 with $\tau=2$ and $C$
being selected from a set of grid points ranging from 1 to 3 based on the GCV
criteria. We compare the proposed methods with the traditional tensor product
spline method (Stone, 1994) and thin plate spline smoothing (Wood, 2003). We
then implement the tensor product spline and the thin plate spline smoothing
using the function gam in R package mgcv.

In all simulation examples in this section, the data are generated from the
following model: $Y_{i}=m(\boldsymbol{v}_{i})+\epsilon_{i}$,
$\boldsymbol{v}_{i}\in\Omega$, $i=1,\ldots,n$. The random noises,
$\epsilon_{i}$’s, are independent and identically distributed and generated
from $N(0,\sigma^{2})$, where $\sigma$ is chosen according to the peak signal-
to-noise ratios (PSNR) defined as

$\textrm{PSNR}=10\log_{10}\left[\left\\{\max_{i}m(\boldsymbol{v}_{i})\right\\}^{2}/\sigma^{2}\right].$

In the simulation studies, we set $\mathrm{PSNR}=5$ and $10$, representing
scenarios of high and moderate noise levels, respectively. Depending on the
nature of the design points, we conduct numerical experiments with both random
and fixed designs. In a fixed design, the locations $\boldsymbol{v}_{i}$’s are
deterministic, and this is commonly used in imaging analysis; on the other
hand, the case of random design is more towards the statistical learning
setup.

4.1. Random design with complete data

In this example, we consider a random design in which the observations are
randomly generated over the entire domain. To be more specific, we set the
sample size $n=20,000$ and $50,000$. To mimic some complicated scenarios in
practice, we consider two different domains: (i) a cube with a hole inside
($\Omega_{1}$); (ii) a 3D horseshoe ($\Omega_{2}$). These domains are
illustrated in Figure 4.1. To compare the TPST and the ATPST, we consider two
types of functions with different degrees of variation for each domain. The
true functions shown in the first column of Figure 4.1 are smooth functions
with less volatility, while the functions shown in the second column are more
waggling. Furthermore, in this simulation study, we investigate the effect of
sample size, degree of spline polynomial, and tetrahedral partition.

|  |  |
---|---|---|---
(a) True $m_{1}$ | (b) True $m_{2}$ | (c) $\triangle_{1}$ for $\Omega_{1}$ | (d) $\triangle_{2}$ for $\Omega_{1}$
|  |  |
(e) True $m_{3}$ | (f) True $m_{4}$ | (g) $\triangle_{1}$ for $\Omega_{2}$ | (h) $\triangle_{2}$ for $\Omega_{2}$
Figure 4.1: True functions and tetrahedral partitions over domains
$\Omega_{1}$ and $\Omega_{2}$.

To evaluate the estimation and prediction performance of each method, we
calculate the out-of-sample mean integrated squared error (MISE). Figure 4.2
presents the average of the MISEs over 200 replications for all different
scenarios. Based on Figure 4.2, one can observe that as the sample size or the
PSNR increases, the estimation and prediction accuracy improves for all the
methods. Regardless of the simulation scenarios, the proposed TPST and ATPST
methods outperform the other two traditional methods since the traditional
methods cannot handle the “leakage” problem over the complex domain.

For the ATPST, as discussed in Section 3.2, we adaptively assign weights for
each penalty term based on the total variation within each tetrahedron $T$;
see Figure 3.1. Comparing the TPST and the ATPST, one can see that when the
actual underlying function is relatively smooth, the performance of the TPST
and the ATPST are similar to each other. However, when the actual underlying
function is fluctuant, the ATPST improves the estimation and prediction
accuracy of the TPST for most of the simulation scenarios.

---
Figure 4.2: Plots of average MISEs based on different methods under random
design.

To evaluate the effect of tetrahedral partitions on the TPST and ATPST, we
consider a relatively coarse mesh $\triangle_{1}$ and a fine mesh
$\triangle_{2}$ for each domain. An illustration of these partitions is given
in the third and the fourth columns of Figure 4.1, and the number of vertices
and the number of tetrahedrons are summarized in Table 4.1. From Figure 4.2,
we can see that for each domain, both TPST and ATPST estimators have similar
performance based on the two partitions. Moreover, for the proposed penalized
splines, when the number of the tetrahedron is sufficiently large to capture
the pattern and features of the data, finer tetrahedral partitions will not
further benefit the estimation.

Table 4.1: Number of tetrahedrons (vertices) in each partition.  $\Omega_{1}-\triangle_{1}$ | $\Omega_{1}-\triangle_{2}$ | $\Omega_{2}-\triangle_{1}$ | $\Omega_{2}-\triangle_{2}$
---|---|---|---
240 (120) | 456 (180) | 504 (207) | 1236 (426)

To study the effect of the degree of spline on the proposed TPST and ATPST.
Based on Figure 4.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 biased but with larger variance and more computationally
intensive. In this example, $d=3$ is preferred as they are more
computationally efficient. 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. The speed is comparable to tensor product and thin plate
spline estimators. The ATPST needs a little more time to select the best
weight functions, and it usually takes less than five minutes to obtain the
estimates under all scenarios.

4.2. Fixed design with missing data

This example further investigates the proposed method on a fixed design at two
different scales/resolutions. Noticing that missing areas or voxels may occur
within the domain, we mimic three types of missing schemes in this example,
including (i) complete data with no missing, (ii) missing at random, and (iii)
missing in a contiguous block as well as at random. Figures 4.3–4.5 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 Figures 4.4 (b) and (j) contains 12% of
the data, and thus the lowest missing rate in the second missing type is 12%.

Based on simulation example 1, we can see that the performance for various
methods is relatively consistent. Thus, in this example, we only consider the
first domain, $\Omega_{1}$, and try the same tetrahedral partitions as in
example 1. We consider two fixed resolutions/scales here: the lower
resolution/scale is $60\times 20\times 20$ with 22,160 voxels falling within
the domain, and the higher resolution/scale is $75\times 25\times 25$ with
42,600 voxels inside.

Similar to example 1, we calculate the average MISEs over 200 replications
with different missing types and missing rates and summarize them in Tables
4.2–4.3, while demonstrating them in Figure 4.6. In addition, Figures 4.3 and
4.4 show a typical simulation example for $n=50,000$ and PSNR $=5$ based on
different methods under different missing schemes. Based on these tables and
figures, one can clearly see that the prediction accuracy improves for all of
the methods as the missing rate decreases. The proposed methods outperform the
two traditional methods regardless of the type of missing scheme and the
missing rate.

Furthermore, the type of missing does not affect the two proposed methods very
much and the MISEs are very comparable for similar missing rates for the two
types of missing. However, the thin plate spline smoothing usually performs
better for missing at random, while the tensor product spline is better when
missing a contiguous block and at random.

Missing at Random
---
|  |  |
(a) True $m_{1}$ | (c) TPST 1 | (e) ATPST 1 | (g) Thin Plate
|  |  |
(b) Missing Data | (d) TPST 2 | (f) ATPST 2 | (h) Tensor Product
|  |  |
(i) True $m_{2}$ | (k) TPST 1 | (m) ATPST 1 | (o) Thin Plate
|  |  |
(j) Missing Data | (l) TPST 2 | (n) ATPST 2 | (p) Tensor Product
Figure 4.3: True and estimated functions when missing at random based on
different methods for a typical replication with $n=50,000$, missing rate
$=30\%$, and PSNR $=5$. Missing a Contiguous Block
---
|  |  |
(a) True $m_{1}$ | (c) TPST 1 | (e) ATPST 1 | (g) Thin Plate
|  |  |
(b) Missing Data | (d) TPST 2 | (f) ATPST 2 | (h) Tensor Product
|  |  |
(i) True $m_{2}$ | (k) TPST 1 | (m) ATPST 1 | (o) Thin Plate
|  |  |
(j) Missing Data | (l) TPST 2 | (n) ATPST 2 | (p) Tensor Product
Figure 4.4: True and estimated functions when missing a contiguous block based
on different methods for a typical replication with $n=50,000$, and PSNR $=5$.
Missing in Contiguous Block & at Random
---
|  |  |
(a) True $m_{1}$ | (c) TPST 1 | (e) ATPST 1 | (g) Thin Plate
|  |  |
(b) Missing Data | (d) TPST 2 | (f) ATPST 2 | (h) Tensor Product
|  |  |
(i) True $m_{2}$ | (k) TPST 1 | (m) ATPST 1 | (o) Thin Plate
|  |  |
(j) Missing Data | (l) TPST 2 | (n) ATPST 2 | (p) Tensor Product
Figure 4.5: True and estimated functions when missing a contiguous block & at
random based on different methods for a typical replication with $n=50,000$,
missing rate $=30\%$, and PSNR $=5$. Random Missing
---
Tube+Random Missing
Figure 4.6: Plots of average MISEs based on different methods under a fixed
design. Table 4.2: MISEs of different methods under fixed design with missing
at random.

$m$ | PSNR | Grid | Missing | TPST | ATPST | Thin | Tensor
---|---|---|---|---|---|---|---
Rate | $\triangle_{1}$ | $\triangle_{2}$ | $\triangle_{1}$ | $\triangle_{2}$ | Plate | Product
$m_{1}$ | 5 | 1 | 0 | 0.1370 | 0.1341 | 0.1230 | 0.1220 | 0.8397 | 0.9077
0.1 | 0.1453 | 0.1424 | 0.1313 | 0.1300 | 0.8465 | 0.9186
0.2 | 0.1555 | 0.1528 | 0.1414 | 0.1405 | 0.8554 | 0.9356
0.3 | 0.1688 | 0.1662 | 0.1545 | 0.1537 | 0.8672 | 0.9533
0.4 | 0.1864 | 0.1839 | 0.1719 | 0.1713 | 0.8831 | 0.9807
0.5 | 0.2048 | 0.2027 | 0.1909 | 0.1907 | 0.9027 | 1.0141
2 | 0 | 0.0944 | 0.0911 | 0.0828 | 0.0816 | 0.8091 | 0.8525
0.1 | 0.1002 | 0.0969 | 0.0884 | 0.0870 | 0.8138 | 0.8607
0.2 | 0.1068 | 0.1038 | 0.0947 | 0.0936 | 0.8191 | 0.8650
0.3 | 0.1161 | 0.1133 | 0.1037 | 0.1026 | 0.8260 | 0.8788
0.4 | 0.1270 | 0.1243 | 0.1143 | 0.1133 | 0.8350 | 0.8936
0.5 | 0.1410 | 0.1384 | 0.1280 | 0.1272 | 0.8468 | 0.9122
10 | 1 | 0 | 0.0601 | 0.0587 | 0.0613 | 0.0593 | 0.7925 | 0.8301
0.1 | 0.0639 | 0.0625 | 0.0652 | 0.0632 | 0.7952 | 0.8337
0.2 | 0.0684 | 0.0672 | 0.0700 | 0.0680 | 0.7988 | 0.8415
0.3 | 0.0745 | 0.0734 | 0.0763 | 0.0743 | 0.8035 | 0.8470
0.4 | 0.0821 | 0.0812 | 0.0839 | 0.0822 | 0.8098 | 0.8588
0.5 | 0.0913 | 0.0906 | 0.0931 | 0.0916 | 0.8177 | 0.8720
2 | 0 | 0.0412 | 0.0396 | 0.0416 | 0.0397 | 0.7807 | 0.8086
0.1 | 0.0438 | 0.0422 | 0.0443 | 0.0423 | 0.7826 | 0.8124
0.2 | 0.0469 | 0.0453 | 0.0475 | 0.0456 | 0.7847 | 0.8116
0.3 | 0.0510 | 0.0494 | 0.0518 | 0.0498 | 0.7874 | 0.8181
0.4 | 0.0559 | 0.0545 | 0.0569 | 0.0549 | 0.7910 | 0.8237
0.5 | 0.0623 | 0.0610 | 0.0636 | 0.0617 | 0.7958 | 0.8303
$m_{2}$ | 5 | 1 | 0 | 1.4608 | 1.4360 | 1.2021 | 1.1922 | 2.3709 | 3.3827
0.1 | 1.5668 | 1.5372 | 1.3097 | 1.2929 | 2.4189 | 3.0703
0.2 | 1.6951 | 1.6606 | 1.4232 | 1.4041 | 2.4760 | 3.2162
0.3 | 1.8524 | 1.8124 | 1.5607 | 1.5373 | 2.5584 | 3.0512
0.4 | 2.0584 | 2.0118 | 1.7274 | 1.7040 | 2.6645 | 3.7881
0.5 | 2.3183 | 2.2649 | 1.9632 | 1.9358 | 2.7985 | 3.5030
2 | 0 | 0.9754 | 0.9715 | 0.7414 | 0.7414 | 2.1619 | 3.0866
0.1 | 1.0415 | 1.0326 | 0.7953 | 0.7937 | 2.1910 | 2.9860
0.2 | 1.1235 | 1.1123 | 0.8655 | 0.8629 | 2.2263 | 2.8162
0.3 | 1.2248 | 1.2080 | 0.9575 | 0.9537 | 2.2689 | 2.9415
0.4 | 1.3558 | 1.3329 | 1.0845 | 1.0801 | 2.3296 | 3.1860
0.5 | 1.5255 | 1.4955 | 1.2532 | 1.2417 | 2.4061 | 2.9414
10 | 1 | 0 | 0.6116 | 0.6144 | 0.5656 | 0.5696 | 2.0510 | 2.9641
0.1 | 0.6506 | 0.6519 | 0.6023 | 0.6037 | 2.0688 | 2.6100
0.2 | 0.6978 | 0.6975 | 0.6434 | 0.6433 | 2.0907 | 2.7000
0.3 | 0.7572 | 0.7552 | 0.6948 | 0.6936 | 2.1213 | 2.4602
0.4 | 0.8361 | 0.8319 | 0.7668 | 0.7636 | 2.1615 | 3.0989
0.5 | 0.9410 | 0.9344 | 0.8655 | 0.8596 | 2.2123 | 2.6836
2 | 0 | 0.4374 | 0.4465 | 0.4046 | 0.4133 | 1.9770 | 2.8668
0.1 | 0.4603 | 0.4683 | 0.4244 | 0.4323 | 1.9879 | 2.7417
0.2 | 0.4891 | 0.4961 | 0.4495 | 0.4569 | 2.0011 | 2.5429
0.3 | 0.5255 | 0.5309 | 0.4823 | 0.4887 | 2.0168 | 2.6098
0.4 | 0.5721 | 0.5761 | 0.5246 | 0.5298 | 2.0386 | 2.8269
0.5 | 0.6490 | 0.6501 | 0.5956 | 0.5974 | 2.0753 | 2.5329

Table 4.3: MISEs of different methods under fixed design with missing both in
a contiguous block and at random.

$m$ | PSNR | Grid | Missing | TPST | ATPST | Thin | Tensor
---|---|---|---|---|---|---|---
Rate | $\triangle_{1}$ | $\triangle_{2}$ | $\triangle_{1}$ | $\triangle_{2}$ | Plate | Product
$m_{1}$ | 5 | 1 | 0.12 | 0.1643 | 0.1628 | 0.1561 | 0.1554 | 0.9323 | 0.8481
0.2 | 0.1710 | 0.1691 | 0.1625 | 0.1618 | 0.9439 | 0.8536
0.3 | 0.1812 | 0.1798 | 0.1737 | 0.1731 | 0.9598 | 0.8626
0.4 | 0.1985 | 0.1972 | 0.1895 | 0.1893 | 0.9852 | 0.8760
0.5 | 0.2171 | 0.2158 | 0.2085 | 0.2082 | 1.0117 | 0.8915
2 | 0.12 | 0.1168 | 0.1132 | 0.1081 | 0.1072 | 0.8728 | 0.8152
0.2 | 0.1224 | 0.1193 | 0.1138 | 0.1129 | 0.8769 | 0.8193
0.3 | 0.1303 | 0.1282 | 0.1226 | 0.1215 | 0.8874 | 0.8257
0.4 | 0.1397 | 0.1373 | 0.1312 | 0.1304 | 0.9004 | 0.8334
0.5 | 0.1537 | 0.1521 | 0.1453 | 0.1446 | 0.9189 | 0.8438
10 | 1 | 0.12 | 0.0794 | 0.0792 | 0.0837 | 0.0825 | 0.8435 | 0.7969
0.2 | 0.0827 | 0.0824 | 0.0870 | 0.0858 | 0.8499 | 0.7992
0.3 | 0.0883 | 0.0882 | 0.0930 | 0.0915 | 0.8560 | 0.8028
0.4 | 0.0963 | 0.0963 | 0.1010 | 0.0996 | 0.8656 | 0.8081
0.5 | 0.1057 | 0.1058 | 0.1102 | 0.1088 | 0.8758 | 0.8143
2 | 0.12 | 0.0555 | 0.0550 | 0.0582 | 0.0571 | 0.8211 | 0.7843
0.2 | 0.0584 | 0.0579 | 0.0613 | 0.0600 | 0.8210 | 0.7860
0.3 | 0.0628 | 0.0623 | 0.0659 | 0.0646 | 0.8250 | 0.7887
0.4 | 0.0676 | 0.0672 | 0.0709 | 0.0696 | 0.8298 | 0.7916
0.5 | 0.0746 | 0.0743 | 0.0781 | 0.0768 | 0.8370 | 0.7959
$m_{2}$ | 5 | 1 | 0.12 | 1.7386 | 1.7060 | 1.2766 | 1.2741 | 2.6896 | 2.4587
0.2 | 1.8429 | 1.8063 | 1.3653 | 1.3625 | 2.8285 | 2.4906
0.3 | 1.9842 | 1.9453 | 1.5069 | 1.4934 | 2.9021 | 2.5590
0.4 | 2.1575 | 2.1107 | 1.6649 | 1.6452 | 2.9690 | 2.6538
0.5 | 2.4060 | 2.3554 | 1.8783 | 1.8608 | 3.2190 | 2.7878
2 | 0.12 | 1.1926 | 1.1747 | 0.7875 | 0.7875 | 2.3785 | 2.2406
0.2 | 1.2536 | 1.2322 | 0.8332 | 0.8321 | 2.3984 | 2.2651
0.3 | 1.3544 | 1.3262 | 0.9096 | 0.9064 | 2.4597 | 2.3055
0.4 | 1.4779 | 1.4502 | 1.0130 | 1.0103 | 2.5459 | 2.3520
0.5 | 1.6444 | 1.6110 | 1.1742 | 1.1695 | 2.6820 | 2.4205
10 | 1 | 0.12 | 0.7024 | 0.7038 | 0.6167 | 0.6224 | 2.1975 | 2.1215
0.2 | 0.7394 | 0.7398 | 0.6505 | 0.6539 | 2.2987 | 2.1310
0.3 | 0.7904 | 0.7901 | 0.6970 | 0.6986 | 2.3045 | 2.1567
0.4 | 0.8546 | 0.8514 | 0.7534 | 0.7533 | 2.2841 | 2.1929
0.5 | 0.9470 | 0.9423 | 0.8374 | 0.8360 | 2.3999 | 2.2423
2 | 0.12 | 0.5122 | 0.5192 | 0.4437 | 0.4538 | 2.1131 | 2.0442
0.2 | 0.5328 | 0.5388 | 0.4618 | 0.4709 | 2.1114 | 2.0534
0.3 | 0.5670 | 0.5709 | 0.4899 | 0.4979 | 2.1341 | 2.0682
0.4 | 0.6093 | 0.6135 | 0.5270 | 0.5342 | 2.1744 | 2.0855
0.5 | 0.6673 | 0.6696 | 0.5804 | 0.5879 | 2.2360 | 2.1112

5\. Conclusions and Discussion

Challenges in 3D visualization and modern computer-aided design motivated the
advanced statistical analysis of 3D geometric data. Conventional smoothing
methods sometimes perform poorly over complex domains. The proposed approaches
are able to effectively denoise or deblur the collected data while preserving
inherent geometric features or spatial structures and provide multi-resolution
reconstruction. The experimental results demonstrate the effectiveness of the
proposed approaches in comparison to existing smoothing techniques.

Advancements in information technology have dramatically enhanced the
resolution of images. Even for a single high-quality image, the number of
observations is very high, and thus, the computational burden could also be
high. To address the “big” data issue, the divide and conquer algorithm is an
appealing method that recursively breaks down a problem into two or more sub-
problems of the same or related type until these become simple enough to be
solved directly. It then aggregates solutions to the sub-problems to produce a
global solution. Unlike tensor product spline or thin-plate spline, one unique
feature of the proposed TPST is its great scalability in computing, because
the spline basis function is generated restricted to each triangle without any
overlap, while the constraints on the spline coefficients achieve smoothness.
Based on this feature, a divide and conquer algorithm could be considered in
future works.

Acknowledgements

Li Wang’s research was partially supported by National Science Foundation
award DMS-1916204.

Supplemental Materials for “Spline Smoothing of 3D Geometric Data”

A. Detailed Introduction of Tetrahedral Partitions

A.1. Tetrahedral partitions

To recall you memory, we use $T$ to denote a tetrahedron, that is, a convex
hull of four noncoplanar points in $\mathbb{R}^{3}$.

######  Definition A.1 (Definition 16.5 in Lai and Schumaker (2007)).

A collection $\triangle=\\{T_{1},\ldots,T_{N}\\}$ of $N$ tetrahedra is called
a tetrahedral partition of $\Omega=\cup_{h=1}^{N}T_{h}$ provided that any pair
of tetrahedra in $\triangle$ intersect at most at a common vertex, along a
common edge, or along a common triangular face.

An example of tetrahedral partition is illustrated in Figure A.1(a), where
tetrahedron
$T_{1}=\langle\boldsymbol{v}_{5},\boldsymbol{v}_{1},\boldsymbol{v}_{4},\boldsymbol{v}_{3}\rangle$
and tetrahedron
$T_{2}=\langle\boldsymbol{v}_{2},\boldsymbol{v}_{1},\boldsymbol{v}_{3},\boldsymbol{v}_{4}\rangle$
form a tetrahedral partition of
$\Omega=\langle\boldsymbol{v}_{2},\boldsymbol{v}_{5},\boldsymbol{v}_{3},\boldsymbol{v}_{4}\rangle$.
In contrast, the partition in Figure A.1 (b) is not a tetrahedral partition,
in which tetrahedra
$T_{4}=\langle\boldsymbol{v}_{2},\boldsymbol{v}_{1},\boldsymbol{v}_{6},\boldsymbol{v}_{4}\rangle$
and
$T_{3}=\langle\boldsymbol{v}_{2},\boldsymbol{v}_{1},\boldsymbol{v}_{3},\boldsymbol{v}_{6}\rangle$
form a tetrahedral partition of
$\langle\boldsymbol{v}_{2},\boldsymbol{v}_{1},\boldsymbol{v}_{3},\boldsymbol{v}_{4}\rangle$,
but tetrahedra $T_{1}$ and $T_{3}$ and $T_{4}$ do not form a tetrahedral
partition of $\Omega$. In general, any type of polygon shapes can be used for
the partition of $\Omega$.

(a) A tetrahedral partition

(b) A non-tetrahedral partition

(c) Shape parameters

(d) Non-$\beta$-quasi-uniform tetrahedral partition

(e) Two tetrahedra sharing a common face

Figure A.1: (a) and (b) provides an illustration of an example and a
counterexample of tetrahedral partition; (c) shows an illustration of the
shape parameters; (d) shows an example of a non-$\beta$-quasi-uniform
tetrahedral partition; and (e) provides an example of two tetrahedra sharing a
common face.

######  Definition A.2 (Definition 15.4 in Lai and Schumaker (2007)).

A tetrahedron
$T=\langle\boldsymbol{v}_{1},\boldsymbol{v}_{2},\boldsymbol{v}_{3},\boldsymbol{v}_{4}\rangle$
with $\\{\boldsymbol{v}_{i}:=(x_{i},y_{i},z_{i})\\}_{i=1}^{4}$ is
nondegenerate provided that it has nonzero volume, and the vertices of $T$ are
in canonical order provided that if we rotate and translate $T$ so that the
triangular face
$\langle\boldsymbol{v}_{1},\boldsymbol{v}_{2},\boldsymbol{v}_{3}\rangle$ lies
in the $xy$-plane with
$\boldsymbol{v}_{1},\boldsymbol{v}_{2},\boldsymbol{v}_{3}$ in counterclockwise
order, then $z_{4}>0$.

######  Definition A.3 (Definition 16.1 in Lai and Schumaker (2007)).

Let $|T|$ be the length of the longest edge of $T$, and $\varrho_{T}$ be the
radius of the largest ball which can be inscribed in ${T}$(inradius of $T$),
then the ratio $\beta_{T}:={|T|}/{\varrho_{T}}$ is called the shape parameter
of $T$.

######  Remark A.1.

The shape parameter $\beta_{T}$ describes the shape of $T$. If $T$ is a
tetrahedron whose six edges are all of the same length, then
$\beta_{T}=2\sqrt{6}$. Such a tetrahedron is called a regular tetrahedron. For
any other tetrahedron, $\varrho_{T}$ is larger, and the larger $\varrho_{T}$
becomes, the flatter the tetrahedron $T$ becomes. Another way to describe the
shape of a tetrahedron is in terms of certain angles at the vertices of $T$.

######  Definition A.4.

A tetrahedral partition $\triangle$ is $\beta$-quasi-uniform if there is a
positive value $\beta$ such that $\triangle$ satisfies

$\displaystyle|\triangle|/\varrho_{T}\leq\beta<\infty,\quad\mbox{ for all
}T\in\triangle,$ (A.1)

where $|\triangle|:=\max\\{|T|,T\in\triangle\\}$ is referred as the size of
$\triangle$ , i.e., the length of the longest edge of $\triangle$.

Given a domain $\Omega$, we seek a tetrahedral partition whose minimum angle
is as large as possible. When $\beta$ is small, the tetrahedra are relatively
uniform in the sense that all angles of tetrahedra in the tetrahedral
partition $\triangle$ are relatively the same. Thus, the tetrahedral partition
looks more like a uniform tetrahedral partition which has the largest minimum
angle. Figure A.1 (d) gives an example of non-$\beta$-quasi-uniform
tetrahedral partition

$\displaystyle\triangle=\\{$ $\displaystyle
T_{1}=\langle\boldsymbol{v}_{2},\boldsymbol{v}_{1},\boldsymbol{c}_{1},\boldsymbol{v}_{4}\rangle,T_{2}=\langle\boldsymbol{v}_{1},\boldsymbol{c}_{2},\boldsymbol{c}_{1},\boldsymbol{v}_{4}\rangle,T_{3}=\langle\boldsymbol{c}_{1},\boldsymbol{v}_{1},\boldsymbol{c}_{3},\boldsymbol{c}_{2}\rangle,$
$\displaystyle
T_{4}=\langle\boldsymbol{c}_{3},\boldsymbol{v}_{1},\boldsymbol{c}_{4},\boldsymbol{c}_{2}\rangle,T_{5}=\langle\boldsymbol{c}_{3},\boldsymbol{v}_{1},\boldsymbol{v}_{3},\boldsymbol{c}_{4}\rangle,\ldots\\},$

where the edges $\langle\boldsymbol{v}_{4},\boldsymbol{c}_{1}\rangle$ and
$\langle\boldsymbol{c}_{2},\boldsymbol{c}_{3}\rangle$ are perpendicular to the
edge $\langle\boldsymbol{v}_{2},\boldsymbol{v}_{3}\rangle$, and
$\langle\boldsymbol{c}_{1},\boldsymbol{c}_{2}\rangle$ and
$\langle\boldsymbol{c}_{3},\boldsymbol{c}_{4}\rangle$ are perpendicular to
$\langle\boldsymbol{v}_{3},\boldsymbol{v}_{4}\rangle$.

Let $N$ be the number of the tetrahedra in the polygonal domain $\Omega$. From
(A.1), we have $N\leq(4\pi|\triangle|^{3})^{-1}3V_{\Omega}\beta^{3}$, where
$V_{\Omega}$ denotes the volume of $\Omega$.

A.2. Computer representation of tetrahedra partitions

A widely adopted data structure for efficiently storing a tetrahedral mesh is
using two separate matrices for nodes and tetrahedra elements in the mesh. We
first use a 3-column matrix to store all the nodes in the partition, where
each row corresponds to the Cartesian coordinates $(x,y,z)$ of one node. Then,
each tetrahedron within the partition can be represented concisely by
specifying the indices of the four nodes that define the tetrahedron.
Therefore, the list of tetrahedra elements can be simplified to a matrix with
four columns.

######  Example A.1.

Two adjacent tetrahedra
$T=\langle\boldsymbol{v}_{2},\boldsymbol{v}_{1},\boldsymbol{v}_{3},\boldsymbol{v}_{4}\rangle$
and
$\widetilde{T}=\langle\boldsymbol{v}_{5},\boldsymbol{v}_{1},\boldsymbol{v}_{4},\boldsymbol{v}_{3}\rangle$
share a common triangular face
$F=\langle\boldsymbol{v}_{1},\boldsymbol{v}_{3},\boldsymbol{v}_{4}\rangle$, as
illustrated in Figure A.1 (e), where the Cartesian coordinates of the five
vertices are $\boldsymbol{v}_{1}=(0,0,0)$, $\boldsymbol{v}_{2}=(1,0,0)$,
$\boldsymbol{v}_{3}=(0,1,0)$, $\boldsymbol{v}_{4}=(0,0,1)$,
$\boldsymbol{v}_{5}=(-1,0,0)$, respectively.

A.3. Delaunay tetrahedral partitions

Delaunay triangulation is one of the central topics in mesh generation due to
its theoretical properties that have been well-studied in 2D space (Chew,
1989; Lee and Schachter, 1980) and can be extended to 3D or even higher. In
this section, we mainly focus on Delaunay tetrahedral partition in the 3D
Euclidean space. Given a set of vertices $\mathcal{V}$, a tetrahedral
partition $\triangle$ is called Delaunay tetrahedral partition if for every
tetrahedron $T$ in $\triangle$, there is no vertex
$\boldsymbol{v}\in\mathcal{V}$ lying in the interior of the ball passing
through the vertices of $T$. The 3D Delaunay tetrahedral partition maximizes
the minimal solid angle appearing in the tetrahedra in $\triangle$. Chapter 16
in Lai and Schumaker (2007) provides a detailed discussion of various
theoretical properties for tetrahedral partitions.

|
---|---

Figure A.2: Tetrahedral partitions that can/cannot form a Delaunay partition.

Due to the mathematical guarantees of mesh quality, there have been extensive
efforts to develop Delaunay-based algorithms that can efficiently handle
arbitrarily complicated geometry for 2D domains (Ruppert, 1995; Shewchuk,
1996). However, extending these algorithms to 3D is challenging. In the past
three decades, many studies have been conducted regarding the construction of
tetrahedral partitions given complex domains (Shewchuk, 1998; Shewchuk and Si,
2014; Cheng et al., 2005; Jamin et al., 2015). Different in the ways of
handling domain boundaries and improving mesh quality, the tetrahedralization
methodologies can be roughly characterized into three groups: boundary-
constrained, Delaunay refinement, and variational methods; see Si (2015) for
more detailed descriptions.

The Delaunay-based mesh generation algorithms for 3D domains have been
embedded in computational softwares, for example, the delaunayTriangulaion and
delaunay functions in MATLAB. It has also been implemented in various software
packages and toolboxes, including QualMesh (Cheng et al., 2005), the CGAL
library (Jamin et al., 2015; The CGAL Project, 2020) for C++, TetGen
program(Si, 2015), iso2mesh (Fang and Boas, 2009) built on CGAL and TetGen for
MATLAB/Octave, and some MATLAB functions such as distmesh/distmeshnd (Persson
and Strang, 2004). These implementations differ in their methodologies of
handling domain boundaries, generating quality meshes, and input formats, and
therefore are suitable for various applications. For example, the functions
distmesh and distmeshnd take distance functions as representations of boundary
geometry, while delaunayTriangulation generates Delaunay 2D/3D triangulation
based on given point sets. If mathematical representations of boundary
geometry or the vertex sets for Delaunay triangulations are not immediately
available for complex real 3D objects, one can also use TetGen which requires
triangular surface mesh on domain boundaries as input.

We illustrate the overall procedure of mesh generation for a general 3D object
with the above tools using the example of a cow in Figure 2.1. Given a 3D
object, we first identify the boundary of the domain of interest and generate
a triangular surface mesh of the boundary using CGALmesh or iso2mesh. Then, a
tetrahedral partition for the complex object can be generated using TetGen or
iso2mesh. As one of the most commonly used programs for generating quality
tetrahedral meshes, TetGen is based on a mixture of the boundary constrained
method (George et al., 1991) and the Delaunay refinement method (Ruppert,
1995; Shewchuk, 1998). It takes triangular surface meshes as inputs and then
conducts Delaunay tetrahedralization followed by constrained Delaunay
tetrahedralization to preserve some desired constraints, such as domain
boundaries and special geometrical features. After these steps, it further
applies the Delaunay refinement techniques in boundary constrained methods to
improve the quality of tetrahedralizations based on user-specified criteria.
The introduced programs and toolboxes have greatly facilitated the
construction of tetrahedral meshes used in a wide range of applications, such
as medical imaging analysis, engineering, computer science, and others. Figure
1.1 presents a few mesh examples we generated with these tools.

B. 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 B.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 B.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.

B.1. Barycentric coordinates and Bernstein basis polynomials

Given a tetrahedron
$T=\langle\boldsymbol{v}_{1},\boldsymbol{v}_{2},\boldsymbol{v}_{3},\boldsymbol{v}_{4}\rangle$,
any fixed point $\boldsymbol{v}:=(x,y,z)\in\mathbb{R}^{3}$ has a unique
representation in terms of
$\langle\boldsymbol{v}_{1},\boldsymbol{v}_{2},\boldsymbol{v}_{3},\boldsymbol{v}_{4}\rangle$,

$\boldsymbol{v}=b_{1}\boldsymbol{v}_{1}+b_{2}\boldsymbol{v}_{2}+b_{3}\boldsymbol{v}_{3}+b_{4}\boldsymbol{v}_{4},~{}~{}\text{
with }b_{1}+b_{2}+b_{3}+b_{4}=1,$

where $(b_{1},b_{2},b_{3},b_{4})$ are called the barycentric coordinates of
$\boldsymbol{v}$ relative to the tetrahedron $T$. When the point
$\boldsymbol{v}$ is inside or on the faces of $T$, all $b_{1}$, $b_{2}$,
$b_{3}$ and $b_{4}$ are nonnegative. By Cramer’s rule, the barycentric
coordinate corresponding to vertex $\boldsymbol{v}_{i}$, satisfies
$b_{i}={\mathrm{det}(\mathbf{M}_{i})}/{\mathrm{det}(\mathbf{M})}$,
$i=1,\ldots,4$, where

$\mathbf{M}:=\begin{pmatrix}1&1&1&1\\\ x_{1}&x_{2}&x_{3}&x_{4}\\\
y_{1}&y_{2}&y_{3}&y_{4}\\\ z_{1}&z_{2}&z_{3}&z_{4}\end{pmatrix}$

and $\mathbf{M}_{i}$ replaces $\mathbf{M}$’s $i$-th column with $(1\ x\ y\
z)^{\top}$.

The barycentric coordinates $(b_{1},b_{2},b_{3},b_{4})$ also have an
interesting geometric interpretation. As shown in Figure B.1, for any
$\boldsymbol{v}\in T$, it divides the tetrahedron $T$ to four sub-tetrahedra,
$T_{1}$, $T_{2}$, $T_{3}$ and $T_{4}$. Notice that
$\mathrm{det}(\mathbf{M})=6V_{T}$, then $b_{i}=V_{T_{i}}/V_{T}$, where $T_{i}$
replace the vertex $\boldsymbol{v}_{i}$ in $V_{T}$ with $\boldsymbol{v}$,
$i=1,\ldots,4$.

|
---|---

Figure B.1: Illustration of barycentric coordinates of a point
$\boldsymbol{v}$ in a tetrahedron $T$.

For a nondegenerate tetrahedron $T$ and a point $\boldsymbol{v}\in 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

$B_{ijkl}^{d,T}(\boldsymbol{v}):=\frac{d!}{i!j!k!l!}b_{1}^{i}b_{2}^{j}b_{3}^{k}b_{4}^{l},\text{
with }i+j+k+l=d.$

For any positive integer $d$ and tetrahedron $T$, let $\mathcal{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 $\mathcal{P}_{d}(T)$ is
$\binom{d+3}{3}$.

According to Theorem 15.8 in Lai and Schumaker (2007) and Lemma D.1 in the
Appendix, the set of Bernstein basis polynomials

$\boldsymbol{B}^{d}_{T}(\boldsymbol{v}):=\\{B_{ijkl}^{d,T}(\boldsymbol{v}):i,j,k,l\geq
0,i+j+k+l=d\\}$

forms a basis for the space of polynomials $\mathcal{P}_{d}(T)$. In addition,
Bernstein basis functions $\\{B_{ijkl}^{d,T}(\boldsymbol{v})\\}_{i+j+k+l=d}$
have the following properties:

  1. 1.

$\\{B_{ijkl}^{d,T}\\}$ form a partition of unity, i.e., for all
$\boldsymbol{v}\in T$, $\sum_{i+j+k+l=d}B_{ijkl}^{d,T}(\boldsymbol{v})=1$;

  2. 2.

for all $\boldsymbol{v}\in T$, $0\leq B_{ijkl}^{d,T}(\boldsymbol{v})\leq 1$;

  3. 3.

$B_{ijkl}^{d,T}$ has a unique maximum at the point
$d^{-1}(i\boldsymbol{v}_{1}+j\boldsymbol{v}_{2}+k\boldsymbol{v}_{3}+l\boldsymbol{v}_{4})$.

######  Remark B.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.

To further illustrate the Bernstein basis functions, we present an example of
$\\{B_{ijkl}^{d,T}(\boldsymbol{v})\\}_{i+j+k+l=d}$ for $d=4$ in Figure B.2. In
this example, there are $\binom{7}{3}=35$ basis functions in total. Note that
all the function values vary between $0$ and $1$, and the colors scale to the
quantiles of function values.

 |   |   |   |   |
---|---|---|---|---|---
$B_{4000}^{4,T}=b_{1}^{4}$  |  $B_{0400}^{4,T}=b_{2}^{4}$  |  $B_{0040}^{4,T}=b_{3}^{4}$  |  $B_{0004}^{4,T}=b_{4}^{4}$  |  $B_{1111}^{4,T}=24b_{1}b_{2}b_{3}b_{4}$  |  Color bar
 |   |   |   |   |
$B_{3100}^{4,T}=4b_{1}^{3}b_{2}$  |  $B_{3010}^{4,T}=4b_{1}^{3}b_{3}$  |  $B_{3001}^{4,T}=4b_{1}^{3}b_{4}$  |  $B_{1300}^{4,T}=4b_{1}b_{2}^{3}$  |  $B_{1030}^{4,T}=4b_{1}b_{3}^{3}$  |  $B_{1003}^{4,T}=4b_{1}b_{4}^{3}$
 |   |   |   |   |
$B_{0310}^{4,T}=4b_{2}^{3}b_{3}$  |  $B_{0301}^{4,T}=4b_{2}^{3}b_{4}$  |  $B_{0130}^{4,T}=4b_{2}b_{3}^{3}$  |  $B_{0103}^{4,T}=4b_{2}b_{4}^{3}$  |  $B_{0031}^{4,T}=4b_{3}^{3}b_{4}$  |  $B_{0013}^{4,T}=4b_{3}b_{4}^{3}$
 |   |   |   |   |
$B_{2200}^{4,T}=6b_{1}^{2}b_{2}^{2}$  |  $B_{2020}^{4,T}=6b_{1}^{2}b_{3}^{2}$  |  $B_{2002}^{4,T}=6b_{1}^{2}b_{4}^{2}$  |  $B_{0220}^{4,T}=6b_{2}^{2}b_{3}^{2}$  |  $B_{0202}^{4,T}=6b_{2}^{2}b_{4}^{2}$  |  $B_{0022}^{4,T}=6b_{3}^{2}b_{4}^{2}$
 |   |   |   |   |
$B_{2110}^{4,T}=12b_{1}^{2}b_{2}b_{3}$  |  $B_{2101}^{4,T}=12b_{1}^{2}b_{2}b_{4}$  |  $B_{2011}^{4,T}=12b_{1}^{2}b_{3}b_{4}$  |  $B_{1210}^{4,T}=12b_{1}b_{2}^{2}b_{3}$  |  $B_{1201}^{4,T}=12b_{1}b_{2}^{2}b_{4}$  |  $B_{1120}^{4,T}=12b_{1}b_{2}b_{3}^{2}$
 |   |   |   |   |
$B_{1102}^{4,T}=12b_{1}b_{2}b_{4}^{2}$  |  $B_{1021}^{4,T}=12b_{1}b_{3}^{2}b_{4}$  |  $B_{1012}^{4,T}=12b_{1}b_{3}b_{4}^{2}$  |  $B_{0211}^{4,T}=12b_{2}^{2}b_{3}b_{4}$  |  $B_{0121}^{4,T}=12b_{2}b_{3}^{2}b_{4}$  |  $B_{0112}^{4,T}=12b_{2}b_{3}b_{4}^{2}$

Figure B.2: Illustrations of Bernstein basis functions
$\\{B_{ijkl}^{4,T}(\boldsymbol{v})\\}$.

Thus, given Bernstein basis functions
$\\{B_{ijkl}^{d,T}(\boldsymbol{v})\\}_{i+j+k+l=d}$, any polynomial
$p(\boldsymbol{v})\in\mathcal{P}_{d}(T)$ can be written uniquely as B-form,

$p(\boldsymbol{v})|_{T}=\sum_{i+j+k+l=d}\gamma_{T;ijkl}B_{ijkl}^{d,T}(\boldsymbol{v})=\boldsymbol{B}^{d}_{T}(\boldsymbol{v})^{\top}\boldsymbol{\gamma}_{T},$
(B.1)

where the coefficients
$\boldsymbol{\gamma}_{T}=\\{\gamma_{T;ijkl}\\}_{i+j+k+l=d}$ are called
B-coefficients of $p$. For the purpose of computer implementation, in this
paper, we employ the lexicographical order for ordering of the coefficients
$\boldsymbol{\gamma}_{T}$. To be specific, $\gamma_{T;ijkl}$ orders ahead of
$\gamma_{T;i^{\prime}j^{\prime}k^{\prime}l^{\prime}}$ either (i)
$i>i^{\prime}$, or (ii) $i=i^{\prime}$ and $j>j^{\prime}$, or (iii)
$i=i^{\prime}$, $j=j^{\prime}$ and $k>k^{\prime}$, or (iv) $i=i^{\prime}$,
$j=j^{\prime}$, $k=k^{\prime}$ and $l>l^{\prime}$. Consequently, we can
express

$\boldsymbol{B}^{d}_{T}(\boldsymbol{v})=\left(B_{d,0,0,0}^{d,T}(\boldsymbol{v}),B_{d-1,1,0,0}^{d,T}(\boldsymbol{v}),B_{d-1,0,1,0}^{d,T}(\boldsymbol{v}),\ldots,B_{0,0,1,d-1}^{d,T}(\boldsymbol{v}),B_{0,0,0,d}^{d,T}(\boldsymbol{v})\right)^{\top}$

and

$\boldsymbol{\gamma}_{T}=\left(\gamma_{T;d,0,0,0},\gamma_{T;d-1,1,0,0},\gamma_{T;d-1,0,1,0},\ldots,\gamma_{T;0,0,1,d-1},\gamma_{T;0,0,0,d}\right)^{\top}.$
(B.2)

Accordingly, in (B.2), the index of the element $\gamma_{T;ijkl}$ in the
vector $\boldsymbol{\gamma}_{T}$ is:

$\sum_{m=0}^{d-i}\frac{(m+1)m}{2}+\sum_{n=0}^{d-i-j}(n+1)-k.$

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. Consider Example A.1 illustrated in Figure A.1(e), for two
tetrahedra
$T=\langle\boldsymbol{v}_{2},\boldsymbol{v}_{1},\boldsymbol{v}_{3},\boldsymbol{v}_{4}\rangle$
and
$\widetilde{T}=\langle\boldsymbol{v}_{5},\boldsymbol{v}_{1},\boldsymbol{v}_{4},\boldsymbol{v}_{3}\rangle$
sharing a common triangular face
$F=\langle\boldsymbol{v}_{1},\boldsymbol{v}_{3},\boldsymbol{v}_{4}\rangle$,
assume two sets of Bernstein polynomial basis
$\\{B_{ijkl}^{d,T}(\boldsymbol{v})\\}_{i+j+k+l=d}$ and
$\\{\widetilde{B}_{ijkl}^{d,\widetilde{T}}(\boldsymbol{v})\\}_{i+j+k+l=d}$ are
defined on $T$ and $\widetilde{T}$ using the barycentric coordinates,
respectively. Consider two degree-$d$ polynomials $p(\boldsymbol{v})$ and
$\widetilde{p}(\widetilde{\boldsymbol{v}})$ defined on $T$ and
$\widetilde{T}$, respectively, with B-forms

$\displaystyle
p(\boldsymbol{v})=\sum_{i+j+k+l=d}\gamma_{ijkl}B_{ijkl}^{d,T}(\boldsymbol{v}),~{}~{}~{}~{}\widetilde{p}(\widetilde{\boldsymbol{v}})=\sum_{i+j+k+l=d}\widetilde{\gamma}_{ijkl}\widetilde{B}_{ijkl}^{d,\widetilde{T}}(\widetilde{\boldsymbol{v}}).$

For $\boldsymbol{w}\in F$, the barycentric coordinates with respect to $T$ and
$\widetilde{T}$ are $(0,b_{1},b_{3},1-b_{1}-b_{3})$ and
$(0,b_{1},1-b_{1}-b_{3},b_{3})$, respectively. Accordingly, we have

$\displaystyle p(\boldsymbol{w})=$
$\displaystyle\sum_{j+k+l=d}\gamma_{0jkl}\frac{d!}{j!k!l!}b_{1}^{j}b_{3}^{k}(1-b_{1}-b_{3})^{l},$
$\displaystyle\widetilde{p}(\boldsymbol{w})=$
$\displaystyle\sum_{j+k+l=d}\widetilde{\gamma}_{0jkl}\frac{d!}{j!k!l!}b_{1}^{j}(1-b_{1}-b_{3})^{k}b_{3}^{l}.$
(B.3)

Therefore, $p$ and $\widetilde{p}$ are continuous on $F$ if and only if

$\gamma_{0jkl}=\widetilde{\gamma}_{0jlk}$ (B.4)

for $j,k,l\geq 0$ and $j+k+l=d$.

B.2. Directional derivatives and smoothness

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 $p$, the
directional derivative at point $\boldsymbol{v}$ with respect to direction
$\boldsymbol{u}$ is defined as

$D_{\boldsymbol{u}}p(\boldsymbol{v}):=\frac{\partial}{\partial
t}p(\boldsymbol{v}+t\boldsymbol{u})\Big{|}_{t=0}=\lim_{t\rightarrow
0}\frac{p(\boldsymbol{v}+t\boldsymbol{u})-p(\boldsymbol{v})}{t}.$

Accordingly, for vector $\boldsymbol{u}:=(u_{x},u_{y},u_{z})\in\mathbb{R}^{3}$
and trivariate function $p$, the directional derivative at
$\boldsymbol{v}=(x,y,z)$ is

$D_{\boldsymbol{u}}p(x,y,z):=\frac{\partial}{\partial
t}p(x+tu_{x},y+tu_{y},z+tu_{z})\Big{|}_{t=0}.$

######  Remark B.2.

If $p$ is a polynomial of degree $d$, by calculus,
$D_{\boldsymbol{u}}p(x,y,z)=u_{x}D_{x}p(x,y,z)+u_{y}D_{y}p(x,y,z)+u_{z}D_{z}p(x,y,z)$,
so $D_{\boldsymbol{u}}p$ is a polynomial of $d-1$.

Consider direction $\boldsymbol{u}=\boldsymbol{w}_{1}-\boldsymbol{w}_{2}$,
where for $i=1,2$, $\boldsymbol{w}_{i}\in\mathbb{R}^{3}$ have barycentric
coordinates $(w_{i1},w_{i2},w_{i3},w_{i4})$ with respect to $T$. Then
$\boldsymbol{u}$ is uniquely described by the directional coordinates
$\boldsymbol{a}=(a_{1},a_{2},a_{3},a_{4})=(w_{11}-w_{21},w_{12}-w_{22},w_{13}-w_{23},w_{14}-w_{24})$.
Obviously, $a_{1}+a_{2}+a_{3}+a_{4}=0$. Direct calculation gives the
directional derivative of the Bernstein basis polynomial $B_{ijkl}^{d}$.

######  Lemma B.1 (Lemma 15.12 in Lai and Schumaker (2007)).

Consider direction $\boldsymbol{u}$ with directional coordinates
$\boldsymbol{a}=(a_{1},a_{2},a_{3},a_{4})$. Then

$\displaystyle D_{\boldsymbol{u}}B_{ijkl}^{d}(\boldsymbol{v})=$ $\displaystyle
d\left\\{a_{1}B_{i-1,j,k,l}^{d-1}(\boldsymbol{v})+a_{2}B_{i,j-1,k,l}^{d-1}(\boldsymbol{v})\right.$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left.+a_{3}B_{i,j,k-1,l}^{d-1}(\boldsymbol{v})+a_{4}B_{i,j,k,l-1}^{d-1}(\boldsymbol{v})\right\\},$
(B.5)

for any $\boldsymbol{v}\in T$ and $i+j+k+l=d$.

Consequently, one can obtain the directional derivative for any trivariate
polynomial $p$.

######  Theorem B.1 (Theorems 15.13 and 15.14 in Lai and Schumaker (2007)).

Consider direction $\boldsymbol{u}$ with directional coordinates
$\boldsymbol{a}=(a_{1},a_{2},a_{3},a_{4})$. Then for any trivariate polynomial
$p$ with B-form in (B.1), the directional derivative is

$D_{\boldsymbol{u}}p(\boldsymbol{v})=d\sum_{i+j+k+l=d-1}\gamma_{ijkl}^{(1)}(\boldsymbol{a})B_{ijkl}^{d-1}(\boldsymbol{v}),$
(B.6)

where
$\gamma_{ijkl}^{(1)}(\boldsymbol{a})=a_{1}\gamma_{i+1,j,k,l}+a_{2}\gamma_{i,j+1,k,l}+a_{3}\gamma_{i,j,k+1,l}+a_{4}\gamma_{i,j,k,l+1}$.

In general, given $\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{m}$ with
associated directional coordinates
$\boldsymbol{a}^{(i)}=(a_{1}^{(i)},a_{2}^{(i)},a_{3}^{(i)},a_{4}^{(i)})$,
$i=1,\ldots,m$,

$D_{\boldsymbol{u}_{m}}\cdots
D_{\boldsymbol{u}_{1}}p(\boldsymbol{v})=\frac{d!}{(d-m)!}\sum_{i+j+k+l=d-m}\gamma_{ijkl}^{(m)}(\boldsymbol{a}^{(1)},\ldots,\boldsymbol{a}^{(m)})B_{ijkl}^{d-m}(\boldsymbol{v}),$
(B.7)

where the coefficients are defined recursively as follows:

$\displaystyle\gamma_{ijkl}^{(m)}(\boldsymbol{a}^{(1)},\ldots,\boldsymbol{a}^{(m)})=$
$\displaystyle
a_{1}^{(m)}\gamma_{i+1,j,k,l}^{(m-1)}(\boldsymbol{a}^{(1)},\ldots,\boldsymbol{a}^{(m-1)})+a_{2}^{(m)}\gamma_{i,j+1,k,l}^{(m-1)}(\boldsymbol{a}^{(1)},\ldots,\boldsymbol{a}^{(m-1)})$
$\displaystyle+a_{3}^{(m)}\gamma_{i,j,k+1,l}^{(m-1)}(\boldsymbol{a}^{(1)},\ldots,\boldsymbol{a}^{(m-1)})+a_{4}^{(m)}\gamma_{i,j,k,l+1}^{(m-1)}(\boldsymbol{a}^{(1)},\ldots,\boldsymbol{a}^{(m-1)}),$

for $i=1,\ldots,m$, $m=1,\ldots,d$, with
$\gamma_{ijkl}^{(0)}(\boldsymbol{a})=\gamma_{ijkl}$.

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.1 (e).

######  Theorem B.2.

Suppose $\\{\gamma_{ijkl}\\}$ and $\\{\widetilde{\gamma}_{ijkl}\\}$ are
B-coefficients of $p$ and $\widetilde{p}$ relative to two tetrahedra
$T=\langle\boldsymbol{v}_{2},\boldsymbol{v}_{1},\boldsymbol{v}_{3},\boldsymbol{v}_{4}\rangle$
and
$\widetilde{T}=\langle\boldsymbol{v}_{5},\boldsymbol{v}_{1},\boldsymbol{v}_{4},\boldsymbol{v}_{3}\rangle$,
respectively, where $T$ and $\widetilde{T}$ share a common face
$F=\langle\boldsymbol{v}_{1},\boldsymbol{v}_{3},\boldsymbol{v}_{4}\rangle$.
Then the following statements are equivalent:

  1. 1.

$p$ and $\widetilde{p}$ join together with $\mathcal{C}^{r}$ continuity across
the face $F$;

  2. 2.

For all $\boldsymbol{w}\in F$, $m=0,\ldots,r$ and for all directions
$\boldsymbol{u}$,

$D_{\boldsymbol{u}}^{m}p(\boldsymbol{w})=D_{\boldsymbol{u}}^{m}\widetilde{p}(\boldsymbol{w});$

  3. 3.

For $p$ and $\widetilde{p}$ with B-forms in (B.3), for $i+j+k=d-m$,
$m=0,\ldots,r$,

$\widetilde{\gamma}_{mijk}=\sum_{\nu+\mu+\kappa+\delta=m}\gamma_{\nu,i+\mu,k+\kappa,j+\delta}B_{\nu\mu\kappa\delta}^{m}(\boldsymbol{v}_{5}).$
(B.8)

Consider the case of $d=2$ piecewise polynomial. There are in total 10
Bernstein basis polynomials with coefficients

$\displaystyle\\{\gamma_{ijkl}\\}$
$\displaystyle=(\gamma_{2000},\gamma_{1100},\gamma_{1010},\gamma_{1001},\gamma_{0200},\gamma_{0110},\gamma_{0101},\gamma_{0020},\gamma_{0011},\gamma_{0002})^{\top},$
$\displaystyle\\{\widetilde{\gamma}_{ijkl}\\}$
$\displaystyle=(\widetilde{\gamma}_{2000},\widetilde{\gamma}_{1100},\widetilde{\gamma}_{1010},\widetilde{\gamma}_{1001},\widetilde{\gamma}_{0200},\widetilde{\gamma}_{0110},\widetilde{\gamma}_{0101},\widetilde{\gamma}_{0020},\widetilde{\gamma}_{0011},\widetilde{\gamma}_{0002})^{\top}.$
(B.9)

Note that the barycentric coordinate of $\boldsymbol{v}_{5}$ with respect to
$T$ is $(2,-1,0,0)$.

If the trivariate polynomial is continuous over the whole region, then
applying (B.8) for $m=r=0$ generates

$\widetilde{\gamma}_{0ijk}=\gamma_{0ikj}B_{0000}^{0}(\boldsymbol{v}_{5})=\gamma_{0ikj},~{}~{}~{}~{}i+j+k=d,$
(B.10)

which matches the conclusion in (B.4).

Based on the de Casteljau Algorithm (Theorem D.1 in Appendix), we propose the
following computationally efficient algorithm to calculate the derivatives of
$p(\boldsymbol{v})$.

Inputs: Polynomial with B-form
$p(\boldsymbol{v})=\sum_{i+j+k+l=d}\gamma_{ijkl}B_{ijkl}^{d}(\boldsymbol{v})$,
directions $\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{m}$ with associated
directional coordinates
$\boldsymbol{a}^{(ii)}=\left(a_{1}^{(ii)},a_{2}^{(ii)},a_{3}^{(ii)},a_{4}^{(ii)}\right),ii=1,\ldots,m$

Initialize: $ii:=0$, $\gamma_{ijkl}^{(0)}:=\gamma_{ijkl}$

Outputs: $D_{\boldsymbol{u}_{m}}\ldots
D_{\boldsymbol{u}_{1}}p(\boldsymbol{v})$

for  _$ii=1,\ldots,m$ _ do

for  _$i+j+k+l=d-ii$_ do

$\displaystyle\gamma_{ijkl}^{(ii)}$
$\displaystyle(\boldsymbol{a}^{(1)},\ldots,\boldsymbol{a}^{(m)})$
$\displaystyle=$ $\displaystyle
a_{1}^{(ii)}\gamma_{i+1,j,k,l}^{(ii-1)}(\boldsymbol{a}^{(1)},\ldots,\boldsymbol{a}^{(ii-1)})+a_{2}^{(ii)}\gamma_{i,j+1,k,l}^{(ii-1)}(\boldsymbol{a}^{(1)},\ldots,\boldsymbol{a}^{(ii-1)})$
$\displaystyle+a_{3}^{(ii)}\gamma_{i,j,k+1,l}^{(ii-1)}(\boldsymbol{a}^{(1)},\ldots,\boldsymbol{a}^{(ii-1)})+a_{4}^{(ii)}\gamma_{i,j,k,l+1}^{(ii-1)}(\boldsymbol{a}^{(1)},\ldots,\boldsymbol{a}^{(ii-1)}),$

end for

end for

$D_{\boldsymbol{u}_{m}}\cdots
D_{\boldsymbol{u}_{1}}p(\boldsymbol{v})=\frac{d!}{(d-m)!}\sum_{i+j+k+l=d-m}\gamma_{ijkl}^{(m)}(\boldsymbol{a}^{(1)},\ldots,\boldsymbol{a}^{(m)})B_{ijkl}^{d-m}(\boldsymbol{v})$

Algorithm 1 Algorithm for Derivatives of $p(\boldsymbol{v})$.

B.3. Example of the constraint matrix

We illustrate the construction of the constraint matrix $\mathbf{H}$ using the
Example A.1 presented in Figure A.1 (e).

Continue our discussions in Section B.2 and consider the case of $d=2$
piecewise polynomial. Recall that there are in total 10 Bernstein basis
polynomials, with coefficients $\\{\gamma_{ijkl}\\}$ and
$\\{\widetilde{\gamma}_{ijkl}\\}$ in (B.9), and the barycentric coordinate of
$\boldsymbol{v}_{5}$ with respect to $T$ is $(2,-1,0,0)$.

As shown in (B.10), if the trivariate polynomial is continuous over the whole
region, then $\widetilde{\gamma}_{0ijk}=\gamma_{0ikj}$ for $i+j+k=d$. In this
case, we can write the constraint matrix $\mathbf{H}$ as

$\mathbf{H}=\left(\begin{array}[]{cccccccccccccccccccc}0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!1\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!-1\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\\
0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!1\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!-1\\!&\\!0\\!&\\!0\\!&\\!0\\\
0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!1\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!-1\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\\
0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!1\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!-1\\\
0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!1\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!-1\\!&\\!0\\\
0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!1\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!-1\\!&\\!0\\!&\\!0\\\
\end{array}\right).$ (B.11)

Furthermore, if the trivariate polynomial has continuous first derivatives
over the whole region, then (B.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$,

$\displaystyle\widetilde{\gamma}_{1ijk}=$
$\displaystyle\gamma_{1ikj}B_{1000}^{1}(\boldsymbol{v}_{5})+\gamma_{0,i+1,j}B_{0100}^{1}(\boldsymbol{v}_{5})+\gamma_{0,i,k+1,j}B_{0010}^{1}(\boldsymbol{v}_{5})+\gamma_{0,i,k,j+1}B_{0001}^{1}(\boldsymbol{v}_{5})$
$\displaystyle=$ $\displaystyle 2\gamma_{1ikj}-\gamma_{0,i+1,k,j}.$ (B.12)

Therefore, for $r=1$ where the trivariate spline has continuous first
derivatives over the whole region, we will obtain the following $\mathbf{H}$

$\mathbf{H}=\left(\begin{array}[]{cccccccccccccccccccc}0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!1\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!\\!-1\\!\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\\
0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!1\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!\\!-1\\!\\!&\\!0\\!&\\!0\\!&\\!0\\\
0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!1\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!\\!-1\\!\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\\
0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!1\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!-1\\\
0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!1\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!\\!-1\\!\\!&\\!0\\\
0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!1\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!\\!-1\\!\\!&\\!0\\!&\\!0\\\
0\\!&\\!2\\!&\\!0\\!&\\!0\\!&\\!\\!-1\\!\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!\\!-1\\!\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\\
0\\!&\\!0\\!&\\!2\\!&\\!0\\!&\\!0\\!&\\!\\!-1\\!\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!\\!-1\\!\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\\
0\\!&\\!0\\!&\\!0\\!&\\!2\\!&\\!0\\!&\\!0\\!&\\!\\!-1\\!\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!\\!-1\\!\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\!&\\!0\\\
\end{array}\right).$

Here, the first six rows match the form in (B.11), and the last three rows
correspond to (B.12).

C. Detailed Implementation of TPST

C.1. Directional derivatives for basis functions

Based on the conclusion in (B.7), the $m$th 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 $\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{m}$,
there exists a $\binom{d+3}{3}\times\binom{d+3-m}{3}$ matrix
$\mathbf{C}^{(m)}_{d}(\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{m})$ such that

$D_{\boldsymbol{u}_{m}}\cdots
D_{\boldsymbol{u}_{1}}\mathbf{B}_{d}(\boldsymbol{v})=\mathbf{C}^{(m)}_{d}(\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{m})\mathbf{B}_{d-m}(\boldsymbol{v}).$

######  Remark C.1.

Suppose $D_{\boldsymbol{u}_{m-1}}\cdots
D_{\boldsymbol{u}_{1}}\mathbf{B}_{d}(\boldsymbol{v})=\mathbf{C}^{(m-1)}_{d}(\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{m-1})\mathbf{B}_{d-m+1}(\boldsymbol{v})$,
then we have

$\displaystyle D_{\boldsymbol{u}_{m}}D_{\boldsymbol{u}_{m-1}}\cdots
D_{\boldsymbol{u}_{1}}\mathbf{B}_{d}(\boldsymbol{v})$
$\displaystyle=\mathbf{C}^{(m-1)}_{d}(\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{m-1})D_{\boldsymbol{u}_{m}}\mathbf{B}_{d-m+1}(\boldsymbol{v})$
$\displaystyle=\mathbf{C}^{(m-1)}_{d}(\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{m-1})\mathbf{C}^{(1)}_{d-m+1}(\boldsymbol{u}_{m})\mathbf{B}_{d-m}(\boldsymbol{v}),$

which implies
$\mathbf{C}^{(m)}_{d}(\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{m})=\mathbf{C}^{(m-1)}_{d}(\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{m-1})\mathbf{C}^{(1)}_{d-m+1}(\boldsymbol{u}_{m})$.
Keep decomposing matrix
$\mathbf{C}^{(m-1)}_{d}(\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{m-1})$, we
can have

$\mathbf{C}^{(m)}_{d}(\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{m})=\mathbf{C}^{(1)}_{d}(\boldsymbol{u}_{1})\cdots\mathbf{C}^{(1)}_{d-m+1}(\boldsymbol{u}_{m}).$

Next, we introduce the explicit form of matrix
$\mathbf{C}^{(1)}_{d}(\boldsymbol{u})$. Recall that
$D_{\boldsymbol{u}}B_{ijkl}^{d}(\boldsymbol{v})$ can be written as

$d\left\\{a_{1}B_{i-1,j,k,l}^{d-1}(\boldsymbol{v})+a_{2}B_{i,j-1,k,l}^{d-1}(\boldsymbol{v})+a_{3}B_{i,j,k-1,l}^{d-1}(\boldsymbol{v})+a_{4}B_{i,j,k,l-1}^{d-1}(\boldsymbol{v})\right\\},$

where $(a_{1},a_{2},a_{3},a_{4})$ is the barycentric coordinate of direction
$\boldsymbol{u}$. Then, the
$(\sum_{m=0}^{d-i}\frac{(m+1)m}{2}+\sum_{n=0}^{d-i-j}(n+1)-k)$th row of matrix
$\mathbf{C}^{(1)}_{d}(\boldsymbol{u})$ is
$d(a_{1}\boldsymbol{e}_{\mathcal{I}_{1}}+a_{2}\boldsymbol{e}_{\mathcal{I}_{2}}+a_{3}\boldsymbol{e}_{\mathcal{I}_{3}}+a_{4}\boldsymbol{e}_{\mathcal{I}_{4}})$,
where $\boldsymbol{e}_{\mathcal{I}}$ is unit vector with $\mathcal{I}$th
element being one, $\mathcal{I}_{1}$, $\mathcal{I}_{2}$, $\mathcal{I}_{3}$ and
$\mathcal{I}_{4}$ are indexes of basis functions
$B_{i-1,j,k,l}^{d-1}(\boldsymbol{v})$, $B_{i,j-1,k,l}^{d-1}(\boldsymbol{v})$,
$B_{i,j,k-1,l}^{d-1}(\boldsymbol{v})$, and
$B_{i,j,k,l-1}^{d-1}(\boldsymbol{v})$.

In the following, we provide a simple example for the first and second order
derivatives. Consider a case on the tetrahedron $T$ in the Figure B.1 with
$d=3$. The barycentric coordinates of
$x=\boldsymbol{v}_{1}-\boldsymbol{v}_{2}$. Then, the first and second order
derivatives of all the Berstein basis function are
$D_{x}\mathbf{B}_{3}(\boldsymbol{v})=\mathbf{C}^{(1)}_{3}(x)\mathbf{B}_{2}(\boldsymbol{v})$
and
$D_{xx}^{2}\mathbf{B}_{3}(\boldsymbol{v})=\mathbf{C}^{(2)}_{3}(x,x)\mathbf{B}_{1}(\boldsymbol{v})=\mathbf{C}^{(1)}_{3}(x)\mathbf{C}^{(1)}_{2}(x)\mathbf{B}_{1}(\boldsymbol{v})$,
where matricex $\mathbf{C}^{(1)}_{2}(x)$ and $\mathbf{C}^{(1)}_{3}(x)$ are

$\mathbf{C}^{(1)}_{2}(x)=\begin{pmatrix}2&0&0&0\\\ -2&2&0&0\\\ 0&0&2&0\\\
0&0&0&2\\\ 0&-2&0&0\\\ 0&0&-2&0\\\ 0&0&0&-2\\\ 0&0&0&0\\\ 0&0&0&0\\\
0&0&0&0\\\
\end{pmatrix},~{}\mathbf{C}^{(1)}_{3}(x)=\begin{pmatrix}3&0&0&0&0&0&0&0&0&0\\\
-3&3&0&0&0&0&0&0&0&0\\\ 0&0&3&0&0&0&0&0&0&0\\\ 0&0&0&3&0&0&0&0&0&0\\\
0&-3&0&0&3&0&0&0&0&0\\\ 0&0&-3&0&0&3&0&0&0&0\\\ 0&0&0&-3&0&0&3&0&0&0\\\
0&0&0&0&0&0&0&3&0&0\\\ 0&0&0&0&0&0&0&0&3&0\\\ 0&0&0&0&0&0&0&0&0&3\\\
0&0&0&0&-3&0&0&0&0&0\\\ 0&0&0&0&0&-3&0&0&0&0\\\ 0&0&0&0&0&0&-3&0&0&0\\\
0&0&0&0&0&0&0&-3&0&0\\\ 0&0&0&0&0&0&0&0&-3&0\\\ 0&0&0&0&0&0&0&0&0&-3\\\
0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0\\\ \end{pmatrix}.$

C.2. Details of constructing penalty matrix

In this section, we introduce the details of constructing penalty $\mathbf{P}$
matrix. Applying the results in Section C.1, the second order derivative of
Bernstein basis functions $\mathbf{B}_{d}(\boldsymbol{v})$ at directions
$\boldsymbol{u}_{1}$ and $\boldsymbol{u}_{2}$ are

$D_{\boldsymbol{u}_{2}}D_{\boldsymbol{u}_{1}}\mathbf{B}_{d}(\boldsymbol{v})=\mathbf{C}^{(2)}(\boldsymbol{u}_{1},\boldsymbol{u}_{2})\mathbf{B}_{d-2}(\boldsymbol{v}).$

Therefore, for directions $g,g^{\prime}=x,y,z$, we have

$\displaystyle\int_{T}\left\\{D^{2}_{gg^{\prime}}\mathbf{B}_{d}(\boldsymbol{v})\right\\}\left\\{D^{2}_{gg^{\prime}}\mathbf{B}_{d}(\boldsymbol{v})\right\\}^{\top}\mathrm{d}\boldsymbol{v}=\mathbf{C}^{(2)}(g,g^{\prime})\int_{T}\mathbf{B}_{d-2}(\boldsymbol{v})\mathbf{B}_{d-2}(\boldsymbol{v})^{\top}\mathrm{d}\boldsymbol{v}\\{\mathbf{C}_{d}^{(2)}(g,g^{\prime})\\}^{\top}.$

Let
$\mathbf{L}_{T}^{d-2}=\int_{T}\mathbf{B}_{d-2}(\boldsymbol{v})\mathbf{B}_{d-2}(\boldsymbol{v})^{\top}\mathrm{d}\boldsymbol{v}$
be the $\binom{d+1}{3}\times\binom{d+1}{3}$ matrix with entries
$\int_{T}B_{\nu\mu\kappa\delta}^{T,d-2}B_{\nu^{\prime}\mu^{\prime}\kappa^{\prime}\delta^{\prime}}^{T,d-2}(\boldsymbol{v})\mathrm{d}\boldsymbol{v}$.
By the Lemma 15.29 in Lai and Schumaker (2007), we have

$\displaystyle\int_{T}B_{\nu\mu\kappa\delta}^{T,d-2}(\boldsymbol{v})B_{\nu^{\prime}\mu^{\prime}\kappa^{\prime}\delta^{\prime}}^{T,d-2}(\boldsymbol{v})\mathrm{d}\boldsymbol{v}$
$\displaystyle=\frac{{\nu+\nu^{\prime}\choose\nu}{\mu+\mu^{\prime}\choose\mu}{\kappa+\kappa^{\prime}\choose\kappa}{\delta+\delta^{\prime}\choose\delta}}{{2d-4\choose
d-2}}\int_{T}B_{\nu+\nu^{\prime},\mu+\mu^{\prime},\kappa+\kappa^{\prime},\delta+\delta^{\prime}}^{T,2d-4}(\boldsymbol{v})\mathrm{d}\boldsymbol{v}$
$\displaystyle=\frac{{\nu+\nu^{\prime}\choose\nu}{\mu+\mu^{\prime}\choose\mu}{\kappa+\kappa^{\prime}\choose\kappa}{\delta+\delta^{\prime}\choose\delta}}{{2d-4\choose
d-2}{2d-1\choose 3}}V_{T},$

recall that $V_{T}$ is the volume of tetrahedron $T$. We finally obtain
$\mathbf{P}_{T}^{gg^{\prime}}=\mathbf{C}_{d}^{(2)}(g,g^{\prime})\mathbf{L}_{T}^{d-2}\\{\mathbf{C}_{d}^{(2)}(g,g^{\prime})\\}^{\top}$
and
$\mathbf{P}_{T}=\sum_{g,g^{\prime}\in\\{x,y,z\\}}\mathbf{P}_{T}^{gg^{\prime}}.$

Recall that $\mathcal{E}(s)=\sum_{T\in\triangle}\mathcal{E}(s_{T})$ and
$\mathcal{E}(s_{T})=\boldsymbol{\gamma}_{T}^{\top}\mathbf{P}_{T}\boldsymbol{\gamma}_{T}.$
Then, one can have
$\mathcal{E}(s)=\boldsymbol{\gamma}^{\top}\mathbf{P}\boldsymbol{\gamma}$,
where
$\boldsymbol{\gamma}=(\boldsymbol{\gamma}_{1}^{\top},\ldots,\boldsymbol{\gamma}_{N}^{\top})^{\top}$,
and $\mathbf{P}=\textrm{diag}(\mathbf{P}_{T},T\in\triangle)$ is the block
diagonal matrix.

C.3. An example of the penalty matrix

In this section, we provide a simple example for $\mathbf{P}_{T}$. Consider a
case on the tetrahedron $T$ in the Figure B.1 with $d=3$. The barycentric
coordinates of $x=\boldsymbol{v}_{1}-\boldsymbol{v}_{2}$,
$y=\boldsymbol{v}_{3}-\boldsymbol{v}_{2}$ and
$z=\boldsymbol{v}_{4}-\boldsymbol{v}_{2}$ are $(1,-1,0,0)$, $(0,-1,1,0)$, and
$(0,-1,0,1)$, respectively. The second directional derivative of these
Bernstein basis polynomials are the linear combination of
$B_{1000}^{T,1}(\boldsymbol{v}),\ldots,B_{0001}^{T,1}(\boldsymbol{v})$. For
example, we have
$D_{x}B_{2100}^{T,3}(\boldsymbol{v})=3B_{1100}^{T,2}(\boldsymbol{v})-3B_{2000}^{T,2}(\boldsymbol{v}),$
and
$D_{xx}^{2}B_{2100}^{T,3}(\boldsymbol{v})=-12B_{1000}^{T,1}(\boldsymbol{v})+6B_{0100}^{T,1}(\boldsymbol{v}).$
According to the results in Section C.1, the coefficient matrices for second
order derivatives are $\mathbf{C}_{3}^{(2)}(x,x)$,
$\mathbf{C}_{3}^{(2)}(y,y)$, $\mathbf{C}_{3}^{(2)}(z,z)$,
$\mathbf{C}_{3}^{(2)}(x,y)$, $\mathbf{C}_{3}^{(2)}(x,z)$, and
$\mathbf{C}_{3}^{(2)}(y,z)$

$\stackrel{{\scriptstyle\mbox{$\mathbf{C}_{3}^{(2)}(x,x)$}}}{{\begin{pmatrix}6&0&0&0\\\
-12&6&0&0\\\ 0&0&6&0\\\ 0&0&0&6\\\ 6&-12&0&0\\\ 0&0&-12&0\\\ 0&0&0&-12\\\
0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\ 0&6&0&0\\\ 0&0&6&0\\\ 0&0&0&6\\\ 0&0&0&0\\\
0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\
\end{pmatrix}}},\stackrel{{\scriptstyle\mbox{$\mathbf{C}_{3}^{(2)}(y,y)$}}}{{\begin{pmatrix}0&0&0&0\\\
0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\ 6&0&0&0\\\ -12&0&0&0\\\ 0&0&0&0\\\ 6&0&0&0\\\
0&0&0&0\\\ 0&0&0&0\\\ 0&6&0&0\\\ 0&-12&6&0\\\ 0&0&0&6\\\ 0&6&-12&0\\\
0&0&0&-12\\\ 0&0&0&0\\\ 0&0&6&0\\\ 0&0&0&6\\\ 0&0&0&0\\\ 0&0&0&0\\\
\end{pmatrix}}},\stackrel{{\scriptstyle\mbox{$\mathbf{C}_{3}^{(2)}(z,z)$}}}{{\begin{pmatrix}0&0&0&0\\\
0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\ 6&0&0&0\\\ 0&0&0&0\\\ -12&0&0&-12\\\
0&0&0&0\\\ 0&0&0&0\\\ 6&0&0&0\\\ 0&6&0&0\\\ 0&0&6&0\\\ 0&-12&0&6\\\ 0&0&0&0\\\
0&0&-12&0\\\ 0&6&0&-12\\\ 0&0&0&0\\\ 0&0&6&0\\\ 0&0&0&6\\\ 0&0&0&0\\\
\end{pmatrix}}},$

$\stackrel{{\scriptstyle\mbox{$\mathbf{C}_{3}^{(2)}(x,y)$}}}{{\begin{pmatrix}0&0&0&0\\\
-6&0&0&0\\\ 6&0&0&0\\\ 0&0&0&0\\\ 6&-6&0&0\\\ -6&6&-6&0\\\ 0&0&0&-6\\\
0&0&6&0\\\ 0&0&0&6\\\ 0&0&0&0\\\ 0&6&0&0\\\ 0&-6&6&0\\\ 0&0&0&6\\\ 0&0&-6&0\\\
0&0&0&-6\\\ 0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\
\end{pmatrix}}},~{}\stackrel{{\scriptstyle\mbox{$\mathbf{C}_{3}^{(2)}(x,z)$}}}{{\begin{pmatrix}0&0&0&0\\\
-6&0&0&0\\\ 0&0&0&0\\\ 6&0&0&0\\\ 6&-6&0&0\\\ 0&0&-6&0\\\ -6&6&0&-6\\\
0&0&0&0\\\ 0&0&6&0\\\ 0&0&0&6\\\ 0&6&0&0\\\ 0&0&6&0\\\ 0&-6&0&6\\\ 0&0&0&0\\\
0&0&-6&0\\\ 0&0&0&-6\\\ 0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\
\end{pmatrix}}},~{}\stackrel{{\scriptstyle\mbox{$\mathbf{C}_{3}^{(2)}(y,z)$}}}{{\begin{pmatrix}0&0&0&0\\\
0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\ 6&0&0&0\\\ -6&0&0&0\\\ -6&0&0&0\\\ 0&0&0&0\\\
6&0&0&0\\\ 0&0&0&0\\\ 0&6&0&0\\\ 0&-6&6&0\\\ 0&-6&0&6\\\ 0&0&-6&0\\\
0&6&-6&-6\\\ 0&0&0&-6\\\ 0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\
\end{pmatrix}}}.$

Applying the results in Section C.2, we have

$\mathbf{L}_{T}^{1}=\frac{1}{12}\left(\begin{array}[]{cccc}2&1&1&1\\\
1&2&1&1\\\ 1&1&2&1\\\ 1&1&1&2\\\ \end{array}\right),$

and

$\mathbf{P}_{T}=\sum_{g,g^{\prime}\in\\{x,y,z\\}}\mathbf{C}_{3}^{(2)}(g,g^{\prime})\mathbf{L}_{T}^{1}\\{\mathbf{C}_{3}^{(2)}(g,g^{\prime})\\}^{\top}.$

C.4. 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 $\lambda$. A large $\lambda$ tends to shrink
the second-order derivative function of splines functions to zero and gives a
smooth function with large fitting errors. On the other hand, a small
$\lambda$ may generate a more wiggly function and result in an overfitting
problem.

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, we adopt the block cross-validation strategy in
Roberts et al. (2017) and Valavi et al. (2019). 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 object data into small 3D blocks. Each tetrahedron
is considered as one single 3D block. Figure C.1 shows an illustration of
block CV using a tetrahedral partition. In Figure C.1, we divide the domain
into 384 equal-size tetrahedra and randomly assign these tetrahedra into five
folds with color indicating different folds.

Figure C.1: Illustration of block CV using tetrahedral partition. Tetrahedra
with same color belong to same folds.

D. Proof of Theorems

D.1. Theoretical Properties of Bernstein Bases

######  Lemma D.1.

For the Bernstein basis functions
$\\{B_{ijkl}^{d,T}(\boldsymbol{v})\\}_{i+j+k+l=d}$, each $B_{ijkl}^{d,T}$ has
a unique maximum at the point
$d^{-1}(i\boldsymbol{v}_{1}+j\boldsymbol{v}_{2}+k\boldsymbol{v}_{3}+l\boldsymbol{v}_{4})$.

To show Lemma D.1, we first need to show Lemma B.1.

###### Proof of Lemma B.1.

Suppose $\boldsymbol{b}=(b_{1},b_{2},b_{3},b_{4})$ are the barycentric
coordinates of $\boldsymbol{v}$. Then the barycentric coordinates of
$\boldsymbol{v}+t\boldsymbol{u}$ are
$(b_{1}+ta_{1},b_{2}+ta_{2},b_{3}+ta_{3},b_{4}+ta_{4})$. Thus, for
$i+j+k+l=d$,

$B_{ijkl}^{d}(\boldsymbol{v}+t\boldsymbol{u})=\frac{d!}{i!j!k!l!}(b_{1}+ta_{1})^{i}(b_{2}+ta_{2})^{j}(b_{3}+ta_{3})^{k}(b_{4}+ta_{4})^{l}.$

Hence,

$\displaystyle
D_{\boldsymbol{u}}B_{ijkl}^{d}(\boldsymbol{v})=\frac{\partial}{\partial
t}p(\boldsymbol{v}+t\boldsymbol{u})\Big{|}_{t=0}$ $\displaystyle=$
$\displaystyle\frac{d!}{i!j!k!l!}\left(ia_{1}b_{1}^{i-1}b_{2}^{j}b_{3}^{k}b_{4}^{l}+ja_{2}b_{1}^{i}b_{2}^{j-1}b_{3}^{k}b_{4}^{l}+ka_{3}b_{1}^{i}b_{2}^{j}b_{3}^{k-1}b_{4}^{l}+la_{4}b_{1}^{i}b_{2}^{j}b_{3}^{k}b_{4}^{l-1}\right)$
$\displaystyle=$ $\displaystyle
d\left\\{a_{1}B_{i-1,j,k,l}^{d-1}(\boldsymbol{v})+a_{2}B_{i,j-1,k,l}^{d-1}(\boldsymbol{v})+a_{3}B_{i,j,k-1,l}^{d-1}(\boldsymbol{v})+a_{4}B_{i,j,k,l-1}^{d-1}(\boldsymbol{v})\right\\}.$

∎

###### Proof of Lemma D.1.

For
$T=\langle\boldsymbol{v}_{1},\boldsymbol{v}_{2},\boldsymbol{v}_{3},\boldsymbol{v}_{4}\rangle$,
the barycentric coordinates for $\boldsymbol{v}_{1}$, $\boldsymbol{v}_{2}$,
$\boldsymbol{v}_{3}$ and $\boldsymbol{v}_{4}$ are $(1,0,0,0)$, $(0,1,0,0)$,
$(0,0,1,0)$ and $(0,0,0,1)$, respectively. Accordingly, the directional
coordinates of $\boldsymbol{u}_{1}=\boldsymbol{v}_{1}-\boldsymbol{v}_{2}$,
$\boldsymbol{u}_{2}=\boldsymbol{v}_{1}-\boldsymbol{v}_{3}$ and
$\boldsymbol{u}_{3}=\boldsymbol{v}_{1}-\boldsymbol{v}_{4}$ are $(1,-1,0,0)$,
$(1,0,-1,0)$ and $(1,0,0,-1)$, respectively. For $\boldsymbol{v}\in T$ with
barycentric coordinates $\boldsymbol{b}=(b_{1},b_{2},b_{3},b_{4})$, consider
derivatives of $B_{ijkl}^{d}(\boldsymbol{v})$ with respect to directions
$\boldsymbol{u}_{1}$, $\boldsymbol{u}_{2}$ and $\boldsymbol{u}_{3}$, then

$\displaystyle D_{\boldsymbol{u}_{1}}B_{ijkl}^{d}(\boldsymbol{v})$
$\displaystyle=B_{ijkl}^{d}(\boldsymbol{v})(ib_{1}^{-1}-jb_{2}^{-1}),~{}~{}~{}~{}D_{\boldsymbol{u}_{2}}B_{ijkl}^{d}(\boldsymbol{v})=B_{ijkl}^{d}(\boldsymbol{v})(ib_{1}^{-1}-kb_{3}^{-1}),$
$\displaystyle D_{\boldsymbol{u}_{3}}B_{ijkl}^{d}(\boldsymbol{v})$
$\displaystyle=B_{ijkl}^{d}(\boldsymbol{v})(ib_{1}^{-1}-lb_{4}^{-1}).$

Set these equations to zero and combine $b_{1}+b_{2}+b_{3}+b_{4}=1$ gives
$(b_{1},b_{2},b_{3},b_{4})=d^{-1}(i,j,k,l)$. ∎

D.2. Proof of Theorems B.1 and B.2

###### Proof of Theorem B.1.

For
$p(\boldsymbol{v})=\sum_{i+j+k+l=d}\gamma_{ijkl}B_{ijkl}^{d}(\boldsymbol{v})$,
by Lemma B.1,

$\displaystyle
D_{\boldsymbol{u}}p(\boldsymbol{v})=\sum_{i+j+k+l=d}\gamma_{ijkl}D_{\boldsymbol{u}}B_{ijkl}^{d}(\boldsymbol{v})$
$\displaystyle=$
$\displaystyle\sum_{i+j+k+l=d}\gamma_{ijkl}d\left\\{a_{1}B_{i-1,j,k,l}^{d-1}(\boldsymbol{v})+a_{2}B_{i,j-1,k,l}^{d-1}(\boldsymbol{v})+a_{3}B_{i,j,k-1,l}^{d-1}(\boldsymbol{v})+a_{4}B_{i,j,k,l-1}^{d-1}(\boldsymbol{v})\right\\}$
$\displaystyle=$ $\displaystyle
d\sum_{i+j+k+l=d-1}\left(a_{1}\gamma_{i+1,j,k,l}+a_{2}\gamma_{i,j+1,k,l}+a_{3}\gamma_{i,j,k+1,l}+a_{4}\gamma_{i,j,k,l+1}\right)B_{ijkl}^{d-1}(\boldsymbol{v}).$

Thus, (B.6) follows. Consequently, one can obtain (B.7) by repeatedly applying
(B.6) for directions $\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{m}$. ∎

To show Theorem B.2, we need the following results.

######  Theorem D.1 (Theorems 15.10 in Lai and Schumaker (2007)).

Suppose $p(\boldsymbol{v})$ is a trivariate polynomial with B-form
$p(\boldsymbol{v})=\sum_{i+j+k+l=d}\gamma_{ijkl}B_{ijkl}^{d}(\boldsymbol{v})$.
Define $\gamma_{ijkl}^{(0)}:=\gamma_{ijkl}$, $i+j+k+l=d$. Suppose
$\boldsymbol{v}$ has barycentric coordinates
$\boldsymbol{b}=(b_{1},b_{2},b_{3},b_{4})$. Then

$p(\boldsymbol{v})=\sum_{i+j+k+l=d-m}\gamma_{ijkl}^{(m)}B_{ijkl}^{d-m}(\boldsymbol{v}),$

where for $m=1,\ldots,d$, $\gamma_{ijkl}^{(m)}$ are computed by the recursion

$\gamma_{ijkl}^{(m)}=b_{1}\gamma_{i+1,j,k,l}^{(m-1)}+b_{2}\gamma_{i,j+1,k,l}^{(m-1)}+b_{3}\gamma_{i,j,k+1,l}^{(m-1)}+b_{4}\gamma_{i,j,k,l+1}^{(m-1)},$
(D.1)

for $i+j+k+l=d-m$.

######  Remark D.1.

The recursive formula (D.1) is also referred as de Casteljau algorithm. See
Section 15.6 in Lai and Schumaker (2007) for more details.

######  Lemma D.2.

The coefficients in the recursive formula of de Casteljau algorithm (D.1) are
given by

$\gamma_{ijkl}^{(m)}=\sum_{i^{\prime}+j^{\prime}+k^{\prime}+l^{\prime}=m}\gamma_{i+i^{\prime},j+j^{\prime},k+k^{\prime},l+l^{\prime}}B_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}^{m}(\boldsymbol{v}),$
(D.2)

where $i+j+k+l=d-m$.

###### Proof.

Define operators $E_{1}\gamma_{ijkl}=\gamma_{i+1,j,k,l}$,
$E_{2}\gamma_{ijkl}=\gamma_{i,j+1,k,l}$,
$E_{3}\gamma_{ijkl}=\gamma_{i,j,k+1,l}$ and
$E_{4}\gamma_{ijkl}=\gamma_{i,j,k,l+1}$. Thus, by (D.1), for $i+j+k+l=d-m$,

$\displaystyle\gamma_{ijkl}^{(m)}=$
$\displaystyle\left(b_{1}E_{1}+b_{2}E_{2}+b_{3}E_{3}+b_{4}E_{4}\right)\gamma_{ijkl}^{(m-1)}=\left(b_{1}E_{1}+b_{2}E_{2}+b_{3}E_{3}+b_{4}E_{4}\right)^{m}\gamma_{ijkl}$
$\displaystyle=$
$\displaystyle\sum_{i^{\prime}+j^{\prime}+k^{\prime}+l^{\prime}=m}B_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}^{m}(\boldsymbol{v})E_{1}^{i^{\prime}}E_{2}^{j^{\prime}}E_{3}^{k^{\prime}}E_{4}^{l^{\prime}}\gamma_{ijkl}=\sum_{i^{\prime}+j^{\prime}+k^{\prime}+l^{\prime}=m}\gamma_{i+i^{\prime},j+j^{\prime},k+k^{\prime},l+l^{\prime}}B_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}^{m}(\boldsymbol{v}).$

∎

######  Lemma D.3.

Suppose $p(\boldsymbol{v})$ is a trivariate polynomial with B-form
$p(\boldsymbol{v})=\sum_{i+j+k+l=d}\gamma_{ijkl}B_{ijkl}^{d}(\boldsymbol{v})$.
Then for any $1\leq n\leq d$, the $n$-th order directional derivative of $p$
with respect to the direction
$\boldsymbol{u}=\boldsymbol{v}_{4}-\boldsymbol{v}_{2}$ is given by

$D_{\boldsymbol{u}}^{n}p(\boldsymbol{v})=\frac{d!}{(d-n)!}\sum_{i+j+k+l=d-n}\gamma_{ijkl}^{(n)}B_{ijkl}^{d-n}(\boldsymbol{v})$

with

$\gamma_{ijkl}^{(n)}=\sum_{m=0}^{n}\binom{n}{m}(-1)^{m}\gamma_{i+m,j,k,l+n-m},~{}~{}~{}~{}i+j+k+l=d-n.$

###### Proof.

The directional coordinates of
$\boldsymbol{u}=\boldsymbol{v}_{4}-\boldsymbol{v}_{2}$ is
$\boldsymbol{a}=(-1,0,0,1)$. By (D.2),

$\displaystyle\gamma_{ijkl}^{(n)}=$
$\displaystyle\sum_{i^{\prime}+j^{\prime}+k^{\prime}+l^{\prime}=n}\gamma_{i+i^{\prime},j+j^{\prime},k+k^{\prime},l+l^{\prime}}B_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}^{n}(\boldsymbol{v})=\sum_{i^{\prime}+l^{\prime}=n}\gamma_{i+i^{\prime},j,k,l+l^{\prime}}\frac{n!}{i^{\prime}!l^{\prime}!}(-1)^{i^{\prime}}$
$\displaystyle=$
$\displaystyle\sum_{m=0}^{n}\binom{n}{m}(-1)^{m}\gamma_{i+m,j,k,l+n-m},~{}~{}~{}~{}~{}~{}~{}~{}i+j+k+l=d-n.$

∎

###### Proof of Theorem B.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
$\boldsymbol{\widetilde{u}}$ along $F$, where the directional coordinates are
$(0,1-\widetilde{b}_{3}-\widetilde{b}_{4},\widetilde{b}_{3},\widetilde{b}_{4})$
and
$(0,1-\widetilde{b}_{3}-\widetilde{b}_{4},\widetilde{b}_{4},\widetilde{b}_{3})$
with respect to the tetrahedron $T$ and $\widetilde{T}$, respectively. Thus,
$p$ and $\widetilde{p}$ join continuously along $F$ if and only if

$\sum_{j+k+l=d}\\!\\!\\!\gamma_{0jkl}\frac{d!}{j!k!l!}(1-\widetilde{b}_{3}-\widetilde{b}_{4})^{j}(\widetilde{b}_{3})^{k}(\widetilde{b}_{4})^{l}\\!=\\!\\!\\!\\!\\!\\!\sum_{j+k+l=d}\\!\\!\\!\widetilde{\gamma}_{0jkl}\frac{d!}{j!k!l!}(1-\widetilde{b}_{3}-\widetilde{b}_{4})^{j}(\widetilde{b}_{4})^{k}(\widetilde{b}_{3})^{l},$

That is, $\gamma_{0jkl}=\widetilde{\gamma}_{0jlk}$, with $j+k+l=d$, which
matches our conclusion in (B.8).

ii) Then we consider for $r>0$. First note that

$D_{\boldsymbol{u}}^{n}p(\boldsymbol{v})=D_{\boldsymbol{u}}^{n}\widetilde{p}(\boldsymbol{\widetilde{v}})$
(D.3)

holds for any $\boldsymbol{v}\in F$, $n=0,\ldots,r$, if and only if (D.3)
holds for the direction
$\boldsymbol{u}=\boldsymbol{v}_{5}-\boldsymbol{v}_{3}$. It is a fact because
by the argument in i), all derivatives of $p$ and $\widetilde{p}$
corresponding to the directions $\boldsymbol{v}_{3}-\boldsymbol{v}_{1}$ and
$\boldsymbol{v}_{4}-\boldsymbol{v}_{1}$, agree at every point on $F$. And
derivatives in all other directions can be written as linear combinations of
$D_{\boldsymbol{u}}$, $D_{\boldsymbol{v}_{3}-\boldsymbol{v}_{1}}$ and
$D_{\boldsymbol{v}_{4}-\boldsymbol{v}_{1}}$.

Let $\boldsymbol{b}=(b_{1},b_{2},b_{3},b_{4})$ be the barycentric coordinates
of $\boldsymbol{v}_{5}$ relative to the tetrahedron $T$. Correspondingly, the
directional coordinates of $\boldsymbol{u}$ are
$\boldsymbol{a}=(b_{1},b_{2},b_{3}-1,b_{4})$ and
$\boldsymbol{\widetilde{a}}=(1,0,0,-1)$ with respect to $T$ and
$\widetilde{T}$, respectively.

By Theorem B.1, for each $0\leq n\leq r$,

$\displaystyle D_{\boldsymbol{u}}^{n}p(\boldsymbol{v})|_{F}=$
$\displaystyle\frac{d!}{(d-n)!}\sum_{j+k+l=d-n}\gamma_{0jkl}^{(n)}(\boldsymbol{a})B_{0jkl}^{d-n}(\boldsymbol{v}),$
$\displaystyle D_{\boldsymbol{u}}^{n}\widetilde{p}(\boldsymbol{v})|_{F}=$
$\displaystyle\frac{d!}{(d-n)!}\sum_{j+k+l=d-n}\gamma_{0jkl}^{(n)}(\boldsymbol{\widetilde{a}})B_{0jkl}^{d-n}(\boldsymbol{v}).$

Since for points $\boldsymbol{v}\in F$,
$\widetilde{B}_{0jkl}^{d-n}(\boldsymbol{v})=B_{0jlk}^{d-n}(\boldsymbol{v})$,
it follows that (D.3) holds if and only if for $j+k+l=d-n$, $n=0,\ldots,r$,

$\widetilde{\gamma}_{0jkl}^{(n)}(\boldsymbol{\widetilde{a}})=\gamma_{0jlk}^{(n)}(\boldsymbol{a}).$
(D.4)

By Lemma D.3, for $j+k+l=d-n$,

$\widetilde{\gamma}_{0jkl}^{(n)}(\boldsymbol{\widetilde{a}})=\sum_{m=0}^{n}(-1)^{n-m}\binom{n}{m}\widetilde{\gamma}_{m,j,k,d-m-
j-k}$ (D.5)

In another direction, following the proof of Lemma D.2, for $j+l+k=d-n$,

$\displaystyle\gamma_{0jlk}^{(n)}(\boldsymbol{a})=$
$\displaystyle\left\\{b_{1}E_{1}+b_{2}E_{2}+(b_{3}-1)E_{3}+b_{4}E_{4}\right\\}^{n}\gamma_{0jlk}$
$\displaystyle=$
$\displaystyle\left(b_{1}E_{1}+b_{2}E_{2}+b_{3}E_{3}+b_{4}E_{4}-E_{3}\right)^{n}\gamma_{0jlk}$
$\displaystyle=$
$\displaystyle\sum_{m=0}^{n}(-1)^{n-m}\binom{n}{m}\left(b_{1}E_{1}+b_{2}E_{2}+b_{3}E_{3}+b_{4}E_{4}\right)^{n}\gamma_{0,j,l+n-m,k}$
$\displaystyle=$
$\displaystyle\sum_{m=0}^{n}(-1)^{n-m}\binom{n}{m}\gamma_{0,j,d-j-
k-m,k}^{(m)}(\boldsymbol{b}).$ (D.6)

Combining (D.5) and (D.6), (D.4) holds if and only if for $j+k+l=d-n$ and
$n=0,\ldots,r$,
$\widetilde{\gamma}_{njkl}=\gamma_{0jlk}^{(n)}(\boldsymbol{b})$. By Lemma D.2,

$\gamma_{0jlk}^{(n)}(\boldsymbol{b})=\sum_{i^{\prime}+j^{\prime}+k^{\prime}+l^{\prime}=n}\gamma_{i^{\prime},j+j^{\prime},k+k^{\prime},l+l^{\prime}}B_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}^{n}(\boldsymbol{v}_{5}).$

Thus, (D.3) holds if and only if

$\widetilde{\gamma}_{njkl}=\sum_{i^{\prime}+j^{\prime}+k^{\prime}+l^{\prime}=n}\gamma_{i^{\prime},j+j^{\prime},k+k^{\prime},l+l^{\prime}}B_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}^{n}(\boldsymbol{v}_{5}),$

that is, (B.8) follows. ∎

## References

  * Chai et al. (2020) Chai, X., Gu, H., Li, F., Duan, H., Hu, X., and Lin, K. (2020), “Deep learning for irregularly and regularly missing data reconstruction,” Scientific Reports, 10, 1–18.
  * Cheng et al. (2005) Cheng, S. W., Dey, T. K., Ramos, E. A., and Ray, T. (2005), “Quality meshing of polyhedra with small angles,” International Journal of Computational Geometry & Applications, 15, 421–461.
  * Chew (1989) Chew, L. P. (1989), “Constrained delaunay triangulations,” Algorithmica, 4, 97–108.
  * Chung et al. (2018) Chung, M. K., Wang, Y., and Wu, G. (2018), “Discrete Heat Kernel Smoothing in Irregular Image Domains,” in 2018 40th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), pp. 5101–5104.
  * Duchon (1976) Duchon, J. (1976), “Fonctions-spline à énergie invariante par rotation,” Tech. Rep. 27, Université de Grenoble.
  * Fang and Boas (2009) Fang, Q. and Boas, D. A. (2009), “Tetrahedral mesh generation from volumetric binary and grayscale images,” in 2009 IEEE International Symposium on Biomedical Imaging: From Nano to Macro, pp. 1142–1145.
  * George et al. (1991) George, P. L., Hecht, F., and Saltel, E. (1991), “Automatic mesh generator with specified boundary,” Computer Methods in Applied Mechanics and Engineering, 92, 269–288.
  * Goldsmith et al. (2014) Goldsmith, J., Huang, L., and Crainiceanu, C. M. (2014), “Smooth scalar-on-image regression via spatial Bayesian variable selection,” Journal of Computational and Graphical Statistics, 23, 46–64.
  * Green and Silverman (1994) Green, P. J. and Silverman, B. W. (1994), Nonparametric regression and generalized linear models: a roughness penalty approach, Chapman and Hall, London.
  * Guo et al. (2010) Guo, J., Duan, X., and Shum, C. (2010), “Non-isotropic Gaussian smoothing and leakage reduction for determining mass changes over land and ocean using GRACE data,” Geophysical Journal International, 181, 290–302.
  * Huo et al. (2020) Huo, S., Liu, G., Zhang, J., and Song, C. (2020), “A smoothed finite element method for octree-based polyhedral meshes with large number of hanging nodes and irregular elements,” Computer Methods in Applied Mechanics and Engineering, 359, 112646.
  * Jamin et al. (2015) Jamin, C., Alliez, P., Yvinec, M., and Boissonnat, J.-D. (2015), “CGALmesh: A Generic Framework for Delaunay Mesh Generation,” ACM Transactions on Mathematical Software (TOMS), 41, 1–24.
  * Kang (2020) Kang, Y. (2020), “Consistent blind image deblurring using jump-preserving extrapolation,” Journal of Computational and Graphical Statistics, 29, 372–382.
  * Kim and Wang (2021) Kim, M. and Wang, L. (2021), “Generalized spatially varying coefficient models,” Journal of Computational and Graphical Statistics, 30, 1–10.
  * Lai (1989) Lai, M. J. (1989), “On construction of bivariate and trivariate vertex splines on arbitrary mixed grid partitions,” PhD disseertation, Texas A&M University.
  * Lai and Schumaker (2007) Lai, M. J. and Schumaker, L. L. (2007), Spline functions on triangulations, Cambridge University Press.
  * Lai and Wang (2013) Lai, M. J. and Wang, L. (2013), “Bivariate penalized splines for regression,” Statistica Sinica, 23, 1399–1417.
  * Lee and Schachter (1980) Lee, D.-T. and Schachter, B. J. (1980), “Two algorithms for constructing a Delaunay triangulation,” International Journal of Computer & Information Sciences, 9, 219–242.
  * Lindgren et al. (2011) Lindgren, F., Rue, H., and Lindström, J. (2011), “An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach,” Journal of the Royal Statistical Society, Series B, 73, 423–498.
  * Morris and Carroll (2006) Morris, J. S. and Carroll, R. J. (2006), “Wavelet-based functional mixed models,” Journal of the Royal Statistical Society, Series B, 68, 179–199.
  * Persson and Strang (2004) Persson, P. O. and Strang, G. (2004), “A simple mesh generator in MATLAB,” SIAM Review, 46, 329–345.
  * Ramsay (2002) Ramsay, T. (2002), “Spline smoothing over difficult regions,” Journal of the Royal Statistical Society, Series B, 64, 307–319.
  * Reiss and Ogden (2010) Reiss, P. T. and Ogden, R. T. (2010), “Functional generalized linear models with images as predictors,” Biometrics, 66, 61–69.
  * Roberts et al. (2017) Roberts, D. R., Bahn, V., Ciuti, S., Boyce, M. S., Elith, J., Guillera-Arroita, G., Hauenstein, S., Lahoz-Monfort, J. J., Schröder, B., Thuiller, W., et al. (2017), “Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure,” Ecography, 40, 913–929.
  * Ruppert (1995) Ruppert, J. (1995), “A Delaunay refinement algorithm for quality 2-dimensional mesh generation,” Journal of Algorithms, 18, 548–585.
  * Shewchuk (1996) Shewchuk, J. R. (1996), “Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator,” in Applied Computational Geometry: Towards Geometric Engineering, eds. Lin, M. C. and Manocha, D., Springer, Berlin, Heidelberg, vol. 1148 of Lecture Notes in Computer Science, pp. 203–222.
  * Shewchuk (1998) — (1998), “Tetrahedral mesh generation by Delaunay refinement,” in Proceedings of the Fourteenth Annual ACM Symposium on Computational Geometry, pp. 86–95.
  * Shewchuk and Si (2014) Shewchuk, J. R. and Si, H. (2014), “Higher-quality tetrahedral mesh generation for domains with small angles by constrained delaunay refinement,” in Proceedings of the Thirtieth Annual Symposium on Computational Geometry, pp. 290–299.
  * Si (2015) Si, H. (2015), “TetGen, a Delaunay-based quality tetrahedral mesh generator,” ACM Transactions on Mathematical Software (TOMS), 41, 1–36.
  * Stone (1994) Stone, C. J. (1994), “The use of polynomial splines and their tensor products in multivariate function estimation,” The Annals of Statistics, 22, 118–171.
  * The CGAL Project (2020) The CGAL Project (2020), CGAL User and Reference Manual, CGAL Editorial Board, 5.1.1 ed.
  * Valavi et al. (2019) Valavi, R., Elith, J., Lahoz-Monfort, J. J., and Guillera-Arroita, G. (2019), “blockCV: An R package for generating spatially or environmentally separated folds for k-fold cross-validation of species distribution models,” Methods in Ecology and Evolution, 10, 225–232.
  * Wang and Ranalli (2007) Wang, H. and Ranalli, M. G. (2007), “Low-rank smoothing splines on complicated domains,” Biometrics, 63, 209–217.
  * Wilhelm et al. (2016) Wilhelm, M., Dedè, L., Sangalli, L. M., and Wilhelm, P. (2016), “IGS: an IsoGeometric approach for Smoothing on surfaces,” Computer Methods in Applied Mechanics and Engineering, 302, 70–89.
  * Wong et al. (2016) Wong, R. K., Lee, T. C., Paul, D., and Peng, J. (2016), “Fiber direction estimation, smoothing and tracking in diffusion MRI,” The Annals of Applied Statistics, 10, 1137–1156.
  * Wood (2003) Wood, S. N. (2003), “Thin plate regression splines,” Journal of the Royal Statistical Society, Series B, 65, 95–114.
  * Yue and Speckman (2010) Yue, Y. and Speckman, P. L. (2010), “Nonstationary spatial Gaussian Markov random fields,” Journal of Computational and Graphical Statistics, 19, 96–116.
  * Zhu et al. (2014) Zhu, H., Fan, J., and Kong, L. (2014), “Spatially varying coefficient model for neuroimaging data with jump discontinuities,” Journal of the American Statistical Association, 109, 1084–1098.