Spatiotemporal Statistical Shape Model Construction for the Observation of Temporal Change in Human Brain Shape

Abstract

This chapter introduces a spatiotemporal statistical shape model (stSSM) using brain MR image which will represent not only the statistical variability of shape but also a temporal change of the statistical variance with time. The proposed method applies expectation-maximization (EM)-based weighted principal component analysis (WPCA) using a temporal weight function, where E-step estimates Eigenvalues of every data using temporal Eigenvectors, and M-step updates Eigenvectors to maximize the variance. The method constructs stSSM whose Eigenvectors change with time. By assigning a predefined weight parameter for each subject according to subjects’ age, it calculates the weighted variance for time-specific stSSM. To validate the method, this study employed 105 adult subjects (age: 30–84 years old with mean ± SD = 60.61 ± 16.97) from OASIS database. stSSM constructed for time point 40–80 with a step of 2. The proposed method allows the characterization of typical deformation patterns and subject-specific shape changes in repeated time-series observations of several subjects where the modeling performance was observed by optimizing variance.

1. Introduction

Quantifying cortical morphological dynamics of brain deformation will help neuroscientists identifying and characterizing brain deformation disorders. To be precise, if physicians could learn to predict the normal cortical shape evolution for healthy adults, they can predict abnormal or early deformation of brain with time as well. The shape variability from statistical shape models has been successfully utilized to perform various deformation-based researches in the field of image analysis in both two-dimensional (2D) and three-dimensional (3D) images. In specific cases, their application for image segmentation in the context of statistical shape models has been well acknowledged. To constract those statistical models, a set of segmentations of the shape of interest is required. At the same time, a set of feature vectors are also needed which can be decidedly defined in each sample shape. There are a number of researchers who have used statistical models containing shape information as initials for segmentation via deformable models [1, 2, 3], deformable registration [4], or shape analysis [5]. Meanwhile, the fundamental problem faced at the time of constructing these models is the reality that they need the determination of point correspondences between the different shapes as the manual identification of such correspondences is a time-consuming and tiresome work. It is specially applicable in 3D where the number of feature vectors required to describe the shape accurately rises dramatically compared to 2D applications. An easy but efficient way of handling this issue for the construction of statistical shape models is the use of distance transformations which has been proposed by a number of authors [2, 6, 7]. For example, in [2], the authors proposed an approach in which the statistical analysis is carried out directly on the signed distance maps of a set of aligned shapes.

However, learning predictive models to trace forth the evolution trajectories of adult cortical shapes remains challenging. There are some studies to analyze the brain deformation in adults. For example, they are studying nonrigid brain shape registration with respect to baseline distribution, brain region segmentation based on fuzzy object model, sulcal curves extraction on the outer cortex, etc. In [8], authors described an automated way in which correspondences between the surfaces of different shapes are established via a nonrigid registration algorithm [9]. A similar approach has been proposed by [10]: there, a modified iterative closest point algorithm [11] is used to calculate correspondences between geometric surface features of the shapes such as crest lines. The difficulty of analyzing the adult brain shape points to the critical selection of parameter optimization.

Recently, Durrleman et al. [12] have proposed the construction method of spatiotemporal statistical shape model but the study is confined to longitudinal data. To fill this critical gap, we propose spatiotemporal statistical shape model (stSSM) for temporal 3D shape change analysis of the adult brain. We are performing a dimensionality reduction analysis directly on a parametric representation of the feature vector calculated from training sample. To match different anatomies, all subjects are aligned using anterior-commissure (AC), posterior-commissure (PC), and mid-plane for location and orientation alignment with a reference of AC as origin. This enables the construction of average models of the brain shape and their statistical variability across a population of subjects. Our approach is closely related to the statistical shape models first proposed by [13, 14], but differs in an important aspect. Rather than performing a classic principal component analysis (PCA) directly, we are creating an stSSM for temporal 3D shape change analysis of the adult brain with respect to cortical surfaces using an expectation-maximization (EM)-based weighted PCA (WPCA) learning framework where the weight function is defined as a Gaussian function, whose center is time point, and variance is a predefined parameter.

2. Preliminary

