Randomness of Shapes and Statistical Inference on Shapes via the Smooth Euler Characteristic Transform

Abstract

In this article, we establish the mathematical foundations for modeling the randomness of shapes and conducting statistical inference on shapes using the smooth Euler characteristic transform. Based on these foundations, we propose two Chi-squared statistic-based algorithms for testing hypotheses on random shapes. Simulation studies are presented to validate our mathematical derivations and to compare our algorithms with state-of-the-art methods to demonstrate the utility of our proposed framework. As real applications, we analyze a dataset of mandibular molars from four genera of primates and show that our algorithms have the power to detect significant shape differences that recapitulate known morphological variation across suborders. Altogether, our discussions bridge the following fields: algebraic and computational topology, probability theory and stochastic processes, Sobolev spaces and functional analysis, analysis of variance for functional data, and geometric morphometrics. Supplementary materials for this article are available online, including a standardized description of the materials available for reproducing the work.

Keywords:

• Functional data analysis
• Karhunen–Loève expansion
• o-minimal structures
• Persistence diagrams
• Reproducing kernel Hilbert spaces

1 Introduction

The quantification of shapes has become an important research direction. It has brought advances to many fields including geometric morphometrics (Boyer et al. 2011; Gao, Kovalsky, and Daubechies 2019; Gao et al. 2019), biophysics and structural biology (Wang et al. 2021; Tang et al. 2022), and radiogenomics (Crawford et al. 2020). When shapes are considered as random variables, their corresponding quantitative summaries are also random, implying that such summaries of random shapes are statistics. The statistical inference on shapes based on these quantitative summaries has been of particular interest (Fasy et al. 2014; Roycraft, Krebs, and Polonik 2023).

In this article, we bring together mathematical and statistical approaches to make three significant contributions to shape statistics: (i) we provide mathematical foundations for the randomness of shapes encountered in applications, bridging algebraic topology (Hatcher 2002) and stochastic processes (Hairer 2009); (ii) we connect the statistical inference on shape-valued data to the well-studied analysis of variance for functional data (fdANOVA, Zhang 2013), bridging topological data analysis (TDA, Edelsbrunner and Harer 2010) and functional data analysis (FDA, Hsing and Eubank 2015); and (iii) our framework does not rely on any assumptions about diffeomorphisms or pre-specified landmarks.

1.1 A Motivating Scientific Question

Through modeling the randomness of shapes, we aim to address the following statistical inference question: Is the observed difference between two groups of shapes statistically significant? For example, the mandibular molars in Figure 1 are from four genera of primates. A pertinent question from a morphological perspective is: In Figure 1, do the molars from genus Tarsius exhibit significant differences from those from the other genera?

Fig. 1 Left: Molars from two suborders of the primates: Haplorhini and Strepsirrhini. The Haplorhini suborder has genera Tarsius (yellow) and Saimiri (grey). The Strepsirrhini suborder has genera Microcebus (blue) and Mirza (green). Right: Relationship between the four primate genera. Tarsier molars exhibit additional high cusps (highlighted in red). A similar figure was published in Wang et al. (2021).

The primary objective of this article is to propose a powerful approach for testing hypotheses on random shapes. This would help address morphology-motivated statistical inference questions like the one raised above. In achieving this objective, we lay down the mathematical foundations that justify our hypothesis testing methods. We take two key steps: In Step 1, we find the appropriate representations of shapes; and in Step 2, we test hypotheses on shapes using these representations. In Section 1.2, we provide a literature review on shape representations and introduce the topological summary employed in this article. Section 1.3 begins by presenting the main theme of our hypothesis testing approach, followed by an overview of our contributions. Since the molars in Figure 1 are diffeomorphic to the two-dimensional unit sphere, some existing diffeomorphism-related methods can be considered for representing the molars (e.g., parameterized surfaces; Kurtek et al. 2011). In contrast, we aim to propose an approach that does not rely on any diffeomorphic assumptions, allowing for a wider range of applications.

1.2 Overview of Shape and Topological Data Analysis

In classical geometric morphometrics, shapes are represented using prespecified points called landmarks (Kendall 1989). The manual landmarking of a collection of shapes requires domain knowledge, can be very labor intensive, and is subject to bias (Boyer et al. 2011). Furthermore, an equal number of landmarks must be selected for each shape in a study in order to make comparisons (e.g., the Procrustes framework discussed in sec. 2.1 of Gao et al. 2019). This necessitates comprehensive information about entire collections of shapes for consistency, which can be difficult to obtain (e.g., landmarking cancer tumors, which can have very different morphology across a population of patients). Unfortunately, many datasets do not come with prespecified landmarks (e.g., Goswami 2015). Although many algorithms can automatically sample reasonable landmarks on shapes when their parameters are fine-tuned (e.g., Gao et al. 2019; Gao, Kovalsky, and Daubechies 2019), using a finite number of landmarks extracted from a continuum inevitably results in the loss of information. Diffeomorphism-based approaches (Dupuis, Grenander, and Miller 1998; Gao et al. 2019) are part of the “computational anatomy” that was historically studied by the “pattern theory school” pioneered by Ulf Grenander (Grenander and Miller 1998). They enable the comparison of (dis-)similarity between shapes with benefit of bypassing the need for landmarks. However, these approaches are based on the assumption that the shapes being compared are diffeomorphic to one another, making them unsuitable for many datasets (e.g., fruit fly wings in Miller 2015). Furthermore, parameterized curves and surfaces (PCS) provide a toolbox for assessing the heterogeneity of shapes with summary statistics that are invariant to reparameterizations (Kurtek et al. 2010, 2011, 2012). Despite their effectiveness in analyzing real data (e.g., DT-MRI brain fibers; Kurtek et al. 2012), PCS are based on assumptions about the diffeomorphism types of the shapes of interest. For example, Kurtek et al. (2011) focuses on surfaces that are diffeomorphic to the two-dimensional unit sphere.

TDA opens the door for landmark-free and diffeomorphism-free representations of shapes. Motivated by differential topology, Turner, Mukherjee, and Boyer (2014) proposed the persistent homology transform (PHT) with the capability to sufficiently encode all information within shapes (Ghrist, Levanger, and Mai 2018). To describe the PHT, we briefly provide some basics of TDA. One common statistical invariant in TDA is the persistence diagram (PD, Edelsbrunner and Harer 2010). When equipped with the Wasserstein distance, the collection of PDs, denoted as $\mathcal{D}$, is a Polish space (Mileyko, Mukherjee, and Harer 2011). Thus, probability measures can be applied, and the randomness of shapes can be represented using the probability measures on $\mathcal{D}$. The PHT takes values in $C({\mathbb{S}}^{d-1};{\mathcal{D}}^d)=\{\text{continuous mapsF}:{\mathbb{S}}^{\text{d‐1}}\to {\mathcal{D}}^{\mathrm{d}}\}$, where ${\mathbb{S}}^{d-1}$ denotes the sphere $\{x\in {\mathbb{R}}^d:|\left|x\right||=1\}$ and ${\mathcal{D}}^d$ is the $d$-fold Cartesian product of $\mathcal{D}$ (Turner, Mukherjee, and Boyer 2014, Lemma 2.1 and Definition 2.1). A single PD does not preserve all information of a shape (Crawford et al. 2020). In contrast, the PHT is injective, which means it preserves all the information of a shape. However, since $\mathcal{D}$ is not a vector space and the distances on $\mathcal{D}$ are abstract (e.g., the Wasserstein and bottleneck distances, Cohen-Steiner et al. 2007), many fundamental statistical concepts do not easily apply to summaries resulting from the PHT. For example, the definition of moments corresponding to probability measures on $\mathcal{D}$ (e.g., means) is highly nontrivial (Mileyko, Mukherjee, and Harer 2011). The difficulty in defining these concepts hinders the application of PHT-based statistical methods in $C({\mathbb{S}}^{d-1};{\mathcal{D}}^d)$.

The smooth Euler characteristic transform (SECT, Crawford et al. 2020) offers an alternative summary statistic for shapes. The SECT not only preserves the information of shapes (Ghrist, Levanger, and Mai 2018, Corollary 1) but also represents shapes using continuous functions instead of PDs. More precisely, the values of the SECT are maps from the sphere ${\mathbb{S}}^{d-1}$ to a separable Banach space $\mathcal{B}\stackrel{\text{def}}{=}C([0,T])$, the collection of continuous functions on a compact interval $[0,T]$ (values of T will be given in (3.1)). Hence, for any shape K, its SECT, denoted as ${\{\text{SECT}(K)(\nu )\}}_{\nu \in {\mathbb{S}}^{d-1}}$, lies in ${\mathcal{B}}^{{\mathbb{S}}^{d-1}}=\{\text{maps }F:{\mathbb{S}}^{d-1}\to \mathcal{B}\}$. Specifically, $\text{SECT}(K)(\nu )$ belongs to $\mathcal{B}$ for each $\nu \in {\mathbb{S}}^{d-1}$. As a result, the randomness of shapes K is represented via the SECT by a collection of $\mathcal{B}$-valued random variables. Probability theory in separable Banach spaces is better developed than in $\mathcal{D}$ (e.g., Hairer 2009). In particular, a $\mathcal{B}$-valued random variable is a stochastic process with its sample paths in $\mathcal{B}$. As we will demonstrate in Section 3, $\mathcal{B}$ here can be replaced with a reproducing kernel Hilbert space (RKHS). The theory of stochastic processes has evolved over a century and FDA is a well-developed branch of statistics. Consequently, a myriad of tools are available to underpin both the randomness of shapes and the statistical inference on shapes.

From an application perspective, Crawford et al. (2020) applied the SECT to magnetic resonance images taken from tumors in a cohort of glioblastoma multiforme (GBM) patients. Using summary statistics derived from the SECT as predictors within Gaussian process regression, the authors demonstrated that the SECT can predict clinical outcomes more effectively than existing tumor shape quantification approaches and common molecular assays. The relative performance of the SECT in the GBM study suggests a promising future for its utility in medical imaging and broader statistical applications related to shape analyses. Similarly, Wang et al. (2021) used derivatives of the Euler characteristic transform (ECT) as predictors in statistical models for subimage analysis. This analysis is akin to variable selection, aiming to identify physical features that are important for distinguishing between two classes of shapes. Lastly, Marsh et al. (2022) highlighted that the SECT outperforms the standard measures employed in organoid morphology.

1.3 Overview of Contributions and Article Organization

