Landmark-free, parametric hypothesis tests regarding two-dimensional contour shapes using coherent point drift registration and statistical parametric mapping

Landmarks univariate (UV) analysis

Landmarks were defined for each dataset as depicted in Fig. 2. Both the number of landmarks (Table 1) and their locations were selected in an ad hoc manner, with the qualitative requirement of readily identifiable, homologous locations. The ultimately selected landmarks arguably span a representative range of landmarking possibilities.

One operator used a mouse to manually digitize the landmarks for each of the 90 shapes (10 shapes for each of 9 datasets). The operator was ignorant of the final shape groupings for the ultimate two-sample tests (see below), implying that the landmarking was performed without grouping bias.

The landmarks were spatially registered using Generalized Procrustes Analysis (GPA) (Gower, 1975), and the resulting registered landmarks were analyzed in a univariate manner, using Procrustes ANOVA (Goodall, 1991)—a method which considers the variance in the Procrustes distance across a dataset. Note that the Procrustes distance is a scalar quantity that summarizes shape difference, and thus that this method is univariate. GPA and Procrustes ANOVA were both conducted using the geomorph package for R (Adams & Otárola-Castillo, 2013).

Landmarks mass-multivariate (mass-MV) analysis

This approach was identical to the Landmarks-UV approach described above, except for statistical analysis. The two-sample Hotelling’s T² statistic was calculated for each landmark according to its definition: (1)${\displaystyle T_i^2=\frac{n_1{n}_2}{n_1+n_2}\left(\right.{\overline{r}}_{1i}-{\overline{r}}_{2i}{\left)\right.}^{\top }{W}_i^{-1}\left(\right.{\overline{r}}_{1i}-{\overline{r}}_{2i}\left)\right.}$

where i indexes landmarks, the subscripts “1” and “2” index the two groups, n is sample size, ${\overline{r}}_i$ is the mean position vector of landmark i, and W_i is the pooled covariance matrix for landmark i: (2)${\displaystyle W_i=\frac{1}{n_1+n_2-2}\left(\sum _{j=1}^{n_1}\left(r_{1ij}-{\overline{r}}_{1i}\right){\left(r_{1ij}-{\overline{r}}_{1i}\right)}^{\top }+\sum _{j=1}^{n_2}\left(r_{2ij}-{\overline{r}}_{2i}\right){\left(r_{2ij}-{\overline{r}}_{2i}\right)}^{\top }\right)}$

where the i index is dropped for convenience in Eq. (2).

Statistical inference was conducted in a mass-multivariate manner, using Statistical Parametric Mapping (SPM) (Friston et al., 2007). SPM bases statistical inferences on the distribution of the maximum T² value $\left(T_{max}^2\right)$, which can be roughly interpreted as the largest landmark effect, and which is defined as: (3)${\displaystyle T_{max}^2\equiv \underset{i\in L}{max}{T}_i^2}$

where L is the number of landmarks.

SPM provides a parametric solution to the distribution of $T_{max}^2$ under the null hypothesis, so significance can be assessed by determining where in this distribution the observed $T_{max}^2$ lies. Classical hypothesis testing involves the calculation of a critical threshold (T²)_critical , defined as the (1 − α)th percentile of this distribution, and all landmarks whose T² values exceed (T²)_critical are deemed significant at a Type I error rate of α. This is a correction for multiple comparisons (i.e., across multiple landmarks) that is ‘mass-multivariate’ in the following sense: ‘mass’ refers to a family of tests, in this case a family of landmarks, and ‘multivariate’ refers to a multivariate dependent variable, in this case is a two-component position vector. This is similar to traditional corrections for multiple comparisons like Bonferroni corrections, with one key exception: rather than using the total number of landmarks L as the basis for the multiple comparisons correction, as the Bonferroni correction does, SPM instead solves the mass-MV problem by assessing the correlation amongst neighboring landmarks or semilandmarks, and using the estimated correlation to provide a less severe correction than the Bonferroni correction, unless there is no correlation, in which case the SPM and Bonferroni corrections are equivalent.

Contours univariate (UV) analysis

Similar to the Landmarks UV approach, this approach ultimately conducted Procrustes ANOVA, but did so on contour data rather than landmark data. This was achieved through two main processing steps: coherent point drift (CPD) point set registration (Fig. 4) and optimum roll correspondence (Fig. 5). Coherent point drift (CPD) (Myronenko & Song, 2010) is a point set registration algorithm that spatially aligns to sets of points that belong to the same or a similar object. Neither an equal number of points nor homologous points are required (Fig. 4), making this approach useful for contours that have an arbitrary number of points.

Since contour points from arbitrary datasets may generally be unordered (Fig. 5A), we started our analyses by randomly ordering all contour points, then applying CPD to the unordered points. We acknowledge that many 2D contour datasets consist of ordered points—including those in the database used for this study (Carlier et al., 2016)—but since 3D surface points are much more likely to be unordered, we regard unordered point support as necessary for showing that the proposed method is generalizable to 3D analyses. Following CPD, we re-ordered the points using parametric surface modeling (Bingol & Krishnamurthy, 2019), which fits a curved line to the contour, and parameterizes the contour using position u, where u ranges from zero to one (Fig. 6). This contour parameterization results in a continuous representation of the contour, from which an arbitrary number of ordered points (Fig. 5B) can be used to discretize the contour of each shape for subsequent analysis. We used NURBS parameterization with B-spline interpolation (Bingol & Krishnamurthy, 2019) to calculate specific contour point locations. We then applied an optimum roll transformation, which found the value of u for one contour that minimized the deformation energy across the two contours (Figs. 5C, 5D).

We repeated contour parameterization, ordering, and optimum roll correspondence across all contour shapes, using the shape with the maximum number of contour points in each dataset as the template shape to which the nine other shapes were registered. Note that this registration procedure is unrelated to the traditional landmark analyses described in ‘Landmark UV analysis’ above, for which an equal number of points is a requirement of registration and analysis. The correspondence analysis step resulted in an equal number of contour points, upon which we conducted Procrustes ANOVA.

Contours mass-multivariate (mass-MV) analysis

This approach was identical to the Contours-UV approach, with the exception of statistical analysis, which we conducted using SPM as outlined above. Unlike the landmark data above, which are generally spatially disparate, contour points are spatially proximal, and neighboring points tend to displace in a correlated manner. For example, if one contour point in a specific shape lies above the mean point location, its immediate neighbors also tend to lie above the mean location). SPM leverages this correlation to reduce the severity of the multiple comparisons correction, and SPM solutions converge to a common (T²)_critical regardless of the number of contour points, provided the number of contour points is sufficiently large to embody the spatial frequencies of empirical interest, as outlined in classical signal processing theory (Nyquist, 1928).

As SPM uses parametric inference to calculate the critical T² threshold, and Procrustes ANOVA uses nonparametric inference, we also conduct Contours Mass-MV analysis using statistical non-parametric mapping (Nichols & Holmes, 2002), which uses permutation to numerically build the $T_{max}^2$ distribution under the null hypothesis. This permutation approach converges to the parametric solution when the residuals are normally distributed (i.e., point location variance follows an approximately bivariate Gaussian distribution). All SPM analyses were conducted in spm1d (Pataky, 2012); note that one-dimensional SPM is sufficient because the contour domain (U) is one-dimensional (Fig. 6).