Basically, SSM construction is a method of extracting the average shape and a number of modes of variation from a collection of training samples. The methods are strongly dependent on the chosen shape of representation. A vital necessity for building shape models with statistical variability is that the features of all training samples need to be located at corresponding positions and sharing same origin. With an assumption of a one-to-one correspondence of anatomical structures across subjects, the alignment of images between different subjects with chamfer distance calculation provides a dense set of correspondences. In this study, we performed an expectation-maximization algorithm-based wPCA on a compact parameterization of the feature vector which is required to construct our proposed model. In the following subsections, we will describe the proposed approach in detail.

2.1. Distance mask calculation

For distance mask calculation, level sets were introduced by [15] and made popular for computer vision and image analysis by [16]. They basically feature an implicit shape representation and can be utilized with regional or edge-based features. The authors of [17] presented a method of embedding the distance maps into the linear space, which could solve the modeling problems and Cremers et al. [18] have given an overview of statistical approaches to level set segmentation including prior shape knowledge.

In order to construct stSSM, firstly, brain region was segmented from MR images. The segmentation is done automatically using FSL [19]. In the next step, chamfer distances were calculated to extract shape features through a level set algorithm where negative values are assigned inside the brain region and positive value outside the brain region.

2.2. Dimensionality reduction

After calculation of chamfer distance, the next step is to reduce the dimensionality of the training set, that is, to find a small set of modes that best describes the observed variation. This is usually accomplished using principal component analysis (PCA) [20]. PCA is a very powerful tool to analyze data by creating a custom set of “principal component” eigenvectors that are optimized to describe the most data variance with the fewest number of components. Theoretically, the steps are simple: the principal components {Φ_k} of a dataset are simply the eigenvectors of the covariance of that dataset, sorted by their descending eigenvalues which results to a new observation, C = μ + ∑ᴪ_kΦ_k where μ is the mean of the initial dataset and ᴪ_i is the reconstruction coefficient for eigenvector Φ_i [21]. One of the limitation of classic PCA is that it does not distinguish between variance due to measurement noise vs. variance due to genuine underlying signal variations. Even when an estimation of the measurement variance is available, this information is not included while constructing the eigenvectors [22]. In order to overcome this problem, Bailey [22] introduces PCA based on EM-PCA.

In our method, we used weighted PCA with EM-PCA. EM is an iterative technique for solving parameters to maximize a likelihood function for models with unknown latent variables which basically consists of E-step and M-step. E-step estimates eigenvalues of every data using temporal eigenvectors, and M-step updates eigenvectors so that it maximizes the variance. E- and M-steps are iterated until convergence of updating vectors. The detail is described below.

[E-step] Let the current eigenvector be $\overrightarrow{\mathrm{\Phi }}$ . Principal component (PC) score of data j is calculated by:

$$c_j=\overrightarrow{x_j}\cdot \overrightarrow{\mathrm{\Phi }}$$ E1

where $\overrightarrow{x_j}$ is the feature vector of data j.

[M-Step] Update the eigenvector by maximizing the weighted variance of PC score. It is calculated by:

$$\overrightarrow{\mathrm{\Phi }}={\sum }_j{w}_j{c}_j\overrightarrow{x_j}$$ E2

The eigenvectors are normalized after each M-step. E- and M-steps are iterated till the updated value of eigenvector is converged.

3. Method

To construct stSSM, firstly, brain region was segmented from MR images, and chamfer distances were calculated to extract shape features as shown in Figure 1.

Figure 1.

Preprocessing steps of brain MR image for feature extraction. (a) Raw MR image, (b) Segmented image, and (c) Signed distance.

After that, weighted PCA is applied to the shape features. The proposed method assigns a weight parameter for each subject according to subjects’ age, and calculates the weighted variance. Let t_i be a time point, and the method can construct SSM at the time point. The weight function was defined as a Gaussian function, whose center is t_i, and variance is a predefined parameter.

$$w_j=\text{exp}\left(-\frac{{\left(t-t_j\right)}^2}{{\sigma }^2}\right)$$ E3

The shape of temporal weight function is illustrated in Figure 2. That is, subjects near t_i are dominant to decide the Eigenvectors. By shifting the t_i at a short interval from the minimum to the maximum age, the method constructs the stSSM whose eigenvectors change with growing or aging. Preliminarily, the eigenvector at the first time point is initialized randomly. Then, the obtained eigenvectors are used at the stSSM construction.

Figure 2.

Temporal weight function.