Our goal is to address the hypothesis testing question posed in Section 1.1 by employing a landmark-free and diffeomorphism-free approach, which opens up possibilities for further applications in the future. We formulate the question more generically here. Let ${\mathbb{P}}^{(1)}$ and ${\mathbb{P}}^{(2)}$ be two distributions that generate two collections of random shapes, ${\left\{K_i^{(1)}\right\}}_{i=1}^n$ and ${\left\{K_i^{(2)}\right\}}_{i=1}^n$. Detecting whether there is a significant difference between ${\left\{K_i^{(1)}\right\}}_{i=1}^n$ and ${\left\{K_i^{(2)}\right\}}_{i=1}^n$ is equivalent to rejecting the hypothesis ${\mathbb{P}}^{(1)}={\mathbb{P}}^{(2)}$. Since each shape $K_i^{(j)}$ is random, $\text{SECT}(K_i^{(j)})$ is a random variable taking values in a vector space (as discussed in Section 1.2) and can be decomposed as follows (see Theorem 5.1 for a rigorous version) (1.1) $$\text{SECT}(K_i^{(j)})=m^{(j)}+\text{random terms},\text{ for j}\in \{1,2\},$$(1.1) where $m^{(j)}$ denotes the mean of $\text{SECT}(K_i^{(j)})$ with respect to the distribution ${\mathbb{P}}^{(j)}$. The random terms in (1.1) can be characterized by the Karhunen–Loève (KL) expansion (Hsing and Eubank 2015, sec. 7.3). To reject ${\mathbb{P}}^{(1)}={\mathbb{P}}^{(2)}$, it suffices to reject $m^{(1)}=m^{(2)}$. That is, the question posed in Section 1.1 can be addressed by testing for the equality of two means. The important component of the test is the variance represented by the random terms in (1.1). In Section 5, we formulate this test as a two-sample problem in the fdANOVA literature (Zhang 2013, sec. 5.2). In addition, using the KL expansion, we provide a ${\chi }^2$-statistic in Section 5 to test the hypothesis. Throughout the article, our focus is on the two-sample problem. However, one may also consider employing the one-way fdANOVA to compare the means of three or more groups of shapes. The theoretical foundation and numerical experiments for this aspect are left for future research.

To develop our framework, we have to address the following mathematical foundation related questions: (i) What underlying probability spaces allow the randomness of shapes and their corresponding SECT? and (ii) Are the conditions of the KL expansion satisfied in our setting? We answer these questions in Sections 3 and 4—we model the randomness of shapes via the SECT using RKHS-valued random fields. The “theory of random sets” is a well-established framework for characterizing set-valued random variables (Molchanov 2005). However, its application to persistent homology-based statistics (e.g., the SECT) remains underexplored. In this article, we introduce a new probability space to characterize the randomness of shapes in a manner compatible with the SECT.

We first propose a collection of shapes as our sample space on which the SECT is well-defined. We then demonstrate that every shape in this collection has its SECT in $C({\mathbb{S}}^{d-1};\mathcal{H})=\{\text{continuous mapsF}:{\mathbb{S}}^{\text{d‐1}}\to \mathcal{H}\}$, where $\mathcal{H}=H_0^1([0,T])$ is not only a Sobolev space (Brezis 2011) but also an RKHS (reasons for using $[0,T]$ instead of $(0,T)$ for $H_0^1([0,T])$ are in Appendix A.1). Importantly, $C({\mathbb{S}}^{d-1};\mathcal{H})$ is a separable Banach space (Theorem C.1) and, hence, a Polish space. It helps construct a probability space to characterize the distributions of shapes. Building on this probability space, we define the mean and covariance of the SECT. Using the Sobolev embedding, we present some properties of the mean and covariance, which pave the way for the KL expansion of the SECT.

Traditionally, the statistical inference on shapes in TDA is conducted in the persistence diagram space $\mathcal{D}$, which is unsuitable for exponential family-based distributions and requires any corresponding statistical inference to be highly nonparametric (Fasy et al. 2014; Robinson and Turner 2017). The PHT-based statistical inference in $C({\mathbb{S}}^{d-1};{\mathcal{D}}^d)$ is even more difficult. With the KL expansion of the SECT, we propose a ${\chi }^2$-statistic for testing hypotheses on shapes. Beyond the mathematical foundations, we also provide simulation studies to illustrate the performance of our proposed hypothesis testing method. Lastly, we apply our proposed framework to answer the motivating question raised in Section 1.1.

We organize this article as follows. In Section 2, we provide the mathematical preparations. In Section 3, we define the SECT for a specific collection of shapes, highlighting its properties. In Section 4, we construct a probability space to model shape distributions. In Section 5, we propose the KL expansion of the SECT, leading to a statistic for hypothesis testing. In Section 6, we conduct simulation studies to evaluate our method. In Section 7, we apply our method to real data. In Section 8, we conclude the article. The Appendix provides the proofs of theorems, further data analysis, and future research topics.

2 Notations and Mathematical Preparations

To model the shapes discussed in our motivating question from Section 1.1, we need certain preparations regarding (i) topology and (ii) function spaces.

Topology

The first question we must address is: What are the “shapes” in our framework? Ghrist, Levanger, and Mai (2018) and Curry, Mukherjee, and Turner (2022) applied o-minimal structures (van den Dries 1998) to prove the injectivity of the PHT. Subsequent to this, o-minimal structures have been applied in many TDA studies to model shapes (e.g., Jiang, Kurtek, and Needham 2020; Kirveslahti and Mukherjee 2023). To stay consistent with the existing literature, we also model shapes using o-minimal structures. An o-minimal structure is a sequence $\mathcal{S}={\left\{{\mathcal{S}}_n\right\}}_{n\ge 1}$ of subset collections ${\mathcal{S}}_n\subseteq 2^{{\mathbb{R}}^n}$ satisfying six set-theoretical axioms, where $2^{{\mathbb{R}}^n}$ denotes the power set of ${\mathbb{R}}^n$. The precise definition of o-minimal structures is available in van den Dries (1998) and is provided in Appendix A.3 for the reader’s convenience.

A typical example of o-minimal structures is the collection of semialgebraic sets. Specifically, a set $K\subseteq {\mathbb{R}}^n$ is semialgebraic if it can be expressed as a finite union of sets of the form $\{x\in {\mathbb{R}}^n|p(x)=0,q_1(x)>0,\dots ,q_k(x)>0\}$, where $p,q_1,\dots ,q_k$ are polynomial functions on ${\mathbb{R}}^n$. If we define ${\mathcal{S}}_n$ as the collection of semialgebraic subsets of ${\mathbb{R}}^n$, then $\mathcal{S}={\left\{{\mathcal{S}}_n\right\}}_{n\ge 1}$ is an o-minimal structure (van den Dries 1998, chap. 2). The unit sphere ${\mathbb{S}}^{d-1}$, open ball $B(0,R)=\{x\in {\mathbb{R}}^d|\left|\right|x|{|}^2<R^2\}$ for any R > 0, and all polyhedra (e.g., polygon meshes in computer graphics) are semialgebraic, given that they can be represented using either the polynomial $\left|\right|x|{|}^2$ or affine functions. We assume the following:

Assumption 1.

The o-minimal structure $\mathcal{S}$ of interest contains all semialgebraic sets.

Hereafter, a “shape” refers to a compact set $K\in \underset{n\ge 1}{\cup }{\mathcal{S}}_n$ for a prespecified o-minimal structure $\mathcal{S}={\left\{{\mathcal{S}}_n\right\}}_{n\ge 1}$ satisfying Assumption 1. Assumption 1 incorporates many common shapes (e.g., balls and polyhedra) in our framework. More importantly, it implies the subsequent Theorem 2.1 through the “triangulation theorem” (van den Dries 1998, chap. 8). Although the definition of an o-minimal structure $\mathcal{S}$ is highly abstract (see Appendix A.3), each compact set in $\mathcal{S}$ resembles a polyhedron, which is precisely stated as follows.

Theorem 2.1.

Suppose $\mathcal{S}={\left\{{\mathcal{S}}_n\right\}}_{n\ge 1}$ is an o-minimal structure satisfying Assumption 1 and $K\in \underset{n\ge 1}{\cup }{\mathcal{S}}_n$. If K is compact, there exists a finite simplicial complex S such that the polyhedron $\left|S\right|\stackrel{\text{def}}{=}\underset{s\in S}{\cup }s$ is homeomorphic to K, where each $s\in S$ denotes a simplex.

Herein, a finite simplicial complex S is a finite collection of simplexes. Each face of a simplex $s\in S$ also belongs to S (i.e., S is a so-called “closed complex” referred to in chap. 8 of van den Dries 1998). Theorem 2.1 directly results from the “triangulation theorem” (van den Dries 1998); hence, its proof is omitted. For the dth component ${\mathcal{S}}_d$ of $\mathcal{S}={\left\{{\mathcal{S}}_n\right\}}_{n\ge 1}$, Theorem 2.1 indicates that the compact sets $K\in {\mathcal{S}}_d$ are subsets of ${\mathbb{R}}^d$ that are homeomorphic to polyhedra. Theorem 2.1 also implies that the homology groups of each compact $K\in {\mathcal{S}}_d$ are well-defined and finitely generated; hence, the Betti numbers and Euler characteristic of K are well-defined and finite (Hatcher 2002, chap. 2).

Function Spaces

We apply the following notations throughout this article:

1. For any normed space $\mathcal{V}$, let $\left|\right|·|{|}_{\mathcal{V}}$ denote its norm. Denote $\left|\right|·|{|}_{{\mathbb{R}}^d}$ as $\left|\right|·\left|\right|$ for succinctness.

2. Let X be a compact metric space equipped with metric d_X , and let $\mathcal{V}$ denote a normed space. $C(X;\mathcal{V})$ is the collection of continuous maps from X to $\mathcal{V}$. Furthermore, $C(X;\mathcal{V})$ is a normed space equipped with $\left|\right|f|{|}_{C(X;\mathcal{V})}=\underset{x\in X}{\mathrm{\text{sup}}}||f(x)|{|}_{\mathcal{V}}$. The Hölder space $C^{0,\frac{1}{2}}(X;\mathcal{V})$ is defined as $\left\{f\in C(X;\mathcal{V})|\underset{x,y\in X,\mathrm{ }x\ne y}{\mathrm{\text{sup}}}(\frac{\left|\right|f(x)-f(y)|{|}_{\mathcal{V}}}{\sqrt{d_X(x,y)}})<\infty \right\}$. Here, $C^{0,\frac{1}{2}}(X;\mathcal{V})$ is a normed space equipped with $\left|\right|f|{|}_{C^{0,\frac{1}{2}}(X;\mathcal{V})}=|\left|f\right|{|}_{C(X;\mathcal{V})}+\underset{x,y\in X,\mathrm{ }x\ne y}{\mathrm{\text{sup}}}(\frac{\left|\right|f(x)-f(y)|{|}_{\mathcal{V}}}{\sqrt{d_X(x,y)}})$. Obviously, $C^{0,\frac{1}{2}}(X;\mathcal{V})\subseteq C(X;\mathcal{V})$. For simplicity, we denote $C(X)=C(X;\mathbb{R})$ and $C^{0,\frac{1}{2}}(X)=C^{0,\frac{1}{2}}(X;\mathbb{R})$. For a given T > 0 (e.g., see (3.1)), we denote $C([0,T])$ as $\mathcal{B}$.

