Statistical analysis of locally parameterized shapes

08/18/2021 ∙ by Mohsen Taheri, et al. ∙ University of Stavanger 0

The alignment of shapes has been a crucial step in statistical shape analysis, for example, in calculating mean shape, detecting locational differences between two shape populations, and classification. Procrustes alignment is the most commonly used method and state of the art. In this work, we uncover that alignment might seriously affect the statistical analysis. For example, alignment can induce false shape differences and lead to misleading results and interpretations. We propose a novel hierarchical shape parameterization based on local coordinate systems. The local parameterized shapes are translation and rotation invariant. Thus, the inherent alignment problems from the commonly used global coordinate system for shape representation can be avoided using this parameterization. The new parameterization is also superior for shape deformation and simulation. The method's power is demonstrated on the hypothesis testing of simulated data as well as the left hippocampi of patients with Parkinson's disease and controls.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In statistical shape analysis, besides classifying shapes, detecting and characterizing locational differences between two populations of shapes is a matter of special interest. Particularly in medical applications, shape analysis has the power to shed light on organ deformations, supporting diagnosis and treatment. Therefore, in a considerable number of articles, researchers try to find a set of corresponding

geometric object properties (GOPs) among shapes to explore regional dissimilarities. Corresponding GOPs can be represented in various ways, for example, by a set of landmarks on or inside the objects (dryden2016statistical, Ch.1), a point distribution model (PDM) such as spherical harmonics PDM (SPHARM-PDM) (Styner2006)

, or as a subset of skeletal structures (e.g., a set of internal vectors)

(siddiqi2008medial, ch.8). For instance, (cates2007shape) introduced entropy-based surface sampling PDM to compare brain objects of a group of patients with schizophrenia and a healthy control group (CG). A large variety of studies can be found using PDM models to study human organs, e.g., (alhadidi20123d; oguz2008cortical; achterberg2014hippocampal). For the skeletal structure, (pizer2013nested) introduced skeletal representation (s-rep) and consequently discrete s-rep (ds-rep) where ds-rep can be considered as a penalized version of medial representation (m-rep) (pizer1999segmentation) (see Figure 2). An s-rep reflects the object’s interior and defines a smooth implied boundary. (schulz2016non) proposed a hypothesis test for ds-rep and studied the hippocampal differences between schizophrenia and CG.

Generally, these types of studies for shape analysis share three main steps. First, a preprocessing step, where both groups of objects are aligned based on their corresponding GOPs. The objective of the alignment is to quantify shape differences purely without locational information. Thus, the distance between objects is minimized to make them invariant under the act of Euclidean similarity transformations translation, rotation, and scaling. Usually, alignment takes place by generalized Procrustes analysis (GPA) (dryden2016statistical, Ch.7). Second, a set of partial hypothesis tests on GOPs to verify significant locational differences between the two groups. Third, multiple testing methods are applied to control false positives such as family wise error rate (FWER), e.g., (bonferroni1936teoria) or false discovery rate (FDR), e.g., (benjamini1995controlling).

There are fundamental concerns with the alignment during the preprocessing step. Figure 1 illustrates some of the concerns with a simple example. Figure 0(a) shows two ellipsoidal shapes. The blue shape is an ellipsoid, and the red shape is like a boomerang. By visual inspection, we notice that they are different. However, if we think of it as an open arm (blue) and closed arm (red), where each arm consists of three separate parts, namely the upper arm, elbow, and forearm, we observe only a difference around the elbow. The upper and forearm are remaining unchanged. Our visual cortex tries to understand the two objects as a whole, i.e., in a global coordinate system (GCS) , even though the red shape in Figure 0(a) is only a locally deformed version of the blue shape. In other words, both shapes are identical at the top and bottom parts but different only at the middle. Now, let us assume both shapes are sampled at 24 positions. By adding independent random noise to each sample point, we simulated 20 PDMs as depicted in Figure 0(a). In Figure 0(b), shapes are aligned without scaling by GPA. The distribution of almost all of the corresponding points (e.g., point No. 5) are remarkably separated that leads to a large number of false positives in mean difference hypothesis testing of local distributions. This phenomenon can also be observed in the analysis of brain structures like the hippocampus, where alignment leads to false positives between CG and treatment groups, further discussed in Section 4.3. In Figure 0(c), shapes are aligned based on top and bottom parts with weighted GPA (dryden2016statistical, Sect. 7.6.3). Although, it seems weighted GPA is more reasonable than GPA, defining a suitable covariance structure for weighted GPA is not explicit.

(a) Ellipsoidal PDMs

(b) GPA alignment