4. Experimental results

For evaluation purpose, the proposed method has been applied on both artificial and image data.

4.1. Result of artificial data

In order to numerically evaluate the method, weighted PCA has been applied to point distribution model. Assumed a point of an organ, and the statistical distribution changes temporally. Figure 3 shows the artificially generated 2D data where the points located one direction shown at the first time point (a), and rotates anticlockwise with 10-degree every time point. The time period was 20. Because, conventional SSM evaluates all data simultaneously without time, as shown in Figure 3(c). It cannot be statistical variation of shape.

Figure 3.

Temporal change point distribution. (a) t = 1, (b) t = 5, and (c) all data.

Figure 4 shows the angle between x-axis and the first (second) principal axis obtained by the proposed method. The full-width-half-maximum (FWHM) of the weight Gaussian function was 3. The first principal axis rotates anticlockwise with 10-degree every time point. Results show that the proposed method can extract the statistical variability of point distribution in time. Note that the result around t = 1 has errors because the data existed only in the positive side of Gaussian function.

Figure 4.

Temporal change of principal axes.

4.2. Image data

To train and test the proposed model with respect to weighted variance optimization, we used 105 adult subjects (age: 30–84 years old with mean ± SD = 60.61 ± 16.97) from publicly available imaging database called OASIS [23]. These subjects were selected from a larger database of individuals who had participated in MRI studies at Washington University, were all right-handed, and older adults had a recent clinical evaluation. The representative MR imaging acquisition parameters were repetition time (TR) of 9.7 ms, echo time (TE) of 4.0 ms, image resolution (voxel) of 256 by 256, flip angle (FA) of 10°.

4.3. Results

The method was applied to the adult brain MR image. We first segmented the brain region from MR images, and calculated the signed distance map to extract shape features. The proposed method assigns a weight parameter for each subject according to subjects’ age, and calculates the weighted variance. Let t_i be a time point, and the method can construct SSM at the time point. According to the weight function defined in Eq. (3), the center is t_i, and variance is a predefined parameter.

From Figure 2, we can understand the shape of temporal weight function which implies that subjects near t_i are dominant to decide the Eigenvectors. By shifting the t_i at a short interval from the minimum to the maximum age, the method constructs the stSSM whose eigenvectors change with aging. Preliminarily, the eigenvector at the first time point is initialized randomly (usual PCA can be applied to obtain the initial eigenvector). Then, the obtained eigenvectors are used at the following stSSM construction. Figure 5 shows the mean shape of constructed model between age range of 40–80 years old. Figure 6 illustrates stSSM of brain of adult subject at 60 years. The brain shape changes along with temporal domain with reference to variance (σ) can be seen.

Figure 5.

stSSM mean shape (a) 40 years (b) 50 years (c) 60 years (d) 70 years, and (e) 80 years.

Figure 6.

Spatiotemporal statistical shape model of adult brain. The horizontal and vertical axes are the first and the second principal axes.

Next, the stSSM of adult brain between 52 and 58 years and temporal difference with changing variance is shown in Figure 7. While generating stSSM for an age range of 40–80, Gaussian function’s weight parameter has been varied between a range of 5–35.

Figure 7.

Temporal difference with changing variance.

Figure 8 illustrates stSSM of adult subject brain shape at age of 60 years old, where the individual brain shape variety is synthesized within a range of −1σ to +1σ at the first eigenvector (horizontal) and −0.5σ to +0.5σ at the second eigenvector (vertical) with a step of 0.50. σ is the standard deviation (SD) in the training data along each eigenvector. Full-width-half-maximum (FWHM) of the Gaussian weight function used in this case was 15.

Figure 8.

Spatiotemporal statistical shape model of adult brain. The horizontal and vertical axes are the first and the second eigenvectors. The individual brain shape variety is synthesized within a range of −1σ to +1σ at the first eigenvector (horizontal) and −0.50σ to +0.50σ at the second eigenvector (vertical) with a step of 0.50. σ is the standard deviation in the training data along each eigenvector.

5. Discussion

The 3D SSM and stSSM in the medical imaging field are almost exclusively based on imaging modalities such as CT, MRI, that is, the original data representation of the training shapes is not a mesh but rather a segmented volume. Therefore, each shape of the training set must be annotated by features, each designating the same anatomical locus along the set. The set is considered as a collection of shape vectors which after alignment raises a covariance matrix. Figures 1 and 2 show these phases of our model from which we created feature vector matrix from 105 training samples.