3. The inner product of $\mathcal{H}=H_0^1([0,T])=\{f\in L^2([0,1])|f\prime \in L^2([0,T])\text{ and }f(0)=f(T)=0\}$ is defined as $〈f,g〉={\int }_0^{T}f\prime (t)g\prime (t)\mathit{\text{dt}}$ (Brezis 2011, chap. 8.3, Remark 17).

4. Suppose $(Y,d_Y)$ is a metric space (not necessarily compact). Both $\mathcal{B}(Y)$ and $\mathcal{B}(d_Y)$ denote the Borel algebra generated by the metric topology corresponding to d_Y.

5. ${\{F(z)\}}_{z\in Z}$ denotes a function F defined on the set Z.

The following inequalities are useful for deriving many results presented in this article (2.1) $$\left|\right|f|{|}_{\mathcal{B}}\le |\left|f\right|{|}_{C^{0,\frac{1}{2}}([0,T])}\le {\tilde{C}}_T\left|\right|f|{|}_{\mathcal{H}},\text{ for all f}\in \mathcal{H},$$(2.1) where ${\tilde{C}}_T$ is a constant depending only on T. The first inequality in (2.1) results from the definition of $\left|\right|·|{|}_{C^{0,\frac{1}{2}}([0,T])}$, while the second inequality is from Brezis (2011) (Corollary 9.14; also see Appendix L.2). Equation (2.1)(2.1) $$\left|\right|f|{|}_{\mathcal{B}}\le |\left|f\right|{|}_{C^{0,\frac{1}{2}}([0,T])}\le {\tilde{C}}_T\left|\right|f|{|}_{\mathcal{H}},\text{ for all f}\in \mathcal{H},$$(2.1) implies the following Sobolev embedding (2.2) $$H_0^1([0,T])\stackrel{\text{def}}{=}\mathcal{H}\subseteq C^{0,\frac{1}{2}}([0,T])\subseteq \mathcal{B}\stackrel{\text{def}}{=}C([0,T]).$$(2.2)

3 Smooth Euler Characteristic Transform

In this section, we give the background on the SECT and propose corresponding mathematical foundations. Notably, we specify the “sample space”—a collection of shapes on which the SECT is well-defined. The SECT of the shapes in this sample space has properties that are suitable for the probability theory developed in Section 4. The molars in the motivating question from Section 1.1 will be modeled as elements of the sample space.

Suppose an o-minimal structure $\mathcal{S}={\left\{{\mathcal{S}}_n\right\}}_{n\ge 1}$ satisfying Assumption 1 is given, and we focus on shapes in d-dimensional space ${\mathbb{R}}^d$. We assume the shape $K\in {\mathcal{S}}_d$ is compact and $K\subseteq B(0,R)=\{x\in {\mathbb{R}}^d:|\left|x\right||<R\}$, for example, the $K\subseteq {\mathbb{R}}^2$ in Figure 2 or the surfaces of the mandibular molars in ${\mathbb{R}}^3$ as presented by Figure 1. For each direction $\nu \in {\mathbb{S}}^{d-1}$, we define a filtration ${\left\{K_t^{\nu }\right\}}_{t\in [0,T]}$ of sublevel sets by the following (see Figure 2 for an illustration) (3.1) $$\begin{array}{c}{K}_t^{\nu }\stackrel{\text{def}}{=}\left\{x\in K|x·\nu \le t-R\right\},\text{ for all t}\in [0,\mathrm{T}],\\ \hspace{1em}\text{ where }T\stackrel{\text{def}}{=}2R.\end{array}$$(3.1)

Fig. 2 Consider the two-dimensional shape $K\in {\mathcal{S}}_2$ in the left panel. For each pair of ν and t, the equation $x·\nu =t-R$ represents a straight line (or a hyperplane in a high-dimensional space). The subset $K_t^{\nu }$ denotes the region below this line. Let ${\varphi }_{\nu }(x)=x·\nu +R$, then $K_t^{\nu }=\{x\in K|{\varphi }_{\nu }(x)\le t\}$. The right panel presents the function $(\nu ,t)\mapsto \text{SECT}(K)(\nu ,t)$, where $\nu \in {\mathbb{S}}^1$ is identified by $\theta \in [0,2\pi ]$ through $\nu =(\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}\theta ,\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}\theta )$. Procedures for generating the shape K and the right panel are given in Appendix D.1.

We then have the following Euler characteristic curve (ECC, denoted as ${\chi }_t^{\nu }$) in direction ν (3.2) $$\begin{array}{c}{\chi }_t^{\nu }(K)\stackrel{\text{def}}{=}\text{the Euler characteristic of }{\mathrm{K}}_{\mathrm{t}}^{\nu }\text{=}\chi ({\mathrm{K}}_{\mathrm{t}}^{\nu })\\ =\sum _{k=0}^{d-1}{(-1)}^k·{\beta }_k(K_t^{\nu }),\end{array}$$(3.2) for $t\in [0,T]$, where ${\beta }_k(K_t^{\nu })$ is the kth Betti number of $K_t^{\nu }$. The sum in (3.2) ends at d–1 because higher homology groups are trivial (Curry, Mukherjee, and Turner 2022, sec. 4). If $K_t^{\nu }$ is a triangle mesh, $\chi (K_t^{\nu })=\#V-\#E+\#F$, where $\#V,\#E$, and $\#F$ denote the number of vertices, edges, and faces of the mesh, respectively. Due to Theorem 2.1, the compactness of K guarantees that the Betti numbers in (3.2) are well-defined and finite.

The Euler characteristic transform (ECT) defined as $\text{ECT}(K):{\mathbb{S}}^{d-1}\to {\mathbb{Z}}^{[0,T]},\nu \mapsto {\{{\chi }_t^{\nu }(K)\}}_{t\in [0,T]}$ was proposed by Turner, Mukherjee, and Boyer (2014) as an alternative to the PHT. Based on the ECT, Crawford et al. (2020) further proposed the SECT as follows (3.3) $$\begin{array}{cc}\hfill & \text{SECT}(K):{\mathbb{S}}^{d-1}\to {\mathbb{R}}^{[0,T]},\hfill \\ \hfill & \nu \mapsto \text{SECT}(K)(\nu )={\left\{\text{SECT}(K)(\nu ,t)\right\}}_{t\in [0,T]},\hfill \\ \hfill & \text{ where SECT}(K)(\nu ,t)\stackrel{\text{def}}{=}{\int }_0^t{\chi }_{\tau }^{\nu }(K)d\tau -\frac{t}{T}{\int }_0^T{\chi }_{\tau }^{\nu }(K)d\tau .\hfill \end{array}$$(3.3)

A visualization of the function $(\nu ,t)\mapsto \text{SECT}(K)(\nu ,t)$ is presented in Figure 2. The following lemma implies that the Lebesgue integrals in (3.3) are well-defined.

Lemma 3.1.

For any fixed $K\in {\mathcal{S}}_d$ and $\nu \in {\mathbb{S}}^{d-1}$, the function $t\mapsto \chi (K_t^{\nu })$ is piecewise constant with only finitely many discontinuities.

Through the “cell decomposition theorem” (van den Dries 1998, chap. 3), Lemma 3.1 directly follows from either Lemma 3.4 of Curry, Mukherjee, and Turner (2022) or “(2.10) Proposition” in Chapter 4 of van den Dries (1998). Hence, the proof of Lemma 3.1 is omitted.

To investigate the distribution of $\text{SECT}(K)$ over different shapes K, we introduce the following condition to restrict our attention to a subset of ${\mathcal{S}}_d$.

Condition 3.1.