(c) Weighted GPA
Figure 1: Problem of false positives due to alignment. (a) Red and blue indicate two populations of PDMs. Small crosses are the mean centroids. (b,c) Separation of corresponding local distributions.

So far, we have exposed some problems caused by alignment and by understanding shapes in a GCS. We propose a novel hierarchical shape parameterization based on local coordinate systems to overcome the problems. The local parameterization has three major contributions. First, because the representation is translation and rotation invariant, the fundamental issues of alignment are avoided. Second, it understands shapes locally. For example, in Figure 0(a), we would only detect differences in the middle where deformation occurred. Third, it naturally supports the interpretation of potential shapes differences, e.g., as bending or twisting. Forth, it facilitates skeletal deformation and simulation.

Basically, local frames can be defined for different types of object representation. However, a representation that is in particular very suitable is ds-rep. Thus, in this manuscript, we focus on introducing local frames with the application for ds-rep. The paper is structured as follows. Basic notations and amenities of ds-reps are summarized in Section 2. Figure 1(b) shows a ds-rep of a hippocampus.

(a) m-rep

(b) ds-rep
Figure 2: Skeletal structure. (a) 2D m-reps of two ellipsoidal objects. and are corresponding spokes with unit directions and . (b) A fitted ds-rep to a left hippocampus’s mesh. Green, cyan, magenta, and yellow respectively indicate skeletal sheet, up spokes, down spokes, and crest spokes.

Section 2.1.1 introduces the conventional definition of s-rep and ds-rep. The conventional definition is based on GCS with the discussed challenges. In Section 2.1.2 and Section 2.1.3, we introduce the proposed LP-ds-rep. The proposed hierarchical local parameterization of ds-rep, called LP-ds-rep, extracts GOPs with more details and supports sensitive hypothesis testing that is not biased by prior alignment. Section 3 explains euclideanization procedure of spherical data by principal nested spheres (PNS) (jung2012analysis), and adapt the non-parametric hypothesis testing method plus controlling false positive from (Styner2006). Section 4.2 discusses skeletal deformation and simulation by LP-ds-rep. In Section 4.3, we study hippocampal differences between a group of patients with early Parkinson’s disease (PD) and CG given both parameterizations. Besides, we compare the results and show the advantages of our method on simulated data. Finally, we summarize and conclude the work in Section 5.

2 Skeletal representation

The m-rep and its properties have been extensively studied in the literature (pizer1999segmentation; fletcher2004principal; siddiqi2008medial) all the way back to Blum’s medial structure (blum1967transformation). Figure 1(a) illustrates 2D m-reps for the previously discussed shapes with smooth boundaries. Briefly, an m-rep is a sample set of medial axis and spokes. The medial axis of object is where is the boundary of , is the minimum Euclidean distance between and , and is the cardinality sign. In other words, the medial axis is centers of all inscribed spheres of , bi-tangent or multi-tangent to . A spoke is a vector connecting the center of an inscribed sphere to where its length is equal to the sphere’s radius. Thus, an m-rep reflects the interior object properties such as local widths and directions. However, as pointed out in (pizer2013nested), the m-rep is sensitive to boundary noise because every protruding boundary kink results in additional medial branches. This sensitivity affects m-rep correspondence among a population as two versions of the same objects can result in significantly different m-reps. Thus, (pizer2013nested) relaxed the mentioned condition and defined s-rep. As described in (liu2021fitting), for a slabbed-shaped object like the hippocampus, an s-rep has the form , where is a smooth ellipsoidal skeletal sheet and is the field of non-crossing vectors (spokes) on . The tail of each spoke is at a point and its tip is at the object boundary . can be assumed as a combination of three components , , and such that is a fold curve divides into two collocated submanifolds and . and map and respectively to two sides of the object’s boundary considered as northern and southern part. Also, maps to the crest part of the boundary. We call a spoke an up spoke, down spoke, or crest spoke if it belongs to , , or respectively.

The relaxed conditions assure stability in the branching structure and thus good case-to-case correspondence across a population of s-reps. In practice, we sample a finite number of corresponding spokes from an s-rep to obtain a ds-rep. The conventional ds-rep parameterization is understood in a GCS explained in more detail in Section 2.1.1. Afterward, a novel parameterization based on a hierarchical structure of the local frames is introduced. Also, we name the conventional parameterization as globally parameterized ds-rep (GP-ds-rep), and the new parameterization as locally parameterized ds-rep (LP-ds-rep). Further, and denote GP-ds-rep and LP-ds-rep respectively.

2.1 Parameterizations

2.1.1 GP-ds-rep

