Topological Data Analysis and Persistent Homology

Persistent homology is a new area of computational topology that aims to understand the underlying high-dimensional structures from low-dimensional local topological structures. Instead of analyzing shapes and images at a fixed scale, as usually done in traditional approaches, persistent homology observes the changes of topology over different scales and finds the most persistent topological features that are robust under noise perturbations. By 2009, there had been many theoretical advancements in the field, but large-scale real-world applications beyond toy examples to clearly demonstrate its usefulness and power were lacking. Thus, the concept had been a theoretical curiosity in medical imaging fields for many years. We wrote the first paper on the application of persistent homology to MRI using the cortical thickness of the human brain (Chung et al., 2009).

MATLAB codes:
Exact Topological Inference for Network Differences
Persistent and Min-Max Diagrams for Cortical Thickness
Betti Plots for Brain Networks
Topological Network Distances

Large-scale Brain Network Analysis

Diffusion tensor imaging (DTI) regularly produces up to a half million white matter fiber tracts per brain, which correspond to a half million edges and one million nodes in a network. The resulting networks defy many standard network methods due to the sheer computational load. The connectivity matrix corresponding to such data would contain a hundred billion connections, which results in the small-n large-p problem (Chung et al. 2015). Specifically, when the number of voxels (p) are substantially larger than the number of images (n), it produces an under-determined model with infinitely many possible solutions. Our group was first to develop the compressed sensing (CS) technique to thin out the significant number of noise connections in the brain network (Lee et al., 2011).

Any CS or sparse brain network model  is usually parameterized by a tuning parameter  that controls the sparsity of the representation. Increasing the sparse parameter makes the  solution more sparse. Thus, sparse models are inherently multiscale, where the scale of the model is determined by the sparsity. Many existing sparse network models use a fixed sparse parameter that may not be optimal in other datasets or studies. Depending on the choice of the sparse parameter, the final classification and statistical results will be different (Chung et al., 2015). Therefore, there is a need to develop a multiscale network model that provide a consistent analysis results and interpretation regardless of the choice of parameter. Persistent homology offers one possible solution to the multiscale problem. Instead of studying images and networks at a fixed scale, as usually done in traditional approaches, persistent homology summarizes the changes of topological features over different scales and finds the most persistent topological features that are robust under  perturbations. This robust performance under different scales is needed for network models that are parameter and scale dependent. Instead of building networks at one fixed parameter that may not be optimal, I have developed a new persistent homological framework for analyzing the collection of sparse models over every possible sparse parameter (Chung et al., 2015; Chung et al., 2017). The method exploits the topological structure in a family of sparse models in drastically speeding up the computation for building large-scale brain networks with billions of connections.

MATLAB codes:
Transposition Test for Large-scale Network Data
Structural Brain Connectivity from Diffusion Tensor Images
Heritability of Large-scale Brain Networks
 Epsilon Neighbor Structural Brain Network Construction for DTI
Betti Plots for Brain Networks
Sparse Hyper-Network Model
Topological Network Distances

Heat Kernel Smoothing & Diffusion on Manifolds

A lot of anatomical features that characterize brain shape variability can only be measured relative to a cortical surface. It is necessary to develop a surface-based framework for increased sensitivity and specificity of analysis. A significant hurdle in this area had been caused by the lack of surface-based  techniques that can incorporate the non-Euclidean geometric nature of the cortical surface. The concept of smoothing anatomical measurements along the brain tissue boundaries for this purpose (Chung et al. 2001;  Chung et al. 2003). Diffusion smoothing, a concept introduced in Chung et al. (2001), utilizes  isotropic heat flow as a way to smooth brain surface measurements such as cortical thickness, curvatures or functional MRI. The drawbacks of diffusion smoothing are the complexity of numerically solving the PDE and the numerical instability that depends on the smoothness of anatomical boundary. To address the issues, heat kernel smoothing was  introduced (Chung et al., 2005). In heat kernel smoothing, the analytically unknown nonlinear heat kernel is approximated linearly and smoothing is performed as the iterative convolution in a spatially adaptive fashion. These two methods are the standard baseline methods in  building more complex cortical surface models in the field. The methods  can easily detect surface specific tissue growth and atrophy that are more difficult to detect using the traditional 3D volume based methods.

MATLAB codes:
Polynomial Approximation to Heat Diffusion on Manifolds 
  Heat Kernel Smoothing: Iterated Kernel Smoothing Version
Heat kernel smoothing via the Laplace-Beltrami eigenfunctions
Diffusion smoothing on manifolds

Weighted Fourier Series Representation

Weighted Fourier Series (WFS) representation incorporates parameterization, image smoothing, surface registration and statistical inference in a single mathematical framework while reducing the Gibbs phenomenon (ringing artifacts) associated with the traditional Fourier analysis (Chung et al., 2007). In WFS, signals are represented as a weighted linear combination of eigenfunctions of the self-adjoint operator of the underlying manifold. The weights are exponentially decaying and related to the heat kernel of the manifold. On a sphere, the eigenfunctions are spherical harmonics. In the Euclidean space, they are the usual Fourier basis.  We have  recently shown that WFS is in fact mathematically equivalent to the diffusion wavelet transform (Chung et al., 2015). Thus, WFS inherits all the nice localization properties of wavelets. As a special case of WFS, we have developed weighted spherical harmonic (SPHARM) representation. The iterative residual fitting (IRF) algorithm was proposed for estimating more than 20000 spherical harmonic coefficients per brain (Chung et al., 2007) The algorithm sequentially breaks down one gigantic problem into smaller problems in an iterative fashion.

MATLAB codes:
Weighted Spherical Harmonic Representation (Weighed-SPHARM)
Cosine Series Representation of 3D Curves and White Matter Fiber Tracts
Amygdala/hippocampus surface modeling using SPHARM

Tensor-Based & Surface-Based Morphometry

In computational neuroanatomy,  region-of-interest (ROI) free MRI morphometric techniques mainly referred to as tensor-based morphometry (TBM) and deformation-based morphometry (DBM) are popular. In trying to model shape variations in medical images, it is not clear what sort of image features should be analyzed. Building on the framework of fluid dynamics, Chung et al. (2001) modeled and analyzed the distribution of the Jacobian determinant, divergences and vorticities of brain deformation systematically in a unified statistical framework. These quantities are essential in measuring the tissue growth and atrophy at each voxel. Chung et al. (2001) demonstrated the Jacobian determinant is the only meaningful biomarker in tissue growth and atrophy. The Jacobian determinant approach has been used in a wide variety of applications including the localization of the region of brain tissue growth before and after cochlear implant in positron emission tomography (PET) images. The later studies  (Chung et al. 2003) further extended the idea of the Jacobian determinant in the Euclidean space to manifolds in a more generalized tensor formulation and were able to surface specific growth locally in developing children. Currently, the Jacobian determinant is the standard baseline feature for characterizing tissue growth and atrophy in brain imaging.

MATLAB codes:
Amydala Surface Modeling Using Keith Worsley's SurfStat
Mandible Surface Modeling Using the Laplace-Beltrami Eigenfunctions
Voxel-Based Morphometry to Modeling Binary Stroke Lesions in DIffusion Weighted Images

Other MATLAB packages

General linear models (GLM) in imaging

1. How to read huge analyze file into MATLAB
2. Brain surface mesh manipulation tools 

How to read ImageJ ROI file format into MATLAB?
Image format conversion Tips

Cortical mesh MNI format
Corrected p-value on cortex
Functional data analysis on a unit circle
Deformation field manipulation codes

177.