Let $K\in {\mathcal{S}}_d$. The condition is that K satisfies the following inequality (3.4) $$\underset{k\in \{0,\mathrm{ }\dots ,\mathrm{ }d-1\}}{\mathrm{\text{sup}}}\left[\underset{\nu \in {\mathbb{S}}^{d-1}}{\mathrm{\text{sup}}}\left(\#\{\xi \in {\text{Dgm}}_k(K;{\varphi }_{\nu })|\text{pers}(\xi )>0\}\right)\right]\le \frac{M}{d},$$(3.4) where ${\text{Dgm}}_k(K;{\varphi }_{\nu })$ is the PD of K associated with the function ${\varphi }_{\nu }(x)=x·\nu +R$ (also see Figure 2), $\text{pers}(\xi )$ is the persistence of the homology feature ξ, $\#\{·\}$ denotes the cardinality of a multiset, and M > 0 is a sufficiently large prespecified number.

Condition 3.1 involves technicalities from computational topology (Edelsbrunner and Harer 2010). To maintain the flow of the article, we relegate the details of this condition, as well as the definitions of ${\text{Dgm}}_k(K;{\varphi }_{\nu })$ and $\text{pers}(\xi )$, to Appendix B. Heuristically, Condition 3.1 implies the existence of a uniform upper bound on the number of nontrivial homology features of K across all directions ν. Hereafter, we focus on shapes in the following collection $$\begin{array}{c}{\mathcal{S}}_{R,d}^M\stackrel{\text{def}}{=}\{K\in {\mathcal{S}}_d|K\subseteq B(0,R)\text{ is compact and satisfies}\\ \text{Condition} 3.\text{1with fixed M}>0\}.\end{array}$$

Our proposed collection ${\mathcal{S}}_{R,d}^M$ is suitable for modeling shapes in many applications. For example, the surfaces of the molars in Figure 1 are compact subsets of ${\mathbb{R}}^3$, bounded by a common open ball, and can be approximately represented by triangle meshes (hence, modeled by an o-minimal structure satisfying Assumption 1). In addition, the four genera of primates in Figure 1 share a phylogentic relationship which implies that their molars have common baseline features and satisfy Condition 3.1 with a shared upper bound M. In each application, the dimension d and radius R of the ball B(0, R) can easily be determined based on observed shapes. Although our mathematical framework requires the existence of such an M in (3.4), the value of M is not needed for our statistical methodology (see Section 5). Thus, Condition 3.1 does not hinder our proposed statistical methodology.

Lemma 3.1 implies that the function ${\{{\chi }_t^{\nu }(K)\}}_{t\in [0,T]}$ of t belongs to $L^1([0,T])$. Therefore, the function $\text{SECT}(K)(\nu )={\{\text{SECT}(K)(\nu ,t)\}}_{t\in [0,T]}$ of t is absolutely continuous on $[0,T]$. Furthermore, we have the following regularity result of the Sobolev type.

Lemma 3.2.

For any $K\in {\mathcal{S}}_{R,d}^M$ and $\nu \in {\mathbb{S}}^{d-1}$, the function $\text{SECT}(K)(\nu )$ belongs to $\mathcal{H}$.

Lemma 3.2 is a special case of Lemma C.3. It indicates $\text{SECT}({\mathcal{S}}_{R,d}^M)\subseteq {\mathcal{H}}^{{\mathbb{S}}^{d-1}}=\{\text{all mapsF}:{\mathbb{S}}^{\text{d‐1}}\to \mathcal{H}\}$, which is enhanced by the following result.

Theorem 3.2.

For each $K\in {\mathcal{S}}_{R,d}^M$, we have: (i) There exists a constant $C_{M,R,d}^{*}$ depending only on M, R, and d such that the following inequality holds for any directions ${\nu }_1,{\nu }_2\in {\mathbb{S}}^{d-1}$, (3.5) $$\begin{array}{cc}\hfill & \left||\text{SECT}(K)({\nu }_1)-\text{SECT}(K)({\nu }_2)\right|{|}_{\mathcal{H}}\hfill \\ \hfill & \hspace{1em}\le C_{M,R,d}^{*}·\sqrt{\left|\right|{\nu }_1-{\nu }_2\left|\right|+\left|\right|{\nu }_1-{\nu }_2|{|}^2};\hfill \end{array}$$(3.5) and (ii) $\text{SECT}(K)\in C^{0,\frac{1}{2}}({\mathbb{S}}^{d-1};\mathcal{H})$, where ${\mathbb{S}}^{d-1}$ is equipped with the geodesic distance $d_{{\mathbb{S}}^{d-1}}$.

Results complementary to Theorem 3.2 can be found in Theorem C.3, which imply that the function $(\nu ,t)\mapsto \text{SECT}(K)(\nu ,t)$ belongs to $C^{0,\frac{1}{2}}({\mathbb{S}}^{d-1}\times [0,T])$. Theorem 3.2(i) is an analog of Lemma 2.1 in Turner, Mukherjee, and Boyer (2014). Theorem 3.2(ii) implies $\text{SECT}({\mathcal{S}}_{R,d}^M)\subseteq C^{0,\frac{1}{2}}({\mathbb{S}}^{d-1};\mathcal{H})\subseteq C({\mathbb{S}}^{d-1};\mathcal{H})\subseteq {\mathcal{H}}^{{\mathbb{S}}^{d-1}}$. As a result, (3.3) defines the following map (3.6) $$\text{SECT}:{\mathcal{S}}_{R,d}^M\to C({\mathbb{S}}^{d-1};\mathcal{H}),K\mapsto \text{SECT}(K).$$(3.6)

In Appendix D.1, we provide detailed proof-of-concept examples (similar to Figure 2) to visually illustrate the $\text{SECT}$ and support the regularity results in Theorems 3.2 and C.3.

Corollary 1 of Ghrist, Levanger, and Mai (2018) implies the following result, which shows that the SECT preserves all the information of shapes $K\in {\mathcal{S}}_{R,d}^M$.

Theorem 3.3.

The map $\text{SECT}$ defined in (3.6) is injective for all dimensions d.

The map $\text{SECT}:{\mathcal{S}}_{R,d}^M\to C({\mathbb{S}}^{d-1};\mathcal{H})$ is injective, but not surjective. Specifically, Theorem 3.2 suggests that the image of $\text{SECT}$ does not lie outside of $C^{0,\frac{1}{2}}({\mathbb{S}}^{d-1};\mathcal{H})$. An explicit characterization of the image $\text{SECT}({\mathcal{S}}_{R,d}^M)$ remains a topic for future research.

Inspired by Theorem 3.3, one may consider reconstructing a shape $K$ from either the $\text{SECT}(K)$ or the $\text{ECT}(K)$. From a theoretical standpoint, a shape $K$ can be reconstructed using the “Schapira’s inversion formula” (Schapira 1995). Further details are available in Ghrist, Levanger, and Mai (2018). From an algorithmic perspective, the proof of Theorem 3.1 in Turner, Mukherjee, and Boyer (2014) offers an algorithm to reconstruct low-dimensional meshes from their ECT. Nevertheless, effective algorithmic approaches to reconstructing shapes are still underdeveloped. Challenges in reconstructing shapes are extensively discussed in Fasy et al. (2018). A comprehensive exploration of the reconstruction using SECT is also left for future research.

Together with (3.6), Theorem 3.3 allows us to represent each $K\in {\mathcal{S}}_{R,d}^M$ by $\text{SECT}(K)\in C({\mathbb{S}}^{d-1};\mathcal{H})$. This perspective aids us in modeling the randomness of shapes using probability measures on the separable Banach space $C({\mathbb{S}}^{d-1};\mathcal{H})$. Here, we prefer $C({\mathbb{S}}^{d-1};\mathcal{H})$ over $\frac{1}{2}$-Hölder space $C^{0,\frac{1}{2}}({\mathbb{S}}^{d-1};\mathcal{H})$. This is because $\frac{1}{2}$-Hölder spaces are typically not separable (Hairer 2009, Remark 4.21). The separability condition is essential for probability measures on Banach spaces to exhibit non-pathological behavior (Hairer 2009, sec. 4).

4 Probabilistic Distributions over the SECT

To address the motivating question outlined in Section 1.1 using hypothesis testing, we need to view the observed shapes (e.g., the molars in Figure 1) as shape-valued random variables. In this section, we construct a probability space to model the randomness of shapes and make the SECT a random variable (in the measurable sense) taking values in $C({\mathbb{S}}^{d-1};\mathcal{H})$. More importantly, this probability space helps justify the KL expansion of the SECT, which lays the foundations for our hypothesis testing method in Section 5.

Probability Space

Suppose ${\mathcal{S}}_{R,d}^M$ is equipped with a σ-algebra $\mathcal{F}$. A distribution of shapes K across ${\mathcal{S}}_{R,d}^M$ is represented by a probability measure $\mathbb{P}=\mathbb{P}(\mathit{\text{dK}})$ on $\mathcal{F}$. Then, $({\mathcal{S}}_{R,d}^M,\mathcal{F},\mathbb{P})$ is a probability space. For each fixed $(\nu ,t)$, the integer-valued map ${\chi }_t^{\nu }:K\mapsto {\chi }_t^{\nu }(K)$ is defined on ${\mathcal{S}}_{R,d}^M$. Hereafter, we assume the following:

Assumption 2.

For each fixed $(\nu ,t)\in {\mathbb{S}}^{d-1}\times [0,T]$, the map ${\chi }_t^{\nu }:({\mathcal{S}}_{R,d}^M,\mathcal{F})\to (\mathbb{R},\mathcal{B}(\mathbb{R}))$ is a measurable function and, hence, a real-valued random variable.

A σ-algebra $\mathcal{F}$ satisfying Assumption 2 exists. Here, we construct a metric ρ on ${\mathcal{S}}_{R,d}^M$ and show that the Borel algebra $\mathcal{B}(\rho )$ induced by ρ satisfies Assumption 2. We define (4.1) $$\begin{array}{c}\rho (K_1,K_2)\stackrel{\text{def}}{=}\underset{\nu \in {\mathbb{S}}^{d-1}}{\mathrm{\text{sup}}}\left\{{\left({\int }_0^T| {\chi }_{\tau }^{\nu }(K_1)-{\chi }_{\tau }^{\nu }(K_2){|}^{2}d\tau \right)}^{1/2}\right\},\\ \hspace{1em}\text{for all }{\mathrm{K}}_1,{\mathrm{K}}_2\in {\mathcal{S}}_{\mathrm{R},\mathrm{d}}^{\mathrm{M}}.\end{array}$$(4.1)

Theorem 4.1.

The map ρ is a metric on ${\mathcal{S}}_{R,d}^M$. Assumption 2 is satisfied if $\mathcal{F}=\mathcal{B}(\rho )$.

Under Assumption 2, the ECC ${\left\{{\chi }_t^{\nu }\right\}}_{t\in [0,T]}$, for each $\nu \in {\mathbb{S}}^{d-1}$, is a stochastic process defined on the probability space $({\mathcal{S}}_{R,d}^M,\mathcal{F},\mathbb{P})$. Since each sample path ${\{{\chi }_t^{\nu }(K)\}}_{t\in [0,T]}$ has finitely many discontinuities (Lemma 3.1), ${\int }_0^t{\chi }_{\tau }^{\nu }(K)d\tau$ for each $t\in [0,T]$ is a Riemann integral, which is equal to the limit of Riemann sum ${\int }_0^t{\chi }_{\tau }^{\nu }(K)d\tau =\underset{n\to \infty }{\mathrm{\text{lim}}}\left\{\frac{t}{n}{\sum }_{l=1}^n{\chi }_{\frac{\mathit{\text{tl}}}{n}}^{\nu }(K)\right\}$. Given that each ${\chi }_{\frac{\mathit{\text{tl}}}{n}}^{\nu }$ is a random variable under Assumption 2, the limit of the Riemann sum for each $t\in [0,T]$ is a random variable as well. Therefore, for each $\nu \in {\mathbb{S}}^{d-1},{\left\{{\int }_0^t{\chi }_{\tau }^{\nu }d\tau \right\}}_{t\in [0,T]}$ with ${\int }_0^t{\chi }_{\tau }^{\nu }d\tau :K\mapsto {\int }_0^t{\chi }_{\tau }^{\nu }(K)d\tau$ is a stochastic process. Then, under Assumption 2, (3.3) defines the following stochastic process on $({\mathcal{S}}_{R,d}^M,\mathcal{F},\mathbb{P})$ for each $\nu \in {\mathbb{S}}^{d-1}$ (4.2) $$\begin{array}{cc}\hfill & \text{SECT}(\nu )\stackrel{\text{def}}{=}{\left\{{\int }_0^t{\chi }_{\tau }^{\nu }d\tau -\frac{t}{T}{\int }_0^T{\chi }_{\tau }^{\nu }d\tau \stackrel{\text{def}}{=}\text{SECT}(\nu ,t)\right\}}_{t\in [0,T]}.\hfill \end{array}$$(4.2)

Precisely, for each fixed ν, we have the stochastic process $\text{SECT}(\nu ):K\mapsto \text{SECT}(K)(\nu )={\{\text{SECT}(K)(\nu ,t)\}}_{t\in [0,T]}$ defined on $({\mathcal{S}}_{R,d}^M,\mathcal{F},\mathbb{P})$; and, for each fixed $(\nu ,t)$, we have the real-valued random variable $\text{SECT}(\nu ,t):K\mapsto \text{SECT}(K)(\nu ,t)$ defined on $({\mathcal{S}}_{R,d}^M,\mathcal{F},\mathbb{P})$.

Since $\mathcal{H}$ is an RKHS (Appendix A.1), Lemma 3.2 and Theorem 4.2, together with Theorem 7.1.2 of Hsing and Eubank (2015), imply the following. Its proof is omitted.

Theorem 4.2.

(i) For each $\nu \in {\mathbb{S}}^{d-1},\text{ SECT}(\nu )$ is a real-valued stochastic process with sample paths in $\mathcal{H}$. Equivalently, $\text{SECT}(\nu )$ is a random variable taking values in $(\mathcal{H},\mathcal{B}(\mathcal{H}))$. (ii) The map $\text{SECT}$ defined in (3.6) is a random variable taking values in $C({\mathbb{S}}^{d-1};\mathcal{H})$.

Using Theorem 4.2 in conjunction with Theorem 3.3, we can represent random shapes (which model the surfaces of the mandibular molars in Figure 1) as $C({\mathbb{S}}^{d-1},\mathcal{H})$-valued random variables. This representation through the SECT has no loss of information.

In Appendix D.2, we provide proof-of-concept examples to illustrate random shapes and their SECT representations visually. These examples relate the SECT to Fréchet regression (Petersen and Müller 2019), Wasserstein regression (Chen et al. 2023), and manifold learning (Dunson and Wu 2021; Meng and Eloyan 2021; Li, Mukhopadhyay, and Dunson 2022).

Mean and Covariance of the SECT

For deriving the KL expansion in Section 5, we define the mean and covariance of the SECT. To do so, we need the following lemma.

Lemma 4.1.

For any probability measure $\mathbb{P}$ defined on the measurable space $({\mathcal{S}}_{R,d}^M,\mathcal{F})$, we have $\mathbb{E}\{\underset{\nu \in {\mathbb{S}}^{d-1}}{\mathrm{\text{sup}}}\left|\right|\text{SECT}(\nu )|{|}_{\mathcal{H}}^2\}={\int }_{{\mathcal{S}}_{R,d}^M}\{\underset{\nu \in {\mathbb{S}}^{d-1}}{\mathrm{\text{sup}}}\left|\right|\text{SECT}(K)(\nu )|{|}_{\mathcal{H}}^2\} \mathbb{P}(\mathit{\text{dK}})<\infty .$

Lemma 4.1, together with (2.1), implies that $\mathbb{E}|\text{SECT}(\nu ,t){|}^2\le {\tilde{C}}_T^2·\mathbb{E}||\text{SECT}(\nu )|{|}_{\mathcal{H}}^2<\infty$ for all $(\nu ,t)\in {\mathbb{S}}^{d-1}\times [0,T]$. Then, we define the mean and covariance functions as follows (4.3) $$\begin{array}{cc}\hfill & m_{\nu }(t)=\mathbb{E}\left\{\text{SECT}(\nu ,t)\right\}={\int }_{{\mathcal{S}}_{R,d}^M}\text{SECT}(K)(\nu ,t)\mathbb{P}(\mathit{\text{dK}}),\hfill \\ \hfill & {\Xi }_{\nu }(s,t)=\text{cov}(\text{SECT}(\nu ,s),\text{SECT}(\nu ,t)),\text{ for }s,t\in [0,T]\hfill \\ \hfill & \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{ and }\nu \in {\mathbb{S}}^{d-1}.\hfill \end{array}$$(4.3)

Lemma C.4 provides several properties of the mean $m_{\nu }$ and covariance ${\Xi }_{\nu }$ that validate our KL expansion of $\text{SECT}(\nu )$ in Section 5. Additionally, Lemma C.4 demonstrates that the mean $\mathit{m}\stackrel{\text{def}}{=}{\left\{m_{\nu }\right\}}_{\nu \in {\mathbb{S}}^{d-1}}$ of SECT belongs to $C({\mathbb{S}}^{d-1};\mathcal{H})$. A tentative discussion on the “pseudo-inverse” ${\text{SECT}}^{-1}(\mathit{m})$ is provided after Lemma C.4 in Appendix C.

In most shape analysis studies, data are preprocessed by alignment. In Appendix E, we introduce the “ECT alignment” as a preprocessing step before any statistical inference. Throughout the article, we assume that the data have been aligned using this method. The ECT alignment exploits rigid motions, does not rely on landmarks, and is equivalent to the alignment approach outlined in Wang et al. (2021) (Supplementary Section 4). The primary objective of the ECT alignment is to minimize the differences between two shapes that arise from rigid motions. For instance, the molars in Figure 1 were aligned using the ECT alignment. Furthermore, the ECT alignment is compatible with our SECT framework. Appendix E demonstrates that the ECT alignment does not alter the qualitative properties of SECT (e.g., the measurability, Sobolev-regularity, and $\frac{1}{2}$-Hölder continuity).

In applications, it is infeasible to sample infinitely many directions $\nu \in {\mathbb{S}}^{d-1}$ and levels $t\in [0,T]$. For given shapes K, we compute $\text{SECT}(K)(\nu ,t)$ for finitely many directions $\{{\nu }_1,\dots ,{\nu }_{\Gamma }\}\subseteq {\mathbb{S}}^{d-1}$ and levels $\{t_1,\dots ,t_{\Delta }\}\subseteq [0,T]$. To retain information about shapes K, one needs to properly set the numbers of directions and levels (i.e., Γ and Δ). From a theoretical viewpoint, Curry, Mukherjee, and Turner (2022) comprehensively discussed the number Γ of directions needed to recover shapes K from $\text{ECT}(K)$ when K are “piecewise linearly embedded shapes with plausible geometric bounds.” From the numerical perspective, we note the following: (i) Wang et al. (2021) provided detailed simulation studies on the choices of Γ and Δ in sub-image analysis, and a general guideline for setting Γ and Δ in practice was presented in Supplementary Table 1 therein; and (ii) in our Appendix K, we provide detailed numerical experiments on the tradeoffs between the choices of Γ and Δ, the statistical power of our proposed algorithms (Algorithms 1 and 2), and computational cost.

5 Testing Hypotheses on Shapes

In this section, we apply the probabilistic formulation from Section 4 and Lemma C.4 to test hypotheses on shapes. Suppose ${\mathbb{P}}^{(1)}$ and ${\mathbb{P}}^{(2)}$ are two distributions on the measurable space $({\mathcal{S}}_{R,d}^M,\mathcal{F})$. Let ${\mathbb{P}}^{(1)}\otimes {\mathbb{P}}^{(2)}$ be the product probability measure defined on the product σ-algebra $\mathcal{F}\otimes \mathcal{F}$, satisfying ${\mathbb{P}}^{(1)}\otimes {\mathbb{P}}^{(2)}(A\times B)={\mathbb{P}}^{(1)}(A)·{\mathbb{P}}^{(2)}(B)$ for all $A,B\in \mathcal{F}$. To address the motivating question from Section 1.1, we test the following hypotheses (5.1) $$H_0^{*}:{\mathbb{P}}^{(1)}={\mathbb{P}}^{(2)},\mathit{\text{vs}}.H_1^{*}:{\mathbb{P}}^{(1)}\ne {\mathbb{P}}^{(2)},$$(5.1) for example, suppose ${\mathbb{P}}^{(1)}$ and ${\mathbb{P}}^{(2)}$ model the distributions of molars from two genera of primates (Figure 1). Rejecting the $H_0^{*}$ in (5.1) helps distinguish the two genera of primates.

Define $m_{\nu }^{(j)}(t)={\int }_{{\mathcal{S}}_{R,d}^M}\text{SECT}(K)(\nu ,t){\mathbb{P}}^{(j)}(\mathit{\text{dK}})$ for $j\in \{1,2\}$ as the mean functions corresponding to ${\mathbb{P}}^{(1)}$ and ${\mathbb{P}}^{(2)}$. To reject the null $H_0^{*}$ in (5.1) (equivalently, distinguish two collections of shapes), it suffices to reject the null hypothesis H₀ in the following (5.2) $$\begin{array}{cc}\hfill & H_0:m_{\nu }^{(1)}(t)=m_{\nu }^{(2)}(t)\text{ for all }(\nu ,\mathrm{t}),\text{ versus}\hfill \\ \hfill & H_1:m_{\nu }^{(1)}(t)\ne m_{\nu }^{(2)}(t)\text{ for some }(\nu ,\mathrm{t}).\hfill \end{array}$$(5.2)

Analysis of Variance for Functional Data (fdANOVA)

Considering the hypotheses in (5.2) for all directions $\nu \in {\mathbb{S}}^{d-1}$ results in simultaneous multiple-comparisons and inflation of the Type I error. To address this issue, we focus on a specific direction, motivated by the observation that the null hypothesis H₀ in (5.2) is equivalent to $\underset{\nu \in {\mathbb{S}}^{d-1}}{\mathrm{\text{sup}}}\left\{\right||m_{\nu }^{(1)}-m_{\nu }^{(2)}|{|}_{\mathcal{B}}\}=0$. Hence, the direction of interest is defined as (5.3) $${\nu }^{*}\stackrel{\text{def}}{=}\mathrm{\text{arg}}\mathrm{}{\mathrm{\text{max}}}_{\nu \in {\mathbb{S}}^{d-1}}\left\{\left|\right|m_{\nu }^{(1)}-m_{\nu }^{(2)}|{|}_{\mathcal{B}}\right\}.$$(5.3)

Lemma C.4 and (2.2) imply ${\left\{m_{\nu }^{(j)}\right\}}_{j=1}^2\subseteq \mathcal{B}$ for all ν. Lemma C.4, together with (2.1), confirms the existence of a maximizer in (5.3). The maximizer in (5.3) may not be unique. If there are multiple maximizers, we arbitrarily choose one, as this choice does not influence our framework. The null hypothesis $H_0$ in (5.2) is then equivalent to $\left|\right|m_{{\nu }^{*}}^{(1)}-m_{{\nu }^{*}}^{(2)}|{|}_{\mathcal{B}}=0$, where the ${\nu }^{*}$ defined in (5.3) is called a distinguishing direction. Hereafter, we investigate the distribution of $\text{SECT}({\nu }^{*})$.

Based on the discussion above, testing the hypotheses in (5.2) is equivalent to testing $m_{{\nu }^{*}}^{(1)}(t)=m_{{\nu }^{*}}^{(2)}(t)$ for $t\in [0,T]$ using $\text{SECT}({\nu }^{*})$, which is a fdANOVA problem that has been well-studied in the literature (e.g., Zhang 2013, sec. 5.2). However, many state-of-the-art fdANOVA approaches are incompatible with $\text{SECT}({\nu }^{*})$. For example, the Gaussianity of $\text{SECT}({\nu }^{*})$ is not guaranteed (Remark C.1), and the “two-sample problem assumptions” in Section 5.2 of Zhang (2013) may not be satisfied. Besides, the L²-norm-based test (Zhang and Chen 2007) and F-type test (Shen and Faraway 2004) are not preferred when the functional data are not Gaussian (Zhang 2013, chap. 5). Additionally, many fdANOVA methods are time-consuming. For example, tests based on random projections (TRP, Cuesta-Albertos and Febrero-Bande 2010) require the computation of (at least 30) L²-projections for each observed function, followed by the application of appropriate ANOVA tests to these projections. To address the Gaussianity issue and achieve computational efficiency, we propose a method for fdANOVA using the KL expansion. Our test has a foundation that aligns with the probabilistic framework of $\text{SECT}({\nu }^{*})$ in Section 4; it is comparable with the existing methods in terms of size and power (see Appendix J); and it is also computationally efficient, allowing for the permutation test used with our method.

Karhunen–Loève Expansion

Let ${\Xi }_{{\nu }^{*}}^{(j)}(s,t)$ be the covariance function of the stochastic process $\text{SECT}({\nu }^{*})$ corresponding to ${\mathbb{P}}^{(j)}$, for $j\in \{1,2\}$ (see (4.3)). Hereafter, we assume the following, which is true under the null hypothesis $H_0^{*}:{\mathbb{P}}^{(1)}={\mathbb{P}}^{(2)}$ in (5.1).

Assumption 3

(Homoscedasticity). ${\Xi }_{{\nu }^{*}}\stackrel{\text{def}}{=}{\Xi }_{{\nu }^{*}}^{(1)}={\Xi }_{{\nu }^{*}}^{(2)}$, where ${\nu }^{*}$ is defined in (5.3).

This is a standard assumption in the fdANOVA literature (e.g., Zhang 2013, sec. 5.2) and can be tested using the methods proposed by Jia Guo and Zhang (2019).

We define an integral operator on $L^2({[0,T]}^2)$ as $f\mapsto {\int }_0^{T}f(s)·{\Xi }_{{\nu }^{*}}(s,·)\mathit{\text{ds}}$. This operator is compact and self-adjoint (Hsing and Eubank 2015, Theorems 4.6.2 and Example 3.3.4). Moreover, the Hilbert-Schmidt theorem (Reed and Simon 1972, Theorem VI.16) suggests that there is a complete orthonormal basis ${\left\{{\varphi }_l\right\}}_{l=1}^{\infty }$ for $L^2([0,T])$ so that (i) each ${\varphi }_l$ is an eigenfunction with eigenvalue λ_l , (ii) ${\lambda }_1\ge {\lambda }_2\ge \cdots \ge 0$, and (iii) $\underset{l\to \infty }{\mathrm{\text{lim}}}{\lambda }_l=0$. Lemma C.4 and Theorem 7.3.5 of Hsing and Eubank (2015) imply the following KL expansion:

Theorem 5.1

(Karhunen–Loève expansion). (i) For each fixed $j\in \{1,2\}$, we have (5.4) $$\begin{array}{c}\underset{L\to \infty }{\mathrm{\text{lim}}}\underset{t\in [0,T]}{\mathrm{\text{sup}}}{\mathbb{E}}^{(j)}{\left[ \text{SECT}({\nu }^{*},t)-m_{{\nu }^{*}}^{(j)}(t)-\sum _{l=1}^L\sqrt{{\lambda }_l}·Z_l^{(j)}·{\varphi }_l(t)\right]}^2=0,\hfill \end{array}$$(5.4) where $Z_l^{(j)}(K)\stackrel{\text{def}}{=}\frac{1}{\sqrt{{\lambda }_l}}{\int }_0^T\{\text{SECT}(K)({\nu }^{*},t)-m_{{\nu }^{*}}^{(j)}(t)\}·{\varphi }_l(t)\mathit{\text{dt}}$ for $l=1,2,\dots$, and ${\mathbb{E}}^{(j)}$ is the expectation associated with ${\mathbb{P}}^{(j)}$. For each $j\in \{1,2\}$, the random variables ${\left\{Z_l^{(j)}\right\}}_{l=1}^{\infty }$ are defined on the probability space $({\mathcal{S}}_{R,d}^M,\mathcal{F},{\mathbb{P}}^{(j)})$, are mutually uncorrelated, and have mean 0 and variance 1. (ii) There exists $\mathcal{N}\in \mathcal{F}\otimes \mathcal{F}$ so that ${\mathbb{P}}^{(1)}\otimes {\mathbb{P}}^{(2)}(\mathcal{N})=0$ and (5.5) $$\begin{array}{cc}\hfill & {\delta }_l\left(K^{(1)},K^{(2)}\right)\hfill \\ \hfill & \stackrel{\text{def}}{=}\frac{1}{\sqrt{2{\lambda }_l}}{\int }_0^T\left\{\text{SECT}(K^{(1)})({\nu }^{*},t)-\text{SECT}(K^{(2)})({\nu }^{*},t)\right\}·{\varphi }_l(t)\mathit{\text{dt}}\hfill \\ \hfill & ={\theta }_l+\left(\frac{Z_l^{(1)}(K^{(1)})-Z_l^{(2)}(K^{(2)})}{\sqrt{2}}\right),\hfill \\ \hfill & \text{where }{\theta }_l\stackrel{\text{def}}{=}\frac{1}{\sqrt{2{\lambda }_l}}{\int }_0^T\left\{m_{{\nu }^{*}}^{(1)}(t)-m_{{\nu }^{*}}^{(2)}(t)\right\}·{\varphi }_l(t)\mathit{\text{dt}},\hfill \end{array}$$(5.5) for any $(K^{(1)},K^{(2)})\notin \mathcal{N}$ and each fixed $l=1,2,\dots$. The null set $\mathcal{N}$ is allowed to be empty.

Using the KL expansion in (5.4), the random sampling of shapes may be considered, which is discussed in Appendix M.1 and left for future research.

Our Approach

Consider two independent collections of random shapes ${\left\{K_i^{(j)}\right\}}_{i=1}^n\stackrel{\text{iid}}{\sim }{\mathbb{P}}^{(j)}$, for $j\in \{1,2\}$ (i.e., ${\{(K_i^{(1)},K_i^{(2)})\}}_{i=1}^n\stackrel{\text{iid}}{\sim }{\mathbb{P}}^{(1)}\otimes {\mathbb{P}}^{(2)}$). The pairing in $(K_i^{(1)},K_i^{(2)})$ is arbitrary for the following reasons: (i) pairs $(K_i^{(1)},K_i^{(2)})$ and $(K_i^{(1)},K_{i\prime }^{(2)})$ with $i\ne i\prime$ have the same distribution ${\mathbb{P}}^{(1)}\otimes {\mathbb{P}}^{(2)}$, and (ii) numerical experiments in Sections 6 and 7 demonstrate that the performance of our proposed algorithms is numerically invariant to shuffling the index i within each collection ${\left\{K_i^{(j)}\right\}}_{i=1}^n$. Without loss of generality, we assume that all the shapes have been aligned using the “ECT alignment” (Appendix E). Here, we present the theoretical foundation for employing ${\{(K_i^{(1)},K_i^{(2)})\}}_{i=1}^n$ to test the hypotheses in (5.2). This foundation helps address the motivating question from Section 1.1.

Without loss of generality, we assume $(K_i^{(1)},K_i^{(2)})\notin \mathcal{N}$, for all $i=1,2,\dots ,n$, where $\mathcal{N}$ is the null set in Theorem 5.1 satisfying ${\mathbb{P}}^{(1)}\otimes {\mathbb{P}}^{(2)}(\mathcal{N})=0$. Then, we have (5.6) $$\begin{array}{c}{\xi }_{l,i}\stackrel{\text{def}}{=}{\delta }_l\left(K_i^{(1)},K_i^{(2)}\right)={\theta }_l+\left(\frac{Z_l^{(1)}(K_i^{(1)})-Z_l^{(2)}(K_i^{(2)})}{\sqrt{2}}\right),\hfill \end{array}$$(5.6) where δ_l and θ_l are defined in (5.5). Theorem 5.1 implies that, for each fixed l, the random variables ${\left\{{\xi }_{l,i}\right\}}_{i=1}^n$ are iid across $i=1,\dots ,n$ with mean θ_l and variance 1; for each fixed i, the random variables ${\left\{{\xi }_{l,i}\right\}}_{l=1}^{\infty }$ are mutually uncorrelated across $l=1,2,3,\dots$. The following lemma represents the null H₀ in (5.2) using the means ${\left\{{\theta }_l\right\}}_{l=1}^{\infty }$.

Lemma 5.1.

The null H₀ in (5.2) is equivalent to ${\theta }_l=0$ for all positive integers l.

Recall that $\underset{l\to \infty }{\mathrm{\text{lim}}}{\lambda }_l=0$. When eigenvalues λ_l in the denominator of (5.5) are close to zero for large l, the estimated θ_l becomes unstable. Specifically, even if $m_{{\nu }^{*}}^{(1)}(t)\approx m_{{\nu }^{*}}^{(2)}(t)$, an extremely small λ_l can move the corresponding estimated θ_l far away from zero. Using the standard approach in principal component analysis, we focus on ${\left\{{\theta }_l\right\}}_{l=1}^L$ with (5.7) $$L\stackrel{\text{def}}{=}\mathrm{\text{max}}\{1,\tilde{L}\},\text{ where }\tilde{L}\stackrel{\text{def}}{=}\mathrm{\text{min}}\left\{l\in \mathbb{N}|\frac{\sum _{l\prime =1}^l{\lambda }_{l\prime }}{\sum _{l^{″}=1}^{\infty }{\lambda }_{l^{″}}}>0.95\right\}.$$(5.7)

Hence, to test the hypotheses in (5.2) via Lemma 5.1, we test the following (5.8) $$\begin{array}{cc}\hfill & {\widehat{H}}_0:{\theta }_1=\cdots ={\theta }_L=0,\text{ versus }{\widehat{H}}_1:\text{ there exists}\hfill \\ \hfill & \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}l\prime \in \{1,\dots ,L\}\text{such that }{\theta }_{\mathrm{l}\prime }\ne 0.\hfill \end{array}$$(5.8)

Under the null ${\widehat{H}}_0$ in (5.8), for each $l\in \{1,\dots ,L\}$, the central limit theorem indicates that $\frac{1}{\sqrt{n}}{\sum }_{i=1}^n{\xi }_{l,i}$ is asymptotically N(0, 1) when n is large. The mutual uncorrelation in Theorem 5.1 and the asymptotic normality of $\frac{1}{\sqrt{n}}{\sum }_{i=1}^n{\xi }_{l,i}$ provide the asymptotic independence of ${\left\{\frac{1}{\sqrt{n}}{\sum }_{i=1}^n{\xi }_{l,i}\right\}}_{l=1}^L$ across $l=1,\dots ,L$. Then, ${\sum }_{l=1}^L{(\frac{1}{\sqrt{n}}{\sum }_{i=1}^n{\xi }_{l,i})}^2$ is asymptotically ${\chi }_L^2$ under the ${\widehat{H}}_0$ in (5.8). At the asymptotic significance $\alpha \in (0,1)$, we reject the ${\widehat{H}}_0$ if (5.9) $$\begin{array}{c}\sum _{l=1}^L{\left(\frac{1}{\sqrt{n}}\sum _{i=1}^n{\xi }_{l,i}\right)}^2>{\chi }_{L,1-\alpha }^2=\text{ the }1-\alpha \text{ lower quantile}\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{of the }{\chi }_{\mathrm{L}}^2\text{ distribution}.\end{array}$$(5.9)

In applications, neither the mean $m_{\nu }^{(j)}(t)$ nor the covariance ${\Xi }_{\nu }(s,t)$ is known. Hence, the KL expansions in (5.5) cannot be directly used and must be estimated. In Appendix F, we propose a numerical foundation for conducting the asymptotic ${\chi }^2$-test in (5.9) and encapsulate the numerical procedures for the test in Algorithm 1. In all our analyses in Sections 6 and 7, the numerical estimates $\widehat{L}$ (see (F.4) in Appendix F) of the L in (5.7) are smaller than 10. When the $\widehat{L}$ values are large (e.g., several hundred), one may also consider applying the adaptive Neyman test proposed by Fan (1996).

In the simulation studies presented in Tables 1 and J.1, our Algorithm 1 has comparable performance with more than ten existing state-of-the-art fdANOVA methods. Nonetheless, both Algorithm 1 and the existing methods exhibit Type I error inflation (e.g., the rejection rate of Algorithm 1 is 0.118, whereas the significance is 0.05). To mitigate this inflation, we may consider applying the permutation test using one of these methods that is computationally efficient. For example, Górecki and Smaga (2015) proposed a permutation test based on an F-type statistic (FP). Specifically, Górecki and Smaga (2015) approximated each observed function by basis functions via information criteria, and the F-type statistic was approximated by a form conducive to efficiently computing permutation-based p-values. However, the FP also exhibits Type I error inflation (see Tables 1 and J.1). Motivated by the FP, we apply the permutation test to the ${\chi }^2$-statistic defined in (5.9) in the following way: we first apply Algorithm 1 to our original shapes $K_i^{(j)}$ and then repeatedly reapply Algorithm 1 to the shapes with shuffled group labels j. The ${\chi }^2$-test statistic derived from the original data is then compared to that from the shuffled data. A detailed description of our permutation-based approach is presented in Algorithm 2 in Appendix F. Simulations in Section 6 demonstrate that our permutation-based approach eliminates the Type I error inflation encountered by Algorithm 1. The permutation nature of Algorithm 2 is also advantageous for small sample sizes. Note, however, that the power of Algorithm 2 under the alternative is moderately weaker than that of Algorithm 1. Lastly, the runtimes of Algorithms 1 and 2, when applied to simulations, are studied in Appendix K. We present the runtimes when applying the algorithms to real data in Table 2.

Table 1 Rejection rates (from 1000 experiments) for different indices ε (significance $\alpha =0.05$).

Indices ε	0.00	0.01	0.02	0.03	0.04	0.05	0.06	0.08	0.10
Algorithm 1	0.118	0.161	0.315	0.519	0.785	0.910	0.975	0.990	1.000
Algorithm 2	0.046	0.054	0.162	0.343	0.612	0.789	0.931	0.994	1.000
Algorithm 3	0.050	0.050	0.111	0.185	0.335	0.535	0.739	0.983	0.999
FP	0.136	0.153	0.308	0.539	0.810	0.924	0.986	0.997	1.000
TRP-WTPS	0.075	0.091	0.261	0.515	0.790	0.929	0.980	0.997	1.000

NOTE: Appendix J provides a comparison of Algorithms 1, 2, and 3 to other existing fdANOVA methods.

Table 2 P-values of Algorithms 1, 2, and 4 for the dataset of mandibular molars.

	Algorithm 1	Algorithm 2	Algorithm 4
Tarsius vs. Microcebus	$<{10}^{-3}$	$<{10}^{-3}$	$<{10}^{-3}$
Tarsius vs. Mirza	$<{10}^{-3}$	$<{10}^{-3}$	0.001
Tarsius vs. Saimiri	$<{10}^{-3}$	$<{10}^{-3}$	$<{10}^{-3}$
Microcebus vs. Mirza	$<{10}^{-3}$	0.009	0.004
Microcebus vs. Saimiri	$<{10}^{-3}$	$<{10}^{-3}$	$<{10}^{-3}$
Mirza vs. Saimiri	$<{10}^{-3}$	$<{10}^{-3}$	$<{10}^{-3}$
Tarsius vs. Tarsius	0.196 (0.220)	0.496 (0.294)	0.527 (0.273)
Overall runtimes (in hours)	$\approx 3$	$\approx 3$	$\approx 20$

NOTE: In the last row, we present the overall runtime for conducting all hypothesis testing tasks using each algorithm.

6 Experiments Using Simulations

We present simulations showing the performance of our Algorithms 1 and 2. In addition, we compare our algorithms with the “randomization-style null hypothesis significance test (NHST)” (Robinson and Turner 2017), the TRP using Wald-type permutation statistic (TRP-WTPS, Cuesta-Albertos and Febrero-Bande 2010; Pauly, Brunner, and Konietschke 2014), and the FP. Details of the randomization-style NHST are given in Appendix G and referred to as Algorithm 3. The application of the FP and TRP to the SECT is described in Section 5. We implement the FP and TRP-WTPS using the R package fdANOVA with its default parameters as recommended by Górecki and Smaga (2019). Additional simulations comparing our proposed algorithms and other existing fdANOVA methods are presented in Appendix J.

We focus on a family of distributions ${\left\{{\mathbb{P}}^{(\epsilon )}\right\}}_{0\le \epsilon \le 0.1}$ with shapes ${\left\{K_i^{(\epsilon )}\right\}}_{i=1}^n\stackrel{\text{iid}}{\sim }{\mathbb{P}}^{(\epsilon )}$ via (6.1) $$\begin{array}{c}{K}_i^{(\epsilon )}\stackrel{\text{def}}{=}\left\{x\in {\mathbb{R}}^2\left|\underset{y\in S_i^{(\epsilon )}}{\mathrm{\text{inf}}}\right||x-y||\le 0.2\right\},\text{ where}\hfill \end{array}$$(6.1)

$S_i^{(\epsilon )}\stackrel{\text{def}}{=}\left\{\left(\frac{2}{5}+a_{1,i}·\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}t,b_{1,i}·\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}t\right)|\frac{1-\epsilon }{5}\pi \le t\le \frac{9+\epsilon }{5}\pi \right\}\cup \left\{\left(-\frac{2}{5}+a_{2,i}·\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}t,b_{2,i}·\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}t\right)|\frac{6\pi }{5}\le t\le \frac{14\pi }{5}\right\}$ and ${\{a_{1,i},a_{2,i},b_{1,i},b_{2,i}\}}_{i=1}^n\stackrel{\text{iid}}{\sim }N(1,{0.05}^2)$. The ε denotes the dissimilarity between ${\mathbb{P}}^{(\epsilon )}$ and ${\mathbb{P}}^{(0)}$. For each $\epsilon \in [0,0.1]$, through the discussion in Section 5, we test the following hypotheses via fdANOVA methods (i.e., FP, TRP-WTPS, Algorithms 1, and 2) $$\begin{array}{cc}\hfill & H_0:m_{\nu }^{(0)}(t)=m_{\nu }^{(\epsilon )}(t)\text{ for all}(\nu ,\mathrm{t})\in {\mathbb{S}}^{\text{d‐1}}\times [0,\mathrm{T}]\text{ versus}\hfill \\ \hfill & H_1:m_{\nu }^{(0)}(t)\ne m_{\nu }^{(\epsilon )}(t)\text{ for some }(\nu ,\mathrm{t}),\hfill \end{array}$$ where the mean $m_{\nu }^{(\epsilon )}(t)\stackrel{\text{def}}{=}{\int }_{{\mathcal{S}}_{R,d}^M}\text{SECT}(K)(\nu ,t){\mathbb{P}}^{(\epsilon )}(\mathit{\text{dK}})$, and the null hypothesis H₀ is true when $\epsilon =0$. We also test $H_0^{*}:{\mathbb{P}}^{(0)}={\mathbb{P}}^{(ϵ)}\mathit{\text{vs}}.{\mathbb{P}}^{(0)}\ne {\mathbb{P}}^{(ϵ)}$ using Algorithm 3.

We set T = 3, directions ${\nu }_p={(\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}\frac{p-1}{4}\pi ,\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}\frac{p-1}{4}\pi )}^{⊺}$ for $p\in \{1,2,3,4\}$, levels $t_q=\frac{T}{50}q$ for $q\in \{1,\dots ,50\}$ (i.e., Γ = 4 and Δ = 50 in Algorithms 1, 2, and 3), the confidence level 95% (i.e., $\alpha =0.05$), and the number of permutations Π = 1000. For each $\epsilon \in$ {0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08, 0.1}, we independently generate two collections ${\left\{K_i^{(0)}\right\}}_{i=1}^n\stackrel{\text{iid}}{\sim }{\mathbb{P}}^{(0)}$ and ${\left\{K_i^{(\epsilon )}\right\}}_{i=1}^n\stackrel{\text{iid}}{\sim }{\mathbb{P}}^{(\epsilon )}$ through (6.1) with the number of shape pairs set to n = 100, and we compute the SECT of each generated shape in directions ${\left\{{\nu }_p\right\}}_{p=1}^4$ and at levels ${\left\{t_q\right\}}_{q=1}^{50}$. We then implement the fdANOVA methods and Algorithm 3 to these computed SECT statistics and get the corresponding Accept/Reject outputs. We repeat this procedure 1000 times and report the rejection rates across all 1000 replicates for each ε in Table 1. The rejection rates are also visually presented in Figure 3. We choose Γ = 4 as the number of directions in our simulations based on the following observation: in Appendix K, we experiment with all combinations of $\Gamma \in \{2,4,8\},\Delta \in \{25,50,100\}$, and $n\in \{25,50,100\}$. When Δ = 50 and n = 100, the number Γ = 4 is sufficiently large for our Algorithms 1 and 2 to distinguish ${\mathbb{P}}^{(0)}$ from ${\mathbb{P}}^{(\epsilon )}$ with $\epsilon >0$ using the significance level $\alpha =0.05$. Moreover, this choice allows us to demonstrate that even a relatively small number of directions (e.g., Γ = 4) is sufficient for implementing our Algorithms 1 and 2.

Fig. 3 (Left panel) The relationship between ε and the rejection rates computed via Algorithms 1, 2, 3 (see Table 1), and 12 existing fdANOVA methods (see Table J.1 in Appendix J for details on the existing fdANOVA methods). The (red) dashed line presents the significance level $\alpha =0.05$. (Right panel) The shapes in the first row are from ${\mathbb{P}}^{(0)}$, and the shapes in the second row are from ${\mathbb{P}}^{(0.08)}$.

The results in Table 1 and Figure 3 demonstrate that our proposed algorithms are effective at detecting the difference between ${\mathbb{P}}^{(\epsilon )}$ and ${\mathbb{P}}^{(0)}$ in terms of distinguishing their mean functions. Notably, our algorithms (especially Algorithm 2) tend to avoid falsely detecting differences between shape-generating distributions under the null hypothesis (i.e., $\epsilon =0$). As ε increases, ${\mathbb{P}}^{(\epsilon )}$ deviates from ${\mathbb{P}}^{(0)}$, and the power of our algorithms in detecting the deviation increases. When $\epsilon \ge 0.08$, the power of Algorithms 1 and 2 exceeds 0.99. For all the ε, it is difficult to see the deviation of ${\mathbb{P}}^{(\epsilon )}$ from ${\mathbb{P}}^{(0)}$ visually. For instance, by merely observing the shapes in Figure 3, one might find it hard to differentiate between the shape collections generated by ${\mathbb{P}}^{(0)}$ (blue) and ${\mathbb{P}}^{(0.08)}$ (pink). However, in more than 99% of the simulations, our algorithms detect the difference between the two distributions. We also randomly shuffle the index i within each collection ${\left\{K_i^{(\epsilon )}\right\}}_{i=1}^n$ and apply Algorithms 1 and 2 to the shuffled collections. The results obtained from the unshuffled and shuffled shape collections, respectively, are nearly identical. Algorithm 3 performs well in detecting the discrepancy between ${\mathbb{P}}^{(0)}$ and ${\mathbb{P}}^{(\epsilon )}$. However, its power under the alternative hypotheses (i.e., $\epsilon >0$) is weaker than that of our Algorithms 1 and 2. Moreover, Algorithms 1 and 2 exhibit performance comparable to twelve existing state-of-the-art fdANOVA methods (see Table 1, Figure 3, and Table J.1 in Appendix J).

7 Applications

We first apply our proposed Algorithms 1 and 2 to the MPEG-7 shape silhouette database (Sikora 2001) as a toy example. Details of this are provided in Appendix I. This analysis shows that our proposed algorithms can distinguish between shape classes in the silhouette database and do not falsely identify signals when there are no differences between groups.

In this section, we apply our algorithms to address the motivating question in Section 1.1. Specifically, we use Algorithms 1 and 2 to distinguish between the four categories of mandibular molars in Figure 1 that are from four genera of primates. The shapes in Figure 1 come from two suborders of primates: Haplorhini and Strepsirrhini (see Figure 1). In the haplorhine suborder collection, 29 molars came from the genus Tarsius (yellow panels in Figure 1), and 9 molars came from the genus Saimiri (grey panels in Figure 1). In the strepsirrhine collection, 11 molars came from the genus Microcebus (blue panels in Figure 1), and 6 molars came from the genus Mirza (green panels in Figure 1).

Before applying Algorithms 1 and 2, we preprocess the raw triangle mesh data of the surfaces of the molars by aligning them through the ECT alignment approach detailed in Appendix E. The aligned molars are presented in Figure 1. We apply our Algorithms 1 and 2 to the preprocessed molars. For each aligned molar, we compute its SECT for 2918 directions; in each direction, we use 200 sublevel sets. To compare any pair of molar groups, as a proof of concept, we select the smaller size of the two groups as the sample size input n in our algorithms. For example, when comparing the Tarsius and Microcebus groups, we choose n = 11; that is, we compare the first 11 molars of the Tarsius group to all the molars in the Microcebus group. We apply our algorithms to the four groups of molars and present the results in Table 2. The p-values in Table 2 are either ${\chi }^2$-test p-values (Algorithm 1) or permutation-test p-values (Algorithm 2 with 1000 permutations). The small p-values (P < 0.05) in Table 2 show that our proposed algorithms can distinguish the four different genera of primates. Since the genera Microcebus and Mirza belong to the same suborder Strepsirrhini (see Figure 1), the p-value from Algorithm 2 is comparatively large when comparing molars from these two groups. In comparison, although the Tarsius and Saimiri both belong to the suborder Haplorhini, the molars of the two genera look different. Specifically, the paraconids (i.e., the cusp highlighted in red in Figure 1) are only retained by the genus Tarsius and, thus, are a key reason for the small p-values ($P<{10}^{-3}$) when comparing with molars from the Saimiri. Other small p-values ($P<{10}^{-3}$) in our analyses are a result of the corresponding genera belonging to different suborders.