For dimensionality reduction phase, we introduced weighted EM-PCA where every aligned training shape is described by the feature vector matrix. The mean shape has been formed by simply averaging over all samples with corresponding weight parameter where weight has been assigned by subject’s age. From age range of 40–80 years with a step of 2, stSSM model in 3D has been constructed. An eigen decomposition on covariance matrix gives principal modes of variation (eigenvectors) and their respective variances (eigen values). Figure 5 visualizes the mean shape of 40, 50, 60, 70, and 80 years. From the computational point of view, this method is used due to higher numerical stability. The resulting modes of variation are ordered by their variances. In a next step, the model approximated every valid shape by a linear combination of the chosen accumulated variance which presented the shape variation with FWHM Gaussian distribution function. The shape variation at a specific time point is visualized in Figure 6. The horizontal and vertical axes are the first and the second principal axes.

Difference between temporal deformations of all parameters has been shown in Figure 7. We calculated the temporal deformation change between two time points and evaluated the corresponding difference. Positive value has been assigned inside, while negative represents outside with a view that higher temporal deformation difference will show high contrast. Our study shows that when FWHM of Gaussian function is set to a value smaller than 15, it shows larger difference in comparison with FWHM higher than 15 where it is not represented. This observation leads to an optimum value of weight parameter between 10 and 15. Constructed stSSM for specific age is shown in Figure 8.

The proposed method was evaluated by generalization ability according to leave-one-out cross validation (LOOCV) procedure. Generalization ability evaluates the performance of representing new acceptable models. First, data of the evaluation subject were projected into stSSM of the evaluation subject’s age. Next, brain shape of evaluation data was reconstructed, and compared with the original data. Evaluation criteria used in this study was Jaccard Index (J.I.). To restrict the allowed variation to plausible shapes, weight parameter is assigned to a certain interval. Mean ± standard deviation of generalization ability with different Gaussian distribution function has also calculated. Overall, FWHM value 15 achieves the best prediction result which closely follows observation-based optimum value of weight parameter between 10 and 15.

6. Conclusion

Nowadays, statistical shape models have become an exciting robust tool for shape representation of medical images. The use of both 2D and 3D models appeared into the scenario in recent years. In this chapter, we introduced a 3D stSSM for brain MRI data. The presented method allows the characterization of typical deformation patterns and subject-specific shape changes in repeated time-series observations of several subjects. The modeling performance was observed by optimizing variance. This study can be seen as an extension of the usual statistical shape model of scalar measurements to high-dimensional shape or image data. From the analysis of difference between all parameters of temporal deformation, we can draw a conclusion toward optimum value of FWHM.

We believe that the discussed method of this chapter enables the automatic construction of statistical models in 3D and is not limited to the brain but can be applied to other anatomical structures such as the heart or liver. The spatiotemporal statistical shape model has a wide number of possible applications primarily in segmentation and morphometry. Another potential application of this method lies in the use of statistical modes of variation as a priori knowledge for image registration. We are currently working on this idea which uses the modes of variation from predefined weight parameter to provide a more compact parameterization of stSSMs. This may be specifically useful for intersubject analysis tasks where these modes of variation can be learnt from a sample population of subjects as shown in this chapter.

Another potential scope is the morphometric comparison of differences between groups of subjects. Currently available morphometric methods can be classified into voxel-based [24, 25] or deformation-based methods [26, 27]. Mostly, voxel-based methods depend on a global registration between subjects followed by a statistical analysis of tissue differences to differentiate between groups of subjects. On the other hand, deformation-based methods use the information encoded in the deformation to describe the anatomical variability between groups. Future work will include investigation whether statistical shape models can be used as a deformation-based morphometric tool to characterize shape differences between groups of normal and AD subjects.

Ethical approval

All applicable international, national, and/or institutional guidelines for the care and use of animals were followed. This study does not contain uninformed patient data.

Acknowledgments

This work was supported in part by JSPS Grant-in-Aid for Scientific Research on Innovative Areas (Multidisciplinary Computational Anatomy) JSPS KAKENHI Grant Number 17H05304.

Conflict of interest

The authors of this chapter declare that they have no conflict of interest.