There are different ways to fit and parameterize a ds-rep. A current implementation described in detail by (liu2021fitting)

is available under the open-source toolbox SlicerSALT (

http://salt.slicer.org). A GP-ds-rep is a tuple where , and are th spoke’s tail position, direction, and length respectively where is the unit sphere in arbitrary dimension , , and is the number of spokes. Based on the current model fitting, some spokes share a common tail position, so we have where , and . The set forms an configuration matrix representing the skeletal PDM. Assume as identity matrix and as vector of ones. Location and scale can be removed by centering and normalizing skeletal PDM to obtain pre-shape , where is the centering matrix, and is the Euclidean norm. Since , the pre-shape lives on the hypersphere dryden2016statistical, Ch.2; schulz2016non. Thus a GP-ds-rep lives on a manifold as a direct product of Riemannian symmetric spaces, namely where indicates the pre-shape space of the skeletal PDM, is the space of spokes’ directions, and is the space of spokes’ lengths and the scaling factor. Note that spoke positions and directions in this parameterization are in a GCS independent of the skeletal sheet structure. This lead to the discussed challenges in statistical analysis because alignment is necessary.

For m-rep, a semi-local parameterization was proposed by (fletcher2003statistics) based on local frames , where is normal to the medial sheet at , is the bisector direction of two equal-length spokes with common position, , and is the 3D rotation group. Spokes’ directions are defined relative to the local frames by the angle between and the spokes (see Figure 3(a)).

Because the direction of and depends on the spokes’ directions, if then takes an arbitrarily direction that violates the uniqueness and consistency of the fitted frame. Besides, the spokes’ tail positions and frame directions are in GCS. Thus for statistical shape analysis, pre-alignment is still necessary.

Inspired by both, Cartan’s moving frames on space curves (alma991001580569705596), and Fletcher’s semi-local parametrization, we propose a fully local ds-rep parameterization. By utilizing the inherent hierarchical structure of ds-reps, we provide a consistent definition of local frames that is independent of GCS and also avoids arbitrarily frame rotation. This can be done by introducing a leaf-shaped skeletal structure of the skeletal sheet according to the ellipsoidal design of the ds-rep applied in (liu2021fitting) (see Figure 5(Top)).

2.1.2 Local frames

For GP-ds-rep fitting, (liu2021fitting) first deformed an object to an ellipsoid by mean curvature flow. Then by inverse mean curvature flow, warped the ellipsoid’s GP-ds-rep to the target object. Finally, spokes are refined such that the implied boundary becomes as close as possible to the real boundary. As a result, the implied boundary corresponds to the ellipsoid’s boundary, and the skeletal sheet corresponds to the ellipsoid’s skeletal sheet. In other words, the structure of the fitted GP-ds-rep is associated with the ellipsoid’s GP-ds-rep i.e., ellipsoid’s medial axis. By assuming a good correspondence between the ellipsoid’s medial axis and the object’s skeletal sheet, we design a hierarchical structure for the ellipsoid’s medial axis and expand it to the object’s skeletal sheet. Then on the basis of the obtained structure, we define consistent fitted frames in a population of GP-ds-reps.

In this article, we only consider slabular shapes corresponding to an eccentric ellipsoid with principal radii such that . As discussed above, up and down spokes correspond to the ellipsoid’s northern and southern side, while crest spokes correspond to the ellipsoid’s crest.

As illustrated in Figure 3(left), the medial axis of an ellipsoid is an ellipse (i.e., a 2D ellipsoid). The ellipse is symmetric and has a symmetric 2D m-rep. The m-rep consists of medial points located on a straight medial line and a set of spokes. The middle point of the medial line is the ellipsoid’s centroid (i.e., center of gravity ). The boundary of the ellipse is the skeletal sheet’s fold (i.e., skeletal edge) of the ellipsoid. The extension of m-rep spokes connecting the fold to the ellipsoid boundary is a subset of crest spokes. We call the members of this subset crest-midline spokes. For the model fitting, (liu2021fitting) deformed the ellipsoid to the object. After the deformation, the flat ellipsoid’s medial axis transforms to a nonlinear surface as a 2D manifold (i.e., topologically a 2-dimensional disk). Consequently, straight lines on it (e.g., medial line and m-rep spokes) become curves. Also (liu2021fitting) assumed more or less a diffeomorphic transformation. Thus the generated curves do not cross each other.

We call the deformed medial line the spine, and deformed m-rep spokes veins. Thus, veins are a set of non-crossing curves emanating from the spine. During the deformation, the ellipsoid’s centroid moves and ultimately rests in the middle of the object. We assume the displaced centroid as an intrinsic centroid, and call it skeletal centroid or s-centroid. Figure 3 provides an intuition about the ellipsoid’s medial axis deformation.

Figure 3: Ellipsoid’s Medial axis deformation. Left: Medial axis of an eccentric 3D ellipsoid. Right: Object’s s-rep skeletal sheet.

Let be a smooth curve in . We consider as the unit velocity vector tangent to where is the local tangent plane of at with normal . The local frame can be defined as where (see Figure 3(b)).

Note that in practice, the directions of and may choose two opposite directions because is double-sided. So far, the frame directions are in GCS. To have a consistent frame definition independent of GCS, we design a hierarchical structure.

We start with the m-rep of the ellipsoid’s medial axis and prioritize the ellipsoid’s centroid, medial line, and spokes, respectively. Also, we prioritize points on the medial line closer to the ellipsoid’s centroid and points on the spokes closer to the medial line. Analogous to the m-rep, we give priority to the s-centroid, spine, and veins, respectively. Further, based on the geodesic distance on curves, we prioritize spinal points closer to the s-centroid and points on veins closer to the spine. Thus, given a frame at each skeletal point, we introduce a hierarchical frame structure.

Except for the s-centroid frame, each frame has a prior frame called parent frame. In GP-ds-rep, we have a finite number of prioritized frames on the spine or a vein. A vector that connects a frame to its parent frame is called connection. The tip of a connection is at the frame’s origin, and its tail is at the parent’s origin. Therefore, like a spanning tree, each frame has a parent but may have multiple children. Further, we assume that the s-centroid frame is the parent of itself with the connection .

We approximate the direction of at point based on three consecutive frames. Except for the s-centroid frame and two critical endpoints of the spine that we will explain later, each spinal frame has a spinal parent frame and a spinal child frame. Let and be the position of the parent and the child frame of . As illustrated in Figure 4, assume and as connections. Let and be the projection of and on , respectively. We consider where , and . In this sense, is a unit vector tangent to a circle (or a line) crossing , , and .

The endpoints of the spine are critical because their frames have no children on the spine. By construction, the m-rep medial line is a part of the ellipse’s major axis. After deformation, the major axis becomes a curve we call major curve. The major curve contains the spine and two veins. We consider the closest skeletal point (in geodesic sense) on these veins to the spine as the spine’s extension and treat the critical points as any other spinal point. The s-centroid frame has two spinal children. Let and be the position of the children. We define , where , and . Since a vein frame has a parent and a child on the same vein, we consider the same definition for them as discussed for spinal frames. Note we treat a vein frame at the intersection of a vein and the spine as a spinal frame. For the frames on crest spokes’ tails (i.e., on the skeletal fold), we assume the tip of the crest spokes as the position of the child frames. The same procedure is applicable for the ellipsoid’s GP-ds-rep.

(a) m-rep frame

(b) LP-ds-rep frame
Figure 4: Illustration of a local frames. is normal of tangent planes and . (a) and are equal-length spokes with unit directions and , and (b) is a smooth curve on . and are the projection of and on . , , and .

Figure 5 illustrates the hierarchical structure and a fitted LP-ds-rep to a left hippocampus as described in the next section.

Figure 5: LP-ds-rep. Top: Hierarchical structure of the ellipsoid’s medial axis. Black and blue arrows are connections on the medial line and m-rep spokes. Red dot is the ellipsoid’s centroid. Bottom: A fitted LP-ds-rep to a hippocampus. Red arrows indicate spokes. Black and blue arrows are connections on the spine and veins. Orange arrows depict orthogonal local frames. Red dot is the s-centroid.

2.1.3 LP-ds-rep

Given the fitted hierarchical frame structure introduced in the previous section, we are now in the position to define LP-ds-rep. In an LP-ds-rep, spokes and connections are measured based on their local frames, i.e., their tails are located at the origin of a frame. Let and be the th spoke direction and th connection direction in GCS respectively, where and . Consequently, we denote and as spoke and connection directions based on their local frame, i.e. we re-parameterize and to and respectively. Similarly, if be the frame in GCS then denotes ’s vectors, based on its parent frame.

To calculate a vector direction according to a local frame, we use the spherical rotation matrix , where , and . Therefore, transfers to along the shortest geodesic and we have (amaral2007pivotal).

For example, let frame be the parent of , both in GCS. Let , , and be the axes unit vectors of GCS. We align to such that , where , and . Thus, represents in its parent coordinate system. In case we obtain , we adjust the result by because where . Note that frame vectors are orthogonal, so after rotating to the north pole by , the shortest geodesic between and would be on the equator. This preserve the direction of while rotates .

We follow the the same procedure to calculate the spokes’ and connections’ directions based on their local frames . Finally, a LP-ds-rep is given by , where and are th and th spoke direction and connection direction relative to their local frame with lengths and respectively, is the th frame in its parent coordinate system, and and where is the number of spokes and is the number of frames.

Thus, by construction, the LP-ds-rep is invariant under the act of rigid similarity transformation (i.e., rotation and translation). To remove the scale, we define LP-size as . A scaled LP-ds-rep can be expressed by , where , and . Thus, the LP-size of an scaled LP-ds-re is equal to one. Recall, for a GP-ds-rep, the GP-size is defined as the centroid size e.g., centroid size of spokes’ tails. Note the centroid is an extrinsic property, and the centroid size might be a poor measure for the size of an object. Intuitively, by opening or closing an arm, the arm’s volume remains the same despite its centroid size, i.e., the closed arm has a smaller centroid size as its boundary points are closer to the centroid (see Figure 0(a)).

As Section 2.1.1 discussed, the GP-ds-rep space is . In LP-ds-rep, we do not have any pre-shape space. LP-ds-rep GOPs are directions and lengths of the vectors (i.e., spokes, connections, and frames) plus LP-size. Thus the space is , where is the space of directions and is the space of vectors’ lengths and LP-size. If we euclideanize directions as described in the following sections, the space of euclideanized LP-ds-rep is .

2.2 Population mean

Having a population of ds-reps, suitable methods to calculate means are required in order to perform hypothesis tests on mean differences. The corresponding method should incorporate all geometrical components of the model. Both shape spaces, the GP-ds-rep space, and the LP-ds-rep space are composed of several spheres and a real space. This section will first discuss an approach to analyze the spherical parts by PNS. Afterward, approaches to produce GP-ds-rep means and LP-ds-rep means are discussed.

2.2.1 Pns

PNS (jung2012analysis)estimates the joint probability distribution of data on a d-dimensional sphere by a backward view, i.e., in decreasing dimension. Starting with , PNS fits the best lower-dimensional subsphere in each dimension. A subsphere is called great subsphere if its radius is equal to one; otherwise, it is called small subsphere. For the special case , it is called great circle or small circle

respectively. To choose between the great or small subsphere, we use the Kurtosis test from

(kim2020kurtosis).

PNS is designed for spherical distributions and in particular for small sphere distributions as described in (kim2019small). PNS captures the curviness of circular distributions and euclideanize data as residuals. The PNS residuals on

consist of the geodesic distances between the observations and the fitted circle and the minimal arc length between projected data on the fitted circle to the PNS mean. Therefore in many cases, the distribution of euclideanized data (i.e., residuals) is similar to the bivariate normal distribution.

Figure 7 illustrates fitted circle to a cluster of 1000 observations on , and the PNS residuals. Random points are generated from small sphere distribution (kim2019small) where , , , and .

Alternatively, a simpler but faster euclideanization is to map the data on the tangent space. We transform observations to the north pole by , where is the Fréchet mean. Then, we map the transformed data to the tangent space by the Log map , where , and (fletcher2004principal). Note that for concentrated von Mises-Fisher distribution, the distribution of projected data to the tangent space is close to the distribution of PNS residuals.

Figure 6: PNS euclideanization. Left: Small circle distribution and the fitted circle on . Right: Euclideanizated data.
Figure 7: Mean frame by gradient descent with initial frame . Black dots show the movement of toward the Fréchet means of frames’ vectors.

2.2.2 Mean GP-ds-rep

A method to produce means and shape distributions of a population of GP-ds-reps is composite PNS (CPNS) introduced by (pizer2013nested). The method consists of two steps. First, the two spherical parts of the GP-ds-rep shape space are analyzed by PNS. Spokes’ lengths and scaling factor can be mapped to with the . Afterward, all euclideanized variables are concatenated in addition to some scaling factors that make all variables commensurate. The covariance structure of the resulting matrix is investigated with PCA. Consequently, the mean GP-ds-rep is defined as the origin of the CPNS space. This method depends on a proper pre-alignment and is computationally expensive because PNS has to fit sequential high dimensional sub-spheres to .

2.2.3 Mean LP-ds-rep

To formalize the estimation of LP-ds-rep mean, first we define LP-ds-rep distance. Let and be the geodesic and Euclidean distance, respectively. The distance between two scaled LP-ds-reps , and is given by

(1)

where , .

If be a population of scaled LP-ds-reps then mean LP-ds-rep is

(2)

Note in eq. 1, , and , , . To commensurate with we can use linear shift function instead of . However, this change does not influence the result of eq. 2 because right-hand-side terms of eq. 1 are positive and independent. Assume and let

(3)

By assuming the existence of unique solutions for optimization problems (3), and can be estimated as the Fréchet mean of and , respectively, and and as the arithmetic mean of and , respectively. In this sense, the LP-size of would be equal to one (see Result 1 in SUP). In addition, we need to find mean frames. Calculating mean frame based on each frame vectors may violate the orthogonality condition. Thus the aim is to solve the following optimization problem , —s— n_j,b_j,b_j^⟂∈S^2∑_k=1^N(d_g^2(n_j,n^∗_jk)+d_g^2(b_j,b^∗_jk)+d_g^2(b_j^⟂,b^∗⟂_jk))^12, d_g(n_j,b_j)=d_g(b_j,b_j^⟂)=π2. Let where , , and are Fréchet means of , , and respectively. may or may not belong to (i.e., is not necessarily a frame). Thus is the answer of the optimization (2.2.3) without its constrain. We rotate the initial frame to be as close as possible to component wise (see Figure 7). Since rotation preserves the frame orthogonality, the aligned frame approximates the solution. For this purpose we use Algorithm 1

as a gradient descent approach. Note that the PNS mean estimates mean direction based on the euclideanized data, so we consider PNS mean instead of Fréchet mean. Also, we use the geometric mean instead of the arithmetic mean.

To reduce convergence time, choosing an appropriate initial frame for Algorithm 1 is essential. Analogous to GPA, we can consider the initial frame as , i.e., is rotated such that its unit centroid coincides with the unit centroid of . Alternatively, since in practice , and based on the fact that we defined as the cross product of and , we accelerate the convergence by defining the initial frame based on and . Assume the middle point of the shortest geodesic connecting and . We move in opposite directions from toward and by angle to reach two points. We consider these points as and because, , and they have equal distance to and (see SUP).

1.1

Algorithm 1 Frame alignment.
0:  Unit vectors , , , initial frame (e.g., ), step size , and threshold
0:  Aligned frame to ,,
  
  while  do
     for  to 3 do
        
     end for
     
  end while

( is the tangent vector at pointing toward , , and is the spherical rotation matrix.)

2.3 Convert LP-ds-rep to GP-ds-rep

Section 2.1.2 and Section 2.1.3 discuss how to obtain an LP-ds-rep from a GP-ds-rep. For several reasons, e.g., for visualization, we may need to reverse the procedure. For GP-ds-rep visualization, it is sufficient to draw spokes individually. To visualize an LP-ds-rep, we convert it to a GP-ds-rep. We start from as the s-centroid frame. Then, we reconstruct frames by finding the position and orientation of the frame’s children based on . Afterward, we find the information of grandchildren frames based on their parents and so on.

Let frame be in the coordinate system of its parent . To find based on GCS, we rotate by such that . Then is the representation of in GCS. Similarly, we find the direction of connections and spokes in GCS.

Finding the mean shape of a set of objects’ boundaries without an alignment is almost impossible. But we can use LP-ds-rep to estimate the mean boundary without alignment. First, we calculate the mean LP-ds-rep. Then, we convert the mean LP-ds-rep to a GP-ds-rep. Finally, we generate the implied boundary from the GP-ds-rep as demonstrated in (liu2021fitting). Therefore, it is possible to approximate the mean boundary without the alignment, which shows the power of LP-ds-rep.

3 Hypothesis testing

Let and be two groups of either GP-ds-reps or LP-ds-reps of sizes and . Let be the total number of GOPs. To test GOPs’ mean difference, we design partial tests. Let and be the observed sample mean of the th GOP from and respectively. The partial test is versus . Note for GP-ds-rep , and for LP-ds-rep . Also, for GP-ds-rep tests we pre-align the pooled group by GPA.

To test mean differences, we adapted a non-parametric permutation test with minimal assumptions similar to Styner’s approach (Styner2006)

. For the univariate data i.e., vectors’ lengths and shapes’ sizes, the test statistic is t-statistic

where

is the pooled standard deviation. For the multivariate data i.e., euclideanized directions and GP-ds-rep skeletal positions, the test statistic is Hotelling’s T

metric , where

is an unbiased estimate of common covariance matrix

(mardia1982multivariate, ch.3). Given the pooled group , the permutation method randomly partitions times the pooled group into two paired groups of sizes and without replacement, where usually we consider . Afterward, it measures the test statistic between the paired groups. The empirical -value for the th GOP is , where is the th observed test statistics, is the th permutation test statistic, and is the indicator function i.e., if is true, otherwise .

In order to account for the problem of multiple hypothesis testing, one could use Bonferroni’s method (bonferroni1936teoria). Bonferroni’s method tests each hypothesis at level

and guarantees the probability of at least one type I error

be less than significance level . Since Bonferroni’s method is highly conservative we prefer to apply Benjamini-Hochberg (BH) (benjamini1995controlling) as a more moderate approach.

4 Evaluation

4.1 Data

To test our method, we study the hippocampal difference between early PD and CG at baseline. Data are provided by ParkWest (http://parkvest.no), in cooperation with Stavanger University Hospital (https://helse-stavanger.no). At the baseline, we have 182 magnetic resonance (MR) images for PD and 108 for CG with corresponding segmentation of hippocampi. As described in Section 2, GP-ds-reps are fitted to left hippocampi by SlicerSALT (http://salt.slicer.org) and re-parametrize into LP-ds-reps. For the model fitting, we used GP-ds-reps with 122 spokes consisting of 51 up, 51 down, and 20 crest spokes. As up and down spokes share the same tail position, we have in total 71 tail positions. Thus, for LP-ds-rep, we have 122 spokes, 71 local frames, and 71 connections. Before analyzing the Parkinson data, we first study our method based on simulations.

4.2 Simulation

In statistical shape analysis generating random shapes is a matter of interest. Designing simulation based on GP-ds-rep is challenging as we usually need to identify a local frame to bend or twist the object locally. It turned out that LP-ds-rep support naturally skeletal deformations. We can stretch, shrink, bend, and twist the skeletal by manipulating the frames’ orientations and vectors’ lengths. Then, we convert LP-ds-rep to GP-ds-rep to generate the boundary. Consequently, we can add variation to a set of deformed LP-ds-reps’ GOPs to simulate random ds-reps. Figure 8 shows a deformed hippocampus including bending and twisting. The deformation is done by rotating six spinal frames.

For the simulation study, we select the mean LP-ds-rep of CG from the ParkWest data as a template. Based on the template, we generate two LP-ds-rep groups of sizes 150 with different amount of tail bending, i.e. bending in local region. Such bending was observed for example in (pizer2003object) between schizophrenics and controls. Let denotes von Mises-Fisher distribution with mean and concentration parameter on (dhillon2003modeling). For the special case we assume the distribution in radian i.e., if . Given random rotation angle of bending for the first group and for second group, we simulate the orientation of three spinal frames by successively rotating them about their -axis with . This means, the tails in the second group is in average successively bend downward for three consecutive spinal frames. Chosen frames are the closest ones on the hippocampus tail to the s-centroid. Thus, in total, we have a slight downward bending about at the hippocampus tail. Finally, by preserving frame orthogonality, we add noise to all directions by , where for frames’ vectors, spokes, and connections is equal to , , and , respectively. Further we added noise to vectors’ lengths by the truncated normal distribution where is the vector length of the template, and parameters , , and

are heuristically chosen. As a result, we have two groups of random LP-ds-reps, which are approximately similar in most of their GOPs but only different in the orientation of three frames.

Figure 9 illustrates twenty samples of each group in blue and red. Note that LP-ds-reps are not aligned, but since we reconstruct them from the s-centroid frame, shapes have Bookstein’s alignment (dryden2016statistical, Ch. 2) because the s-centroid frames are perfectly aligned.

Figure 8: Skeletal deformation by LP-ds-rep. Left: A ds-rep with its implied boundary in two angles. Middle: Shape bending by spinal frame rotation about and axes. Right: Shape twisting by spinal frames rotation about axis.

Figure 9: ds-rep simulation by LP-ds-rep. Left: Blue and red indicate twenty samples of two groups of simulated ds-reps. Right: Overlaid mean LP-ds-reps.

Hypothesis test on LP-ds-rep from Section 3 correctly detects significant frame directions and label almost all other GOPs as statistically non-significant given a significance level . On the contrary, the test on GP-ds-rep indicates a large number of false positives, i.e., almost all of the positions and directions are statistically significant (see Figure 10). This example confirms our observation from Figure 1 in Section 1 and highlights the fact that GP-ds-rep analysis could be extremely biased.

(a) LP-ds-rep -values
(b) GP-ds-rep -values
(c) Significant frames
Figure 10: LP-ds-rep vs. GP-ds-rep. (a,b) Raw and adjusted -values by BH and Bonferroni. The dotted line indicates significance level . (c) Blue arrows indicate local frames. Red indicates statistically significant frame directions after the BH adjustment using LP-ds-rep analysis.

4.3 Real data analysis

The Parkinson data set described in Section 4.1 was studied earlier by (apostolova2012hippocampal) based on radial distance analysis and parallel slicing and showed some regional atrophy. Since shape correspondence in parallel slicing is controversial, we attempt to reanalyze data by utilizing LP-ds-rep.

First let us compare the shape sizes, see Table 1. Tests on shape size indicate no significant difference. However, volume measurement confirms the LP-size is more compatible with the object volume as for both, the mean object volume and the LP-size of CG are greater than PD. In opposite the mean GP-size of CG is smaller than PD.

Mean CG Mean PD SD CG SD PD p-value
GP-size of spokes’ tips 161.05 162.51 8.97 8.62 0.17
Object volume (mm) 3352.23 3271.44 563.39 616.68 0.26
LP-size 536.31 527.60 36.09 38.86 0.06
Table 1: T-test on shape size.

Figure 11: ds-rep significant GOPs. Red indicate significant GOPs. FDR=0.05 for BH adjustment.

Figure 11 illustrates significant LP-ds-rep and GP-ds-rep GOPs before and after BH adjustment in red. In LP-ds-rep, all the spokes directions are insignificant. In contrast, about 40% of GP-ds-rep spokes’ directions are significant. Also, in LP-ds-rep, there are a few significant connection and frame directions after the adjustment. Based on the LP-ds-rep analysis, it seems the main difference comes from connections’ length on the spine. Figure 11(a) and Figure 11(b) show sorted -values before and after adjustment. Based on Bonferroni adjustment, PD and CG are similar because almost all adjusted -values are greater than 0.05. But based on raw and BH -values, about half of the GP-ds-rep GOPs are significant, while in LP-ds-rep, we have only a few.

(a) LP-ds-rep -values
(b) GP-ds-rep -values
(c) Scaling effect
Figure 12: Test on real data. (a,b) Raw and adjusted -values by BH and Bonferroni. The dotted line indicates significance level 0.05. (c) Bar plot shows the scaling effect on the percentage of .

In addition, we analyzed the shapes without scaling. Detailed results are available in SUP. The general belief is scaling makes shapes more similar. But Figure 11(c) expresses the percentage of significant GP-ds-rep GOPs increases dramatically after the scaling. In other words, scaling increases the number of raw -values less than the level of significance and consequently increases the number of BH adjusted -values less than FDR=0.05. A possible explanation is that GPA tries to make shapes as close as possible by reducing GOPs’ variation. By removing the scale, GPA reduces the variation even more. Hotelling’s T metric is proportional to the inverse common covariance matrix. So by reducing the variation, the test statistic increases, and consequently, the -value decreases. On the contrary, the LP-ds-rep is not sensitive to scaling as it is alignment-independent.

5 Conclusion

Generally, it is common to detect locational dissimilarity between two groups of objects based on the alignment. We showed that alignment-dependent analysis such as ds-rep analysis in GCS could be highly biased. We described GP-ds-rep as a conventional ds-rep parameterization in GCS. We introduced LP-ds-rep as a novel parameterization with a hierarchical structure based on locally fitted frames to overcome inherent challenges caused by alignment. We have defined a mean LP-ds-rep and discussed object deformation and simulation with LP-ds-rep. We explained how to estimate boundary mean shape without alignment. We compared LP-ds-rep with GP-ds-rep to show the advantages of LP-ds-rep. For comparison, we applied simulation and real data analysis. The simulation confirmed that the hypothesis test based on GP-ds-rep for two groups of ds-rep with a slight difference results in many false positives while LP-ds-rep indeed detected the differences. For the real data, we studied left hippocampi of early PD vs. CG. Although hypothesis tests on GP-ds-rep indicated many significant GOPs, tests on LP-ds-rep showed only a few, which seems medically more reasonable. Also, data analysis exposed GP-ds-rep sensitivity to the scaling. We concluded that PD and CG groups are very similar, but the main difference comes from the spinal stretch of the skeletal sheet.

Acknowledgments

This research is funded by the Department of Mathematics and Physics of the University of Stavanger (UiS). Special thanks to Profs. Stephen M. Pizer (UNC), Steve Maron (UNC), James Damon (UNC), and Jan Terje Kvaløy (UiS) for insightful discussions and inspiration for this work. We are indebted to Prof. Guido Alves (UiS) for providing ParkWest data. We also thank Zhiyuan Liu (UNC) for the model fitting toolbox.

SUPPLEMENTARY MATERIALS

Supplementary: SUP materials referenced in this work are available in this pdf. (pdf) R-code: In Supplementary.zip, simulation codes and files are placed. (zip)

References