In addition to testing the difference between genera, we apply our algorithms within the genus Tarsius. Specifically, we focus on the first 28 molars in the Tarsius group. We randomly split the 28 molars into two halves and apply Algorithms 1 and 2 to test the difference between the two halves. We repeat the random splitting procedure 100 times and present the corresponding p-values in Table 2. The results are summarized by their mean and standard deviation (in parentheses). These p-values show that our proposed Algorithm 2 tends to avoid the Type I error for the molars from the genus Tarsius.

Landmark methods are widely used in geometric morphometrics. One state-of-the-art approach is the “Gaussian process landmarking (GPL)” algorithm (Gao et al. 2019; Gao, Kovalsky, and Daubechies 2019) which can automatically sample landmarks on the surfaces of the molars in Figure 1. Gao et al. (2019) showed that these sampled landmarks could induce a continuous Procrustes distance to measure the dissimilarity between molars. A permutation test can be derived using the Procrustes distance induced by the GPL algorithm. This test is detailed in Appendix H and is encapsulated by Algorithm 4. We use the GPL-based Algorithm 4 to differentiate the four collections of molars. For this, we use the MATLAB code from the GitHub repository provided by Gao et al. (2019) to compute the Procrustes distance. Performance of Algorithm 4 is in Table 2, which shows that the GPL-based method and our Algorithm 2 have comparable performance. However, repeatedly computing the Procrustes distance is time-consuming. Hence, Algorithm 2 is more computationally efficient than Algorithm 4 while achieving similar performance (see the last row of Table 2).

We want to note that, in addition to the GPL algorithm, many other existing methods can be applied to measure dissimilarity between molars, including parameterized surfaces (Kurtek et al. 2010, 2011) and the approaches from computational anatomy (Grenander and Miller 1998). Similarly, the parameterized curves (Kurtek et al. 2012) can also be used to analyze the silhouette database in Appendix I. An even more comprehensive comparison of our algorithms with the entire edifice of existing methods is left for future research.

8 Conclusions and Discussions

In this article, we established the mathematical foundations for the randomness of shapes via the SECT. Specifically, (i) $({\mathcal{S}}_{R,d}^M,\mathcal{B}(\rho ),\mathbb{P})$ was constructed as the underlying probability space; (ii) the SECT was modeled as a $C({\mathbb{S}}^{d-1};\mathcal{H})$-valued random variable. We further demonstrated several properties of the SECT ensuring its KL expansion, which led to a ${\chi }^2$-statistic for testing hypotheses on random shapes. We bridged the fdANOVA and TDA. Simulation studies corroborated our mathematical derivations and showed the performance of our hypothesis testing algorithms. Our approach was shown to be powerful in detecting the difference between two shape-generating distributions. We applied our proposed algorithms to silhouette and primate molar datasets. Importantly, our simulations when $\epsilon =0$, together with the applications to the molars and the silhouette database, indicate that our algorithms tend to avoid falsely detecting differences between shape-generating distributions when there are none. Using the molars in Figure 1, we compared the performance of our algorithms to a permutation test based on a state-of-the-art landmarking algorithm (Gao et al. 2019; Gao, Kovalsky, and Daubechies 2019), underscoring the efficiency of our algorithms. We enumerate potential future research areas in Appendix M, for example, the fdANOVA methods can be used for brain connectivity (Chen et al. 2024; Meng and Eloyan 2024) via topological summaries.

Acknowledgments

We are grateful to the Editor, Associate Editor, and three Referees for their thorough review of our article and the insightful suggestions that have tremendously improved its quality. We want to thank Dr. Matthew T. Harrison from the Division of Applied Mathematics at Brown University for useful comments and suggestions. KM wants to thank Mattie Ji from the Department of Mathematics at Brown University for her insightful comments. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Supplementary Materials

The supplementary materials provide the proofs of theorems, further data analysis, and future research topics.

Data Availability Statement

The source code for implementing the simulation studies and applications is publicly available online at https://github.com/JinyuWang123/TDA.git.

Disclosure Statement

The authors report there are no competing interests to declare.

Additional information

Funding

LC would like to acknowledge the support of a David & Lucile Packard Fellowship for Science and Engineering. Research reported in this publication was partially supported by the National Institute On Aging of the National Institutes of Health under Award Number R01AG075511.