Brain white matter hyperintensity lesion characterization in 3D T2 fluid-attenuated inversion recovery magnetic resonance images: Shape, texture, and their correlations with potential growth

Analyses of age-related white matter hyperintensity (WMH) lesions manifested in T₂ fluid-attenuated inversion recovery (FLAIR) magnetic resonance images (MRI) have been mostly on understanding the size and location of the WMH lesions and rarely on the morphological characterization of the lesions. This work extends our prior analyses of the morphological characteristics and texture of WMH from 2D to 3D based on 3D T₂ FLAIR images. 3D Zernike transformation was used to characterize WMH shape; a fuzzy logic method was used to characterize the lesion texture. We then clustered 3D WMH lesions into groups based on their 3D shape and texture features. A potential growth index (PGI) to assess dynamic changes in WMH lesions was developed based on the image texture features of the WMH lesion penumbra. WMH lesions with various sizes were segmented from brain images of 32 cognitively normal older adults. The WMH lesions were divided into two groups based on their size. Analyses of Variance (ANOVAs) showed significant differences in PGI among WMH shape clusters (P = 1.57 × 10^–3 for small lesions; P = 3.14 × 10^–2 for large lesions). Significant differences in PGI were also found among WMH texture group clusters (P = 1.79 × 10^–6). In conclusion, we presented a novel approach to characterize the morphology of 3D WMH lesions and explored the potential to assess the dynamic morphological changes of WMH lesions using PGI.

Introduction

White matter hyperintensities (WMH) on T₂ fluid-attenuated inversion recovery (FLAIR) magnetic resonance brain images (MRI) are commonly observed in older adults over 65 years old with a prevalence rate of ∼ 60–80% in the general population (De Leeuw et al., 2001; Wen and Sachdev, 2004). WMH lesions are even more extensive in those with vascular or Alzheimer’s disease (AD) type of dementia when compared with cognitively normal older adults, suggesting its role in dementia pathogenesis and neurocognitive dysfunction (Bombois et al., 2007; Kloppenborg et al., 2014; Lee et al., 2016). WMH is also frequently observed in patients with multiple sclerosis (MS) (Loizou et al., 2015; Newton et al., 2017). Qualitative and quantitative WMH characterization has been used as a biomarker to assist cerebral small vessel disease diagnosis and to assess treatment effects (Wardlaw et al., 2013). The pathogenic mechanisms of WMH are not well understood, and have been attributed to brain hypoperfusion, white matter demyelization, or both (Greenberg, 2006; Wardlaw et al., 2013). Furthermore, periventricular and subcortical deep WMHs may have different pathogenic mechanisms (Schmidt et al., 2011; Poels et al., 2012; Tseng et al., 2013). The commonly used methods for WMH quantification are to measure its regional or total volume (i.e., the sum of WMH voxel size) within the whole brain based on image tissue segmentation algorithms (DeCarli et al., 2005; Wardlaw et al., 2013). This method, however, neglects entirely the typological or morphological features of WMH lesions which may have important clinical significance as demonstrated in recent studies in patients with MS (Loizou et al., 2015; Newton et al., 2017). In this regard, imaging processing using deep learning may reveal image patterns related to disease progression (Zeng et al., 2021; Li et al., 2022). However, this approach has limitations in that the unique typological or morphological features of WMH and its underlying neurobiological mechanisms cannot be characterized.

White matter hyperintensities shape is a basic morphological feature which can be derived from high-resolution T₂ FLAIR images after tissue segmentation. Shape feature extraction, recognition, and classification can be implemented either in the original or the transformed image space (Khotanzad and Hong, 1990; Mikolajczyk et al., 2003; Carmichael and Hebert, 2004; Tahmasbi et al., 2011). The analysis of 3D shape has been widely applied in the fields of image processing and pattern recognition, such as terrain matching (Rodriguez and Aggarwal, 1990), object retrieval (Körtgen et al., 2003; Novotni and Klein, 2004), anatomical structure analysis (Gerig et al., 2001; Styner et al., 2005; Terriberry et al., 2005; Luders et al., 2006; Gerardin et al., 2009; Shen et al., 2012; Wachinger et al., 2015; Makkinejad et al., 2019), and protein structural similarity retrieval (Mak et al., 2008; Sael et al., 2008). In general, the feature vectors of 3D shapes are first extracted. Then the similarity between the vectors is indexed for comparison, clustering, and recognition. The feature representation of 3D shape is to transform the original space of 3D objects to a high-dimensional feature vector space while preserving the shape information. The resulting feature vector (also known as a shape descriptor) can be used to characterize the unique shape of an object. In (Zhang et al., 2007), the computational techniques used for obtaining shape descriptors were comprehensively reviewed, and categorizations of the approaches were also provided. A shape descriptor, in general, needs to be able to characterize both the global shape contour and the regional topological details (Deng et al., 2016). Additionally, in order to assess the reliability and the accuracy of the descriptor, the descriptor must be able to reconstruct as close to the original object as possible. Due to the complexity of 3D shape feature extraction and the computational instability of numerical feature values, low-order 3D descriptors of objects, especially with the voxel-based approach, were commonly found in current literature (Novotni and Klein, 2004; Venkatraman et al., 2009a). Although a low-order shape descriptor may provide sufficient information for classification of objects at a coarse level, higher-order shape features are required to differentiate subtle regional topological differences in objects with fine structures. Therefore, choosing an appropriate order of shape descriptor is crucial to represent a 3D shape with different morphological features.

Based on literature, three categories of algorithms have applied to study the shapes of 3D objects: (1) surface-based methods using spherical harmonics as the basis functions (Kelemen et al., 1999; Styner et al., 2006), (2) voxel-based methods based on 3D Fourier (Vranic and Saupe, 2001) and Zernike transforms (Novotni and Klein, 2003), and (3) spectrum-based methods by solving the eigenvalues of 2D and 3D Laplace-Beltrami operators on triangular (boundary) and tetrahedral (volume) meshes (Reuter et al., 2006; Lian et al., 2013). All these methods theoretically can characterize 3D object shapes with high fidelity and have been applied to characterize brain structures (Joshi et al., 1997; Gerig et al., 2001; Styner et al., 2006; Wachinger et al., 2015, 2016). However, the surface-based methods are only suitable for smooth shapes with spherical topology, and cannot characterize 3D structures with holes or torus-like surfaces. Spectrum-based methods are isometry invariance, but cannot easily distinguish resembling objects (Reuter et al., 2006; Lian et al., 2013). Among the methods used for shape feature extraction, the voxel-based 3D Zernike transformation has its unique advantages. In contrast to surface-based methods (Zhang et al., 2007), 3D Zernike transform can characterize holes and torus. Compared to the surface-based method based on spherical harmonics, the 3D Zernike transformation combines spherical harmonics with radial polynomials to produce a more compact representation and requires fewer expansion orders (Venkatraman et al., 2009a). The advantages of 3D Zernike descriptor over spherical harmonics have been demonstrated in the benchmark studies involving image retrieval for general 3D objects (Novotni and Klein, 2004) and protein molecules with similar global structures (Sael et al., 2008). Furthermore, Zernike transform is rotational invariance in space, but Fourier transform is not which may lead to the dependence of the derived shape features on the orientation of original objects.

Based on the orthogonality of 2D Zernike moments (Teague, 1980), Canterakis generalized the classical 2D Zernike polynomials to 3D and derived 3D Zernike polynomials and moments. With this theory, Novotni and Klein developed methods for computing the Zernike moments and object reconstruction (Novotni and Klein, 2003, 2004). Similar to the 2D Zernike moments, the magnitudes of complex 3D Zernike moments, named the Zernike descriptor, is rotational invariant.

The 3D Zernike descriptors have been successfully applied to protein structural similarity retrieval (Mak et al., 2008; Sael et al., 2008), protein-protein docking using region-based (Venkatraman et al., 2009b), terrain matching (Ye and Chen, 2012), and amphetamine-type stimulant (ATS) drugs identification (Pratama et al., 2015). The maximum order of 3D Zernike moments used in most, if not all, studies is below 30 and can only represent the rough shape features of 3D objects because the existing computation method developed by Novotni and Klein is time-consuming and computationally instable for calculating higher orders of Zernike moments (Zhang et al., 2007). Recently, we have proposed a new algorithm to calculate high order 3D Zernike moments to characterize objects which have fine structures (Deng and Gwo, 2020).

Texture is also a basic feature of the surface appearance of an object, and one of the important morphological features of images. Textured surfaces are the core of human vision because they are important visual cues about surface characteristics. Texture information is used to identify objects and understand the pre-attentional vision of the scene (Depeursinge et al., 2014). Texture analysis has received widespread attention because of its important role in the field of computer vision and pattern recognition, including facial analysis, industrial inspection, satellite or aerial image analysis, biomedical image analysis, and biometrics object recognition (Depeursinge et al., 2017; Liu et al., 2017). Image texture can be assessed using several quantitative approaches such as structural, spectral transformation, modeling, or statistics-based approaches (Bharati et al., 2004). For image texture analysis, the statistical approaches have the benefits of being rotational, size, and translational invariant in the feature vector spaces. Furthermore, these methods can characterize image intensity distributions directly (Bharati et al., 2004; Castellano et al., 2004). In addition, they require fewer a priori model assumptions, such as basic symbolic image elements and repeated image patterns (Castellano et al., 2004). We used a statistical method based on fuzzy logic to construct the image intensity histogram of WMH lesions, and then cluster lesion features into several groups. Significant differences in the intensity of lesions can be observed in the resulting groups (Gwo et al., 2019). However, the concept of three-dimensional texture is rarely used because textures that exist in more than two dimensions cannot be fully visualized by humans. Computer graphics only provide virtual navigation in multi–planar rendering or translucent visualization and allow observation of two-dimensional projections of opaque textures (Toriwaki and Yoshida, 2009). Depeursinge et al. provided a good review of the challenges and opportunities of 3D texture analysis in biomedical imaging (Depeursinge et al., 2014).

White matter hyperintensities morphological characteristics such as the size, shape, and image texture may change with time which may reflect the progression of underlying dynamic pathophysiological process (Novotni and Klein, 2003; Lian et al., 2013). In this regard, recent studies have shown that the immediate surrounding areas of the defined WMH lesions may be at risk for further tissue damage and conversion to lesions (Reuter et al., 2006; Wachinger et al., 2016). These areas are classified as WMH penumbras (Wachinger et al., 2016). To characterize WMH lesions as well as their penumbras, we developed a seed-based region-growing algorithm to characterize 2D WMH boundaries to explore the potential growth of WMH lesions. We defined this specific WMH boundary characteristic as potential growth index (PGI) and observed that both shape and texture characteristics of 2D WMH are related to PGI (Gwo et al., 2019).

The characterization and quantification of the shape and texture of 3D WMH lesions have not been previously attempted. In this work, we extended 2D to high-order 3D Zernike transform to study the shape of 3D WMH and applied fuzzy logic to the intensity histogram of 3D WMH lesions for texture feature extraction. We also explored the potential growth index to predict of future tissue damage surrounding the 3D WMH lesion, that is, the WMH penumbra. Finally, we evaluated 3D potential growth differences among different lesion shape categories and texture clusters.

Methods and results

Magnetic resonance images acquisition

Full-brain 3D T₂ FLAIR images with voxel size of 1 × 1 × 1 mm³ were collected on a GE Discovery MR 750W 3T MRI scanner (GE Healthcare, Waukesha, WI) with the following parameters: sagittal, time of echo (TE) = 115 ms, time of repetition (TR) = 6.8 s, time of inversion (TI) = 1828 ms, echo train length = 200, receiver bandwidth = 41.67 kHz, fat saturation on, field of view (FOV) = 25.6 × 25.6 cm, slice thickness = 1 mm, number of slices = 176, acquisition matrix size = 256 × 256. All subjects signed informed consent approved by the Institutional Review Boards of the UT Southwestern Medical Center and Texas Health Presbyterian Hospital of Dallas. Thirty-two T₂ FLAIR brain image datasets (15 male, 17 female, 66.7 ± 6.0 years old and normal cognition), which contained clearly identifiable white matter hyperintensity (WMH) lesions with various sizes, were selected from an on-going HIPAC clinical trial (NCT03354143). Inclusion criteria: (1) age 55–79 years; (2) Mini-Mental State Exam (MMSE) > 26 to exclude dementia; (3) normotensive subjects and patiens with hypertension Exclusion criteria: (1) severe cerebrovascular disease such as stroke, transient ischemic attack, traumatic brain injury; (2) clinical diagnosis of dementia or other neurodegenerative diseases; (3) severe depression or other psychopathology; (4) unstable heart disease; (5) chronic kidney diseases with GFR < 40 ml/min; (6) orthostatic hypotension; (7) history of significant autoimmune disorders; (8) history of drug or alcohol abuse within the last 2 years; (9) uncontrolled diabetes mellitus; (10) obstructive sleep apnea; (11) regularly smoking cigarette within the past year; (12) severe obesity with BMI ≥ 45; (13) carotid stent or severe stenosis (> 50%); (14) pacemaker or other medical device of metal that precludes performing MRI; (15) history B12 deficiency or hypothyroidismT₂ FLAIR Image Segmentation.

T₂ FLAIR WMH regions were segmented on each 3D image volume through the lesion prediction algorithm (LPA) implemented in the Lesion Segmentation Toolbox (LST) version 2.0.12 for Statistical Parametric Mapping (SPM12). In LPA, the algorithm is trained using a logistic regression model on T₂ FLAIR brain images from 53 MS patients with severe lesion patterns. LPA was also validated in other patient populations such as older adults with diabetes (Styner et al., 2005). The fitness of a new T₂ FLAIR brain image to this model provides an estimate of lesion probability for each voxel in the image. In this study, we used a threshold of 0.5, as suggested by LST, on the obtained lesion probability maps to identify WMH regions. The segmentation accuracy was further verified through visual inspection. Figure 1 shows an example of the segmentation.

FIGURE 1

Figure 1. (A) An example a slice of T₂ FLAIR image volume from one subject showing multiple white matter hyperintensity (WMH) lesions; (B) the results of WMH segmentation using the lesion prediction algorithm (LPA) showing in red, and (C) a WMH binary mask after tissue segmentation, which was used to form the 3D structure of WMH for feature extraction.

Lesion size distribution

We extracted the WMH3D lesions greater than 30 mm³ in each subject and obtained a total number of 280 lesions. The lesion size distributions of all subjects are shown in Figure 2. The figure shows a wide range of lesion sizes and most of the lesions are relatively small. To explore whether the shape of WMH3D is related to its potential growth, the volume of all WMH3D is required to normalize to the same size. However, the volume scaling process can add or lose 3D shape details, and more so when the size distribution of WMH3D has a wide range as in our subjects (Figure 2). To reduce this scaling issue, this study divided the 280 WMH3D lesions to two groups based on their size. Group S_s had lesion size smaller than or equal to 250 voxels, and Group S_l had lesions larger than 250 voxels The group division generated 206 lesions for S_s and 74 lesions for S_l. To understand the anatomical distribution of the 280 WMH3D lesions, we classified the WMH clusters within or adjacent to the ventricle borders with a 3-mm thickness as periventricular WMH, and the rest as deep WMH. We found that 84.6% of the WMH clusters were at the periventricular region, and they tended to be relatively large with a volume size of 587.3 ± 1660.5 (mm³) with a range of 31–12409 mm³. The deep WMH clusters tended to be small with a volume size of 72.0 ± 44.4 (mm³), with a range of 31–256 mm³.

FIGURE 2

Figure 2. The histogram of WMH3D lesion size distribution for all subjects. The cumulative lesion size distribution is shown in percentage on the right vertical axis.

The rationale of using 250 voxels as the cut-off for the lesion groups was based on both volume and shape characteristics of the lesions. All extracted WMH3D lesions were positioned to 150 × 150 × 150 cubes with their centroids at the centers. Then 3D Zernike transform with an order of 150 was applied to the cubes to obtain the corresponding shape descriptors, and then a K-mean algorithm was used to cluster all 280 lesions into four clusters which are determined by the gap statistic (GAP, described in the later section) (Tibshirani et al., 2001). This process assisted in finding a size classification based on both the volume and shape of the lesion to reduce the influence of the shape changes caused by volume normalization on subsequent analyses. The size distribution results of GAP clustering are shown in Figure 3. The cluster shown in Figure 3A contains 187 lesions. The size of 250 voxels is a proper cutoff and thus was chosen to divide 280 lesions by their sizes into two groups.

FIGURE 3

Figure 3. The four clusters with the histogram of the WMH3D volume size distribution are determined by gap statistic. The numbers of WMH3D in each cluster (A–D) are 187, 59, 25, and 9, respectively.

WMH3D shape feature extraction and classification

WMH3D shape feature extraction using 3D Zernike transformation

The 2D Zernike transformation is based on the Zernike polynomials defined on the unit disc D. This transformation has been extensively applied to imaging shape feature extraction and pattern recognition (Papakostas et al., 2007; Wee and Paramesran, 2007; Gwo et al., 2019). The coefficients of the Zernike polynomial expansion of an object are called Zernike moments (ZMs). The magnitude of the ZMs, which is also named as the Zernike descriptor, is rotational invariant and represents the shape features of the analyzed objects. To define the 3D version of Zernike polynomials, the unit disc D is replaced by a unit ball B. Every point (x,y,z) in the unit ball B can be represented by a spherical coordinates (r, θ, ϕ) as shown in Eq. 1,

$$Z_{nm}=\frac{n+1}{\mathrm{\pi }}{\int }_0^{2\pi }{\int }_0^{1}f(r,\theta ) V_{nm}^{*} (r,\theta ) rdrd\theta ,\left|r\right|\le 1,(1)$$

where

$$\begin{array}{ccc}r\hfill & =\hfill & \sqrt{x^2+y^2+z^2}\hfill \\ \mathrm{\theta }\hfill & =\hfill & {\cos }^{-1}\frac{z}{r}\hfill \\ \varphi \hfill & =\hfill & {\sin }^{-1}\frac{y}{\sqrt{x^2+y^2}}\hfill \end{array}(2)$$

Canterakis introduced the first algorithm to calculate 3D Zernike moments (3DZMs) (Canterakis, 1999), where the 3DZMs were expressed as the linear combination of geometric moments. These 3DZMs were later described as shape descriptors for shape retrieval (Novotni and Klein, 2004). Canterakis’ algorithm has been applied to terrain matching (Wang et al., 2018a,b) and protein–protein interface prediction (Daberdaku and Ferrari, 2018). However, Canterakis’ algorithm could only be used to compute ZM up to the order of 25, due to computational demand and instability. Hosny et al. introduced a fast algorithm using eight ways of (anti-)symmetries (Hosny and Hafez, 2012). To overcome the limitations on computational efficiency and the maximum ZM order that can be computed reliably in previous algorithms, Deng and Gwo proposed a new algorithm based on a recursive approach to calculate 3D Zernike radial polynomials, as described in Eqs 10–13 (Deng and Gwo, 2020). The algorithm used to calculate the 3D Zernike polynomial is briefly described below.

The 3D Zernike polynomial $V_{n\mathrm{\ell }}^{\mathit{m }}(r,\theta ,\varphi )$ is defined as the multiplication of spherical harmonic $Y_{\mathrm{\ell }}^{\mathit{m }}(\theta ,\varphi )$ and radial polynomial R_nℓ (r)as below:

$$V_{n\mathrm{\ell }}^{\mathit{m }}(r,\theta ,\varphi )=Y_{\mathrm{\ell }}^{\mathit{m }}(\theta ,\varphi ){\mathit{ R}}_{n\mathrm{\ell }}\left(r\right)(3)$$

$Y_{\mathrm{\ell }}^{\mathit{m }}(\theta ,\varphi )$ and R_nℓ (r) are computed separately. The spherical harmonic $Y_{\mathrm{\ell }}^{\mathit{m }}(\theta ,\varphi )$ of degree ℓ with order m is given by

$$Y_{\mathrm{\ell }}^{\mathit{m }}(\theta ,\varphi )=\{\begin{array}{cc}\hfill K_{\mathrm{\ell }}^m{P}_{\mathrm{\ell }}^m(\cos \theta )\cos m\varphi \text{if}m\ge 0\hfill & \\ \hfill K_{\mathrm{\ell }}^m{P}_{\mathrm{\ell }}^m(\cos \theta )\sin m\varphi \mathit{\hspace{0.5em}}\text{otherwise}\hfill & \end{array}(4)$$

where $P_{\mathrm{\ell }}^{\mathit{m }}(\cdot )$is the associated Legendre polynomial of degree ℓ, given by

$$P_{\mathrm{\ell }}^{\mathit{m }}\left(x\right)={\left(-1\right)}^m{\left(1-x^2\right)}^{\frac{m}{2}}\frac{d^m}{d{x}^m}\left(P_{\mathrm{\ell }}\left(x\right)\right)(5)$$

and $K_{\ell }^m$ is the normalizing factor given by

$$K_{\mathrm{\ell }}^m={(-1)}^m\sqrt{\frac{{\in }_m\left(2\mathrm{\ell }+1\right)\left(\mathrm{\ell }-m\right)!}{\left(\mathrm{\ell }+m\right)!}}\text{where}{\in }_m=\{\begin{array}{cc}\hfill 1\hfill & \hfill \text{if}m\ne 0\hfill \\ \hfill 2\hfill & \hfill \text{oterwise}\hfill \end{array}(6)$$

Let $\widetilde{P_{\mathrm{\ell }}^{\mathit{m }}}\left(\cos \theta \right)=K_{\mathrm{\ell }}^m{P}_{\mathrm{\ell }}^{\mathit{m }}(\cos \theta )$ be the normalized associated Legendre polynomial. Then Eq. 4 is simplified to

$$Y_{\mathrm{\ell }}^{\mathit{m }}(\theta ,\varphi )=\{\begin{array}{cc}\hfill \widetilde{P_{\mathrm{\ell }}^{\mathit{m }}}\left(\cos \theta \right)\mathrm{ cos}m\varphi \text{if}m\ge 0\hfill & \\ \hfill \widetilde{P_{\mathrm{\ell }}^{\mathit{m }}}\left(\cos \theta \right)\mathrm{ sin}m\varphi \text{otherwise}\hfill & \end{array}(7)$$

The spherical harmonics $Y_{\mathrm{\ell }}^m(\cdot )$form an orthonormal basis for the Hilbert space L²(S²) of the square-integrable functions over the unit sphere S². For any function f of L²(S²), f can be expressed as in Eq. 8 (Szegõ, 1939):

$$f(\theta ,\varphi )=\sum _{\mathrm{\ell }=0}^{\mathrm{\infty }}\sum _{\text{m}=-\mathrm{\ell }}^{\mathrm{\ell }}{C}_{\mathrm{\ell }}^m{Y}_{\mathrm{\ell }}^m(\theta ,\varphi )(8)$$

where $C_{\ell }^m$ are the coefficients; ℓ is a non-negative integer; m is an integer with |m| ℓ. The computation procedures of $\widetilde{P_{\mathrm{\ell }}^{\mathit{m }}}\left(\cos \theta \right)$ for degree ℓ ≤ ℓ_max are summarized to the following (Szegõ, 1939; Deng and Gwo, 2018a):

1. Initialize $\widetilde{P_0^{\mathrm{0 }}}\left(\cos \theta \right)=\sqrt{\frac{1}{4\mathrm{\pi }}}$, which is the normalizing factor for volumetric integration. Then iteratively calculate the following:

2. $\widetilde{P_{\mathrm{\ell }}^{\mathrm{\ell }}}\left(\cos \theta \right)={\text{C}}_3\sin \theta \widetilde{P_{\mathrm{\ell }-1}^{\mathrm{\ell }-\mathrm{1 }}}\left(\cos \theta \right)$ for ℓ1,2,3,…,ℓ_max

3. $\widetilde{P_{\mathrm{\ell }}^{\mathrm{\ell }-\mathrm{1 }}}\left(\cos \theta \right)={\text{C}}_1\cos \theta \widetilde{P_{\mathrm{\ell }-1}^{\mathrm{\ell }-\mathrm{1 }}}\left(\cos \theta \right)$

4. $\widetilde{P_{\mathrm{\ell }}^{\mathit{m }}}\left(\cos \theta \right)={\text{C}}_1\cos \theta \widetilde{P_{\mathrm{\ell }-1}^{\mathit{m }}}\left(\cos \theta \right)-{\text{C}}_2\widetilde{P_{\mathrm{\ell }-\mathrm{2 }}^m}\left(\cos \theta \right)$

$$\begin{array}{ccc}\mathit{f}orm\hfill & =\hfill & 0,\mathrm{\hspace{0.25em}1},..,\mathrm{\ell }-2\hfill \\ \hfill \text{C}& =\hfill & \sqrt{\frac{2\mathrm{\ell }+1}{(\mathrm{\ell }+\mathrm{m})(\mathrm{\ell }-\mathrm{m})}},{\text{C}}_1=\text{C}\sqrt{2\mathrm{\ell }-1}\hfill \\ \hfill {\text{C}}_2& =\hfill & \text{C}\sqrt{\frac{(\mathrm{\ell }+\mathrm{m}-1)(\mathrm{\ell }-\mathrm{m}-1)}{2\mathrm{\ell }-3}},{\text{C}}_3=\sqrt{\frac{2\mathrm{\ell }+1}{2\mathrm{\ell }}}\hfill \end{array}(9)$$

For a Zernike polynomial order n (a non-negative integer), the integer ℓ above needs to ≤ n and n − ℓ = even, and the integer m above needs to |m| ≤ ℓ.

The 3D Zernike radial polynomial R_nℓ(r) in Eq. 3 was originally given in terms of Jacobi polynomials as described in (Szegõ, 1939), but different calculation methods of 3D Zernike radial polynomial have been proposed (Deng and Gwo, 2020). In our work, the R_nℓ is computed recursively, similar to the Kintner’s P-method in the case of 2D Zernike polynomials (Kintner, 1976; Deng and Gwo, 2018b), and is presented in Eq. 10.

$$\begin{array}{c}{R}_{n\mathrm{\ell }}\left(r\right)=\left(K_1{r}^2+K_2\right){\mathit{ R}}_{n-2,\mathrm{\ell }}\left(r\right)+K_3{R}_{n-4,\mathrm{\ell }}\left(r\right)\hfill \\ \hfill \mathit{f}orn=\mathrm{\ell }+4,\mathrm{\ell }+6,\mathrm{\dots },n_{max}\hfill \end{array}(10)$$

where the coefficients K_i are given by the following,

$$\begin{array}{c}\hfill k_0=\left(n-\mathrm{\ell }\right)\left(n+\mathrm{\ell }+1\right)\left(2n-3\right)\hfill \\ \hfill k_1=\left(2n-1\right)\left(2n+1\right)\left(2n-3\right)\hfill \\ \hfill k_2=\left(-2n+1\right)\left(\frac{4{\mathrm{\ell }}^2+4\mathrm{\ell }+1}{2}\right)-\frac{k_1}{2}\hfill \\ \hfill k_3=-\left(n-\mathrm{\ell }-2\right)\left(n+\mathrm{\ell }+1\right)\left(2n+1\right)\hfill \\ \hfill K_1=\frac{k_1}{k_0},K_2=\frac{k_2}{k_0},K_3=\frac{k_3}{k_0}\hfill \end{array}(11)$$

For this recursive formula, the following initial equalities are also required:

$$R_{nn}\left(r\right)=r^{n}forn0,\mathrm{\hspace{0.25em}1},\mathrm{\hspace{0.25em}2},(12)$$

and

$$R_{n,n-2}\left(r\right)=(n+\frac{1}{2})r^n-(n-\frac{1}{2})r^{n-2}forn=2,\mathrm{\hspace{0.25em}3},\mathrm{\hspace{0.25em}4},\mathrm{\dots }(13)$$

Let f (r, θ, ϕ) be a 3D image function within the unit ball B. The 3DZM ${\text{Z}}_{n\mathrm{\ell }}^m$ can be regarded as the inner product of the image function f (r, θ, ϕ) with the basis function $V_{n\mathrm{\ell }}^{\mathit{m }}(r,\theta ,\varphi )$ (Deng and Gwo, 2020), and can be described as

$$Z_{n\ell }^m=(2n+3)\underset{(r,\theta ,\varphi )\in B}{{\displaystyle \iiint }}f(r,\theta ,\varphi )V_{n\ell }^m(r,\theta ,\varphi )r^2\sin \theta drd\theta d\varphi (14)$$

Each moment within Order n corresponds to a (2ℓ+1)-dimensional vector ${\overrightarrow{Z}}_{n\ell }$ as

$$\stackrel{\rightharpoonup }{\mathit{Z}}=(Z_{n\mathrm{\ell }}^{-\mathrm{\ell }},Z_{n\mathrm{\ell }}^{-\mathrm{\ell }+1},\mathrm{\cdots },Z_{n\mathrm{\ell }}^0,\mathrm{\cdots },Z_{n\mathrm{\ell }}^{\mathrm{\ell }-1},Z_{n\mathrm{\ell }}^{\mathrm{\ell }})(15)$$

The l²-norm of ${\overrightarrow{Z}}_{n\ell }$, denoted by

$$||\stackrel{\rightharpoonup }{\mathit{Z}}||=\sqrt{\sum _{m=-\mathrm{\ell }}^{\mathrm{\ell }}{\left|Z_{n\mathrm{\ell }}^m\right|}^2}(16)$$

is rotation invariant, and can be used as the 3D shape descriptor (or Zernike descriptor) of a 3D object. The total number of 3DZMs and the dimension of Zernike descriptor for an expansion up to order n are given by Eqs 17, 18, respectively:

$$\mathrm{Number}\mathrm{of}\mathrm{3DZMs}=\sum _{\mathrm{i}0}^{\mathrm{n}}\lfloor \frac{{\left(\mathrm{i}+2\right)}^2}{4}\rfloor (17)$$

$$\mathrm{Dimension}\mathrm{of}\mathrm{Zernike}\mathrm{descriptor}=\{\begin{array}{cc}\hfill {\left(\frac{n+2}{2}\right)}^2\hspace{1em}\hspace{1em}\hspace{0.5em}\hspace{1em}\mathrm{if}\mathrm{order}n\mathrm{is}\mathrm{even}\hfill & \\ \hfill \frac{\left(n+3\right)\left(n+1\right)}{4}\mathrm{if}\mathrm{order}n\mathrm{is}\mathrm{odd}\hfill & \end{array}(18)$$

The image object function f can be reconstructed with ZM order M as f_M below:

$$f_M(r,\theta ,\varphi )=\sum _{n0}^M\sum _{\mathrm{\ell }}\sum _{m-\mathrm{\ell }}^{\mathrm{\ell }}{Z}_{n\mathrm{\ell }}^m{V}_{n\mathrm{\ell }}^m(r,\theta ,\varphi )(19)$$

When M is large enough, the function f_M can be used to approximate the original image f (Deng and Gwo, 2018b). For a binary shape with the background represented by 0, the error rate ℰ_r between the original image f and the reconstructed f_M can be calculated by

$$\begin{array}{c}{ℰ}_r=\frac{{\sum }_{\forall (x,y,z)}\left(f(x,y,z)\oplus f_M(x,y,z)\right)}{{\sum }_{\forall (x,y,z)}f(x,y,z)}\hfill \\ \text{ where}{f}_M(x,y,z)=\{\begin{array}{cc}\hfill 1\text{if}{f}_M(x,y,z)\ge 0.5\hfill & \\ \hfill \mathrm{0 }\text{otherwise }\hfill & \end{array}\hfill \end{array}(20)$$

where ⊕ is exclusive disjunction and f (x,y,z) 0 or 1. Based on the error rate ℰ_r, an appropriate ZM order M can be chosen.

Overall, the calculation of 3D Zernike moments is summarized as follows: First, the normalized associated polynomial $\widetilde{P_{\mathrm{\ell }}^m}\left(\cos \theta \right)$ of the spherical harmonic function is calculated using Eq. 9. Second, the 3D Zernike radial polynomial R_nℓ is recursively calculated using Eqs 10–13. Then, the 3D Zernike polynomials can be obtained by Eq. 3. Finally, Eq. 14 is used to generate 3D Zernike moments. The 3D spherical harmonics and Radial polynomials are illustrated in Supplementary Figure 1.

For 3D Zernike Transformation of WMH, all lesions were linearly rescaled at the ratio of $\sqrt[3]{v\prime /v}$ in three dimensions, and the intensity of the resulting voxel was calculated by tricubic interpolation (Arata, 1995), where v was the original volume size and the v′ was the volume scaled to. Due to blurring effect in scaling, an appropriate intensity threshold was then chosen so that the scaled volume is closest to v′. To compare the Zernike descriptors at the same scale, the lesions in Group S_s are normalized to about 80 voxels within a 36 × 36 × 36 cube, with the centroid at the center of the cube. Similarly, the lesions in Group S_l are normalized to about 1500 voxels within a 90 × 90 × 90 cube. These two numbers, 80 and 1500, are medium size in their own lesion groups.

As an example, the effect of ZM order on the Zernike transformation of two WMH3D lesions and their reconstruction accuracy is shown in Supplementary Figure 2.

As shown in Supplementary Figure 2 qualitatively and in Figure 4 quantitatively, the order of the Zernike transformation needs to be large enough to preserve the original shape details.

FIGURE 4

Figure 4. The error rates of two lesions calculated from two different normalized sizes at different Zernike reconstruction orders: (A) the error rate of an 80-voxel lesion in a 36 × 36 × 36 voxel cube approached zero at order around 100, and (B) the error rate of a 1500-voxel lesion in a 90 × 90 × 90 voxel cube approached zero at order around 250.

As illustrated in Figure 4, the error rates of WMH3D lesion reconstruction, ℰ_r, would decrease with the reconstruction order. The 3D Zernike transformation with order of 100 for Group S_s and then 250 for Group S_l, were applied to the voxel cubes containing the size-normalized lesions. The error rates of group S_s and S_l were 7.3 × 10^–4± 3.4 × 10^–3 and 7.2 × 10^–4± 1.3 × 10^–3, respectively. Zernike descriptor was obtained for each lesion with 2,601 dimensions for Group S_s and 15,876 for Group S_l, based on Eqs 16, 18. The principal component analysis (PCA) was then used to reduce the large number of dimensions of the Zernike descriptors to reduce computation demand while minimizing information loss. To maintain 99.8% variance of the two lesion size groups, the number of principal component choices for the two groups are 80 and 68, respectively (Figure 5). The Zernike descriptors were projected on the selected principal components (eigenvectors), and the coefficients corresponding to these principal components were used for WMH3D shape clustering and classification.

FIGURE 5

Figure 5. To maintain 99.8% variance of the two lesion size groups, (A) the first 80 principal components and (B) the first 68 principal components were chosen for group S_s and group S_l, respectively.

WMH3D shape classification

A K-means algorithm was used for clustering and classification due to its simplicity and efficiency (Bharati et al., 2004). However, different initial seeds used in the clustering algorithm may generate different clustering results (Hamerly and Elkan, 2002). In this study, we randomly selected the initial clustering seeds from the shape lesion feature space and conducted 1,000 trials to assess the clustering results. We employed a gap statistic method to determine the optimal number of clusters for WMH3D shape clustering (Tibshirani et al., 2001).

1. The sum of the within-cluster dispersion W_k is computed for each choice k clusters (k = 1, 2,…, N).

$$W_k=\sum _{r=1}^k\sum _{x_i\in C_r}{\left(x_i-{\overline{x}}_r\right)}^2(21)$$

where x_i is a data point, C_r denotes cluster r, and ${\overline{x}}_r$ is the vector mean of C_r.

The B reference datasets is uniformly generated by randomly sampling from the bounding rectangle of the original dataset. By Eq. 21, $W_{kb}^{*}$ is computed for each k and b (b = 1, 2, …, B, k = 1, 2, …, N). Then, compute the estimated gap statistic

$$Gap\left(k\right)=\frac{1}{B}\sum _b\text{log}\left(W_{kb}^{*}\right)-\text{log}\left(W_k\right)(22)$$

2. let l = (1/B)∑_b log (W_kb), compute the standard deviation

$$s{d}_k={\left[\frac{1}{B}\sum _b{\left(\text{log}\left(W_{kb}^{*}-l\right)\right)}^2\right]}^{1/2}(23)$$

Let $s_k=s{d}_k\sqrt{\left(1+1/B\right)}$. Choose the optimal number of shape clusters via

$$\widehat{k}=\text{smallest}k\mathrm{such}\mathrm{that}Gap\left(k\right)\ge Gap\left(k+1\right)-s_{k+1}(24)$$

In gap statistic procedure, N is a pre-selected number of clusters such that $\widehat{k}$ can be determined in the range of [1, N], and B is selected to calculate the value of sd_k in a statistical sense. In this study, N and B were set to 20 and 30, respectively.

For the size-normalized lesions, the feature dimension were 80 in Group S_s and 68 in Group S_l, the Gap values were calculated and displayed in Supplementary Figure 3; the optimal numbers of shape clusters were 5 for Group S_s and 4 for Group S_l according to Eq. 24.

Figure 6 shows the WMH3D shape clustering results using the K-means clustering algorithm (Hartigan and Wong, 1979) based on the cluster number of 5 for Group S_s and 4 for Group S_l. The second column shows the number of lesions in each shape cluster. The third column presents the four lesion images closest to their cluster means, as the representative lesion shapes corresponding to their clusters. All lesion images shown in the figure were normalized close to the voxel size of 80 for Group S_l and 1500 for Group S_l. The orientation-adjusted images can be seen in Figure 6 with the significant differences among the shape clusters.

FIGURE 6

Figure 6. A K-means algorithm was performed to categorize WMH3D lesion shapes to different clusters for Group S_s (A) and Group S_l (B).

WMH3D texture feature and classification

Texture feature extraction

To obtain rotation-invariant texture features that are applicable to both small and large-size WMH3D, the fuzzy logic technique that we developed in (Gwo et al., 2019) was extended to WMH3D. Specifically, when segmenting WMH3D, false positives likely occurred at boundaries of a lesion, where signal intensity was usually low and thus led to biased estimation. To reduce chance of false positive, voxels with intensity at the lowest 1% were discarded. To reduce the variations of the signal intensity of individual subjects, a min-max normalization was applied to a WMH3D to normalize its voxel intensity based on Eq. 25.

$$s(x,y,z)=\frac{f(x,y,z)-gMin}{gMax-gMin}\mathrm{for}f(x,y,z)=\text{max}\left((f(x,y,z),gMin)\right)(25)$$

where f (x,y,z) is the intensity of voxel f (x,y,z) and s ∈ [0,1], gMax = maximal voxel intensity of WMH3D and gMin = minimal voxel intensity of WMH3D.

For feature extraction, each voxel intensity was quantized into one of the n bins to create a histogram that represents voxel intensity distribution of a WMH3D. To minimize the interference of quantization to the frequency histogram, we used a fuzzy logic method (Gwo and Wei, 2013) to allocate normalized voxel intensity values to each of the pre-selected bins. Specifically, a normalized voxel intensity s was assigned proportionally two values, called fuzzy values, to the two neighboring bins according its relative positions to the bin centers (Supplementary Figure 4).

The fuzzy logic functions used for assigning voxels to the frequency histogram are presented in Eq. 26 (Gwo and Wei, 2013). The fuzzy value v[j] at bin j is calculated as:

$$\{\begin{array}{c}\hfill \begin{array}{c}v\left[0\right]=1\text{if}s\le \frac{1}{2n}\hfill \\ \begin{array}{c}v\left[j-1\right]=\frac{2j+1}{2}-s\times n\hfill \\ v\left[j\right]=s\times n-\frac{2j-1}{2}\hfill \end{array}\}\text{if}\frac{1}{2n}\le s\le \frac{2j+1}{2n}\hfill \\ \begin{array}{c}v\left[j\right]=\frac{2j+3}{2}-s\times n\hfill \\ v\left[j+1\right]=s\times n-\frac{2j+1}{2}\hfill \end{array}\}\text{if}1-\frac{1}{2n}\ge s>\frac{2j+1}{2n}\hfill \\ v\left[n-1\right]=1\text{if}s\ge 1-\frac{1}{2n}\hfill \end{array}\hfill \end{array}(26)$$

where n = the total number of bins, and j = 0,…, n-1. Since the sizes of WMH3D lesions vary in a wide range (Figure 2), the WMH3D intensity frequency distribution histograms need to be further normalized before they can be compared. Herein, each histogram is normalized to have a total accumulative frequency of 1.

WMH3D texture feature classification

Texture feature classification of individual WMH3D lesion images was conducted using a feature vector clustering method similar to those discussed above in the section of “WMH3D Shape Classification.” Of note, the texture feature vector was based on the frequency histogram presented above using the fuzzy logic method. The influences of different texture feature dimensions (i.e., the number of bins used to construct the intensity histogram) and the numbers of clusters on texture feature classification were explored using the same strategy discussed above for WMH3D shape feature clustering. The sum of within-cluster dispersion W_k value was calculated with the cluster number from 2 to 20 and the feature dimensions from 2 to 15. As illustrated in Supplementary Figure 5, W_k decreased with the increase of the cluster numbers. A noticeable “elbow” phenomenon was seen for a wide range of texture feature dimensions from 2 to 15. In (Hughes, 1968), there are two considerations in choosing an appropriate number of bins: (1) If the number of bins is too large, the fuzzy values accumulated in some bins become sparse, especially for small-sized lesions. Sparsity is problematic for any statistical analysis method. (2) Conversely, if the number of bins is too small, lesion features may not be effectively distinguishable. With these two reasons in mind, in this study, we selected ten bins for texture feature extraction.

The gap statistics discussed above was applied to determine the optimal number of texture feature cluster for pattern recognition based on the K-means algorithm for grouping (Hartigan and Wong, 1979). Supplementary Figure 6 shows that 5 is the optimal number of cluster.

WMH3D potential growth index

Voxel intensity is an important feature in image analysis. In this study, the intensity information of voxels in the penumbra area was used to estimate the likelihood of lesion development. The intuitive assumption is that the higher the intensity value, the higher the probability that the lesion may develop. If the intensity of neighboring voxels (penumbra) around the identified lesions is within a reasonable range discussed below, these will be the voxels of interest for potential growth. The distance between the voxels and the boundary of the lesion should also be considered. It hypothesized that the farther the voxel is from the lesion boundary, the greater the contribution to the PGI.

For each subject, the interesting voxel set P_w in lesion penumbra involves the calculation of PGI, and the corresponding intensity range is defined as follows:

$$P_w=\left\{p(x,y,z)\right|m-f(x,y,z)\le \mathrm{\gamma }\times \mathrm{\sigma }(27)$$

where f (x, y, z) denotes the voxel p(x, y, z) intensity, m is the average intensity of all WMH3Ds in the subject, σ is the corresponding standard deviation, and γ is a user-defined positive real number. Dilation morphological operation is applied to mask image to iteratively generate l layers masks with one-voxel thickness surrounding the lesion, which is the interest area of the penumbra to estimate the PGI of the lesion. The schematic analysis pipeline for calculating WMH3D PGI is shown in Figure 7.

FIGURE 7

Figure 7. Schematic analysis pipeline for calculating WMH3D potential growth index (PGI). (A) The 3D image volume of a subject, (B) the results of WMH segmentation using the lesion prediction algorithm, (C) two enlarged WMH3Ds whose locations are marked with red circles in panels (B,D) the intensities of potential penumbra voxels within the range of interest are displayed in a lighter color with two different viewing directions. The sizes of these two 3D lesions are 91 and 12409 voxels, respectively.

In this study, we chose γ = 2.5, which covers 99.38% of all WMH3D voxel intensities in a subject to demonstrate the presence of potential growth regions of WMH3D lesions. In the l-layer mask of the lesion, we are only interested in voxels with an intensity value greater than the value displayed by the red dotted line in the figure, which is m−2.5 × σ. These apparent layer masks are used to identify the relative location of a growth voxel. A growth voxel at an outer layers of these masks weights more in its contribution to the PGI. Specifically, the weight w_i at ith layer, with total l layers, is given by the following equation:

$$w_i=\frac{i}{{\sum }_{j=1}^{l}j}(28)$$

Once the number of growth voxels at each layers were evaluated, PGI for each WMH3D lesion is calculated below:

$$PGI=\frac{{\sum }_{i=1}^{l}G{V}_i{w}_i}{V_l}(29)$$

where, GV_i = number of “growth voxels” found at the ith layer, and V_l = the total number of voxels in all l layers for a WMH3D. All lesion images were evaluated for their PGIs with l set to 5. The PGIs estimated from the small lesions and the large lesions that are near the ventricle shown in Figure 7 are 0.0569 and 0.0844, respectively.

The relationships between potential growth index and WMH3D shape and texture features

For the shape and the texture clusters classified as shown in Figures 6, 8 above, one-way Analyses of Variance (ANOVAs) were performed to evaluate if there were significant differences in PGI among the shape or the texture clusters. Significant differences were found among both shape (P = 1.57 × 10^–3 for Group S_s and P = 3.14 × 10^–2 for Group S_l shown in Tables 1, 2) and texture (P = 1.79 × 10^–6 shown in Table 3) clusters. Tables 4–6 also show the results after the Bonferroni corrections for multiple comparisons. As a reference, lesion volume size analyses are also included. It is worth noting that among the shape clusters in the S_s group (Table 4 and Figure 6A), cluster 5 is significantly different from the other four clusters in terms of both PGI and lesion size. However, the PGI difference in cluster 5 from other four clusters was not likely driven by the lesion size because this cluster contains size evenly distributed between 30 and 250 voxels. For the S_l group, significant PGI differences were only found between Clusters 1 and 4 and between Clusters 2 and 4 (Table 5 and Figure 6B). However, due to the large lesion size variance in each cluster, there was no significant difference in lesion size among the shape clusters (Table 2). Furthermore, significant PGI differences were found between cluster 4 and the other three clusters (Table 6). Of note, compared with the other three clusters, the average PGI value of cluster 4 is smaller and the texture color is lighter (e.g., high intensity) (Figures 8, 9).

FIGURE 8

Figure 8. The WMH3D lesion images from the 32 subjects were classified into four clusters. The images at three view directions are displayed.

TABLE 1

Table 1. ANOVA analysis of the potential growth indices (PGIs) for shape clusters and corresponding lesion sizes in the S_s group.

TABLE 2

Table 2. ANOVA analysis of the potential growth indices (PGIs) for shape clusters and corresponding lesion sizes in the S_l group.

TABLE 3

Table 3. ANOVA analysis of the potential growth indices (PGIs) for texture clusters and corresponding lesion sizes for all 280 lesions.

TABLE 4

Table 4. PGI and lesion size differences between shape clusters in the S_s group.

TABLE 5

Table 5. Compare the PGI and lesion size between shape clusters in the S_l group.

TABLE 6

Table 6. Compare the PGI between texture clusters for all 280 lesions.

FIGURE 9

Figure 9. The 280 WMH3D lesion images from the 32 subjects were classified into four clusters based on their fuzzy image textures. The label on the y-axis is the probabilities of voxels in the lesion assigned to the bin.

Discussion and conclusion

In this study, we extended our prior work in WMH morphological analyses from 2D to 3D. A total of 280 3D lesions from 32 cognitively normal older adults with a volume size of greater than 30 voxels were used in shape and texture the analysis. When using Zernike transformation to extract shape features of 3D objects, volume normalization is a necessary process. In this regard, excessive scaling may enlarge the boundary details of small lesions while losing the details of large lesions. In this work, before clustering the shape of the lesions, WMH3D lesions were divided into two groups according to the size of the lesions to minimize volume normalization error. The texture features of the lesions used in this study were generated based on the intensity distribution. The fuzzy processing based on image intensity normalization for feature extraction reduced the influences of the intensity quantification. In addition, the intensity distribution was normalized by the size of the lesion, resulting in size independency.

The statistical data analyses showed that regardless of the volume size category of the lesions, PGI had significant differences among the shape clusters. The significant differences were also presented among the texture clusters. T₂ FLAIR WMH lesions were mostly located around the ventricles. The lesions around the ventricle were usually longer in shape and had high voxel intensity in texture, and had lower average PGI values than lesions distant from the ventricles. These observations together suggest that WMH lesion anatomic locations, morphological characteristics, as well as the lesion texture may have impact on lesion progression. Further work with a large sample size and a longitudinal study design would allow us to address these clinically significant questions.

When the number of lesions is sufficient with a large sample size of subjects, the merging of different lesion size groups performed in the present study to reduce the influence of lesion size on the application of Zernike transformation would not be necessary. The etiology of WNH is complex and can be multifactorial (Alber et al., 2019). Given that healthy subjects and patients with hypertension were enrolled this study, we suspect, but cannot prove that WHM lesions observed likely reflect the presence of cerebral small vessel disease (Alber et al., 2019). Whether WMH3D shape and texture characteristics and location are related to different etiology also worth further studies. Finally, studies are also needed to optimize the algorithms and parameters of shape and texture feature extraction, clustering, and PGI estimation with a goal to apply this novel imaging processing method to clinical research.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics statement

The studies involving human participants were reviewed and approved by the Institutional Review Boards of the UT Southwestern Medical Center and Texas Health Presbyterian Hospital of Dallas. The patients/participants provided their written informed consent to participate in this study.

Author contributions

C-YG conducted experiments, wrote code to analyze the data, interpreted the data, and wrote the manuscript. DZ prepared brain images and lesion segmentation. DZ and RZ interpreted the data, participated in the scientific discussions, and provided critical insights. All authors reviewed the manuscript and approved it for publication.

Funding

This study was supported in part by the National Institute of Health (NIH, R01AG057571).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2022.1028929/full#supplementary-material

References

Alber, J., Alladi, S., Bae, H.-J., Barton, D. A., Beckett, L. A., Bell, J. M., et al. (2019). White matter hyperintensities in vascular contributions to cognitive impairment and dementia (VCID): Knowledge gaps and opportunities. Alzheimers Dement. 5, 107–117. doi: 10.1016/j.trci.2019.02.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Arata, L. K. (1995). “Tricubic interpolation,” Graphics gems V, ed. A. W. Paeth (Cambridge, MA: Academic Press), 107–110.

Google Scholar

Bharati, M. H., Liu, J. J., and MacGregor, J. F. (2004). Image texture analysis: Methods and comparisons. Chemometr. Intell. Lab. Syst. 72, 57–71.

Google Scholar

Bombois, S., Debette, S., Delbeuck, X., Bruandet, A., Lepoittevin, S., Delmaire, C., et al. (2007). Prevalence of subcortical vascular lesions and association with executive function in mild cognitive impairment subtypes. Stroke 38, 2595–2597. doi: 10.1161/STROKEAHA.107.486407

PubMed Abstract | CrossRef Full Text | Google Scholar

Canterakis, N. (1999). “3D Zernike moments and Zernike affine invariants for 3D image analysis and recognition,” in Proceedings of the 11th Scandinavian conference on image analysis, Denmark.

Google Scholar

Carmichael, O., and Hebert, M. (2004). Shape-based recognition of wiry objects. IEEE Trans. Pattern Anal. Mach. Intell. 26, 1537–1552. doi: 10.1109/TPAMI.2004.128

PubMed Abstract | CrossRef Full Text | Google Scholar

Castellano, G., Bonilha, L., Li, L. M., and Cendes, F. (2004). Texture analysis of medical images. Clin. Radiol. 59, 1061–1069.

Google Scholar

Daberdaku, S., and Ferrari, C. (2018). Exploring the potential of 3D Zernike descriptors and SVM for protein–protein interface prediction. BMC Bioinformatics 19:35. doi: 10.1186/s12859-018-2043-3

PubMed Abstract | CrossRef Full Text | Google Scholar

De Leeuw, F., de Groot, J. C., Achten, E., Oudkerk, M., Ramos, L., Heijboer, R., et al. (2001). Prevalence of cerebral white matter lesions in elderly people: A population based magnetic resonance imaging study. The rotterdam scan study. J. Neurol. Neurosurg. Psychiatry 70, 9–14.

Google Scholar

DeCarli, C., Fletcher, E., Ramey, V., Harvey, D., and Jagust, W. J. (2005). Anatomical mapping of white matter hyperintensities (WMH). Stroke 36, 50–55.

Google Scholar

Deng, A.-W., and Gwo, C.-Y. (2018a). High degree coefficients of spherical harmonics expansions and image reconstruction. Kyoto: ACEAIT.

Google Scholar

Deng, A.-W., and Gwo, C.-Y. (2018b). Efficient computations for generalized Zernike moments and image recovery. Appl. Math. Comput. 339, 308–322.

Google Scholar

Deng, A.-W., and Gwo, C.-Y. (2020). “A stable algorithm computing 3D Zernike moments and shape reconstructions,” in Proceedings of the 2020 4th international conference on digital signal processing, New York, NY.

Google Scholar

Deng, A.-W., Wei, C.-H., and Gwo, C.-Y. (2016). Stable and fast computation of high-order Zernike moments using a recursive method. Pattern Recogn. 56, 16–25.

Google Scholar

Depeursinge, A., Al-Kadi, O. S., and Mitchell, J. R. (2017). Biomedical texture analysis: Fundamentals, tools and challenges. Cambridge, MA: Academic Press.

Google Scholar

Depeursinge, A., Foncubierta-Rodriguez, A., Van De Ville, D., and Müller, H. (2014). Three-dimensional solid texture analysis in biomedical imaging: Review and opportunities. Med. Image Anal. 18, 176–196. doi: 10.1016/j.media.2013.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerardin, E., Chételat, G., Chupin, M., Cuingnet, R., Desgranges, B., Kim, H.-S., et al. (2009). Multidimensional classification of hippocampal shape features discriminates Alzheimer’s disease and mild cognitive impairment from normal aging. Neuroimage 47, 1476–1486. doi: 10.1016/j.neuroimage.2009.05.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerig, G., Styner, M., Jones, D., Weinberger, D., and Lieberman, J. (2001). “Shape analysis of brain ventricles using spharm,” in Proceedings IEEE workshop on mathematical methods in biomedical image analysis (MMBIA 2001), (Piscataway, NJ: IEEE), 171–178.

Google Scholar

Greenberg, S. M. (2006). Small vessels, big problems. N. Engl. J. Med. 354, 1451–1453.

Google Scholar

Gwo, C.-Y., and Wei, C.-H. (2013). Plant identification through images: Using feature extraction of key points on leaf contours. Appl. Plant Sci. 1:apps.1200005.

Google Scholar

Gwo, C.-Y., Zhu, D. C., and Zhang, R. (2019). Brain white matter hyperintensity lesion characterization in T2 fluid-attenuated inversion recovery magnetic resonance images: Shape, texture, and potential growth. Front. Neurosci. 13:353. doi: 10.3389/fnins.2019.00353

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamerly, G., and Elkan, C. (2002). “Alternatives to the k-means algorithm that find better clusterings,” in Proceedings of the 11th international conference on information and knowledge management, New York, NY, 600–607.

Google Scholar

Hartigan, J. A., and Wong, M. A. (1979). Algorithm AS 136: A k-means clustering algorithm. J. R. Stat. Soc. Ser. C Appl. Stat. 28, 100–108.

Google Scholar

Hosny, K. M., and Hafez, M. A. (2012). An algorithm for fast computation of 3D Zernike moments for volumetric images. Math. Probl. Eng. 2012:353406. doi: 10.1109/TPAMI.2010.139

PubMed Abstract | CrossRef Full Text | Google Scholar

Hughes, G. (1968). On the mean accuracy of statistical pattern recognizers. IEEE Trans. Inf. Theory 14, 55–63.

Google Scholar

Joshi, S. C., Miller, M. I., and Grenander, U. (1997). On the geometry and shape of brain sub-manifolds. Int. J. Pattern Recogn. Artif. Intell. 11, 1317–1343.

Google Scholar

Kelemen, A., Székely, G., and Gerig, G. (1999). Elastic model-based segmentation of 3-D neuroradiological data sets. IEEE Trans. Med. Imaging 18, 828–839. doi: 10.1109/42.811260

PubMed Abstract | CrossRef Full Text | Google Scholar

Khotanzad, A., and Hong, Y. H. (1990). Invariant image recognition by Zernike moments. IEEE Trans. Pattern Anal. Mach. Intell. 12, 489–497.

Google Scholar

Kintner, E. C. (1976). On the mathematical properties of the Zernike polynomials. Opt. Acta 8, 679–680.

Google Scholar

Kloppenborg, R. P., Nederkoorn, P. J., Geerlings, M. I., and van den Berg, E. (2014). Presence and progression of white matter hyperintensities and cognition a meta-analysis. Neurology 82, 2127–2138.

Google Scholar

Körtgen, M., Park, G. J., Novotni, M., and Klein, R. (2003). “3D shape matching with 3D shape contexts,” in Proceedings of the 7th central European seminar on computer graphics, Budmerice, 5–17.

Google Scholar

Lee, S., Viqar, F., Zimmerman, M. E., Narkhede, A., Tosto, G., Benzinger, T. L., et al. (2016). White matter hyperintensities are a core feature of Alzheimer’s disease: Evidence from the dominantly inherited Alzheimer network. Ann. Neurol. 79, 929–939. doi: 10.1002/ana.24647

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., Zeng, N., Wu, P., and Clawson, K. (2022). Cov-net: A computer-aided diagnosis method for recognizing COVID-19 from chest X-ray images via machine vision. Expert Syst. Appl. 207:118029. doi: 10.1016/j.eswa.2022.118029

PubMed Abstract | CrossRef Full Text | Google Scholar

Lian, Z., Godil, A., Bustos, B., Daoudi, M., Hermans, J., Kawamura, S., et al. (2013). A comparison of methods for non-rigid 3D shape retrieval. Pattern Recogn. 46, 449–461.

Google Scholar

Liu, L., Fieguth, P., Guo, Y., Wang, X., and Pietikäinen, M. (2017). Local binary features for texture classification: Taxonomy and experimental study. Pattern Recogn. 62, 135–160.

Google Scholar

Loizou, C., Petroudi, S., Seimenis, I., Pantziaris, M., and Pattichis, C. (2015). Quantitative texture analysis of brain white matter lesions derived from T2-weighted MR images in MS patients with clinically isolated syndrome. J. Neuroradiol. 42, 99–114. doi: 10.1016/j.neurad.2014.05.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Luders, E., Thompson, P. M., Narr, K., Toga, A. W., Jancke, L., and Gaser, C. (2006). A curvature-based approach to estimate local gyrification on the cortical surface. Neuroimage 29, 1224–1230. doi: 10.1016/j.neuroimage.2005.08.049

PubMed Abstract | CrossRef Full Text | Google Scholar

Mak, L., Grandison, S., and Morris, R. J. (2008). An extension of spherical harmonics to region-based rotationally invariant descriptors for molecular shape description and comparison. J. Mol. Graph. Model. 26, 1035–1045. doi: 10.1016/j.jmgm.2007.08.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Makkinejad, N., Schneider, J. A., Yu, J., Leurgans, S. E., Kotrotsou, A., Evia, A. M., et al. (2019). Associations of amygdala volume and shape with transactive response DNA-binding protein 43 (TDP-43) pathology in a community cohort of older adults. Neurobiol. Aging 77, 104–111. doi: 10.1016/j.neurobiolaging.2019.01.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Mikolajczyk, K., Zisserman, A., and Schmid, C. (2003). “Shape recognition with edge-based features,” in Proceedings of the British machine vision conference (BMVC’03), (Durham: The British Machine Vision Association), 779–788.

Google Scholar

Newton, B. D., Wright, K., Winkler, M. D., Bovis, F., Takahashi, M. I, Dimitrov, E., et al. (2017). Three-dimensional shape and surface features distinguish multiple sclerosis lesions from nonspecific white matter disease. J. Neuroimaging 27, 613–619. doi: 10.1111/jon.12449

PubMed Abstract | CrossRef Full Text | Google Scholar

Novotni, M., and Klein, R. (2003). “3D Zernike descriptors for content based shape retrieval,” in Proceedings of the 8th ACM symposium on solid modeling and applications, (New York, NY: ACM), 216–225.

Google Scholar

Novotni, M., and Klein, R. (2004). Shape retrieval using 3D Zernike descriptors. Comput. Aided Des. 36, 1047–1062.

Google Scholar

Papakostas, G. A., Boutalis, Y. S., Karras, D. A., and Mertzios, B. G. (2007). A new class of Zernike moments for computer vision applications. Inf. Sci. 177, 2802–2819.

Google Scholar

Poels, M. M., Zaccai, K., Verwoert, G. C., Vernooij, M. W., Hofman, A., van der Lugt, A., et al. (2012). Arterial stiffness and cerebral small vessel disease. Stroke 43, 2637–2642.

Google Scholar

Pratama, S. F., Muda, A. K., Choo, Y.-H., and Abraham, A. (2015). “A comparative study of 2D UMI and 3D Zernike shape descriptor for ATS drugs identification,” in pattern analysis, intelligent security and the internet of things, eds A. Abraham, A. K. Muda, and Y-H. Choo (Berlin: Springer), 237–249. Available online at: https://link.springer.com/book/10.1007/978-3-319-17398-6

Google Scholar

Reuter, M., Wolter, F.-E., and Peinecke, N. (2006). Laplace–Beltrami spectra as ‘Shape-DNA’of surfaces and solids. Comput. Aided Des. 38, 342–366. doi: 10.1016/j.cad.2009.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Rodriguez, J. J., and Aggarwal, J. (1990). Matching aerial images to 3-D terrain maps. IEEE Trans. Pattern Anal. Mach. Intell. 12, 1138–1149.

Google Scholar

Sael, L., Li, B., La, D., Fang, Y., Ramani, K., Rustamov, R., et al. (2008). Fast protein tertiary structure retrieval based on global surface shape similarity. Proteins 72, 1259–1273. doi: 10.1002/prot.22030

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmidt, R., Schmidt, H., Haybaeck, J., Loitfelder, M., Weis, S., Cavalieri, M., et al. (2011). Heterogeneity in age-related white matter changes. Acta Neuropathol. 122, 171–185.

Google Scholar

Shen, K. K., Fripp, J., Mériaudeau, F., Chételat, G., Salvado, O., Bourgeat, P., et al. (2012). Initiative, detecting global and local hippocampal shape changes in Alzheimer’s disease using statistical shape models. Neuroimage 59, 2155–2166.

Google Scholar

Styner, M., Lieberman, J. A., McClure, R. K., Weinberger, D. R., Jones, D. W., and Gerig, G. (2005). Morphometric analysis of lateral ventricles in schizophrenia and healthy controls regarding genetic and disease-specific factors. Proc. Natl. Acad. Sci. U.S.A. 102, 4872–4877. doi: 10.1073/pnas.0501117102

PubMed Abstract | CrossRef Full Text | Google Scholar

Styner, M., Oguz, I., Xu, S., Brechbühler, C., Pantazis, D., Levitt, J. J., et al. (2006). Framework for the statistical shape analysis of brain structures using SPHARM-PDM. Insight J. 242–250.

Google Scholar

Szegõ, G. (1939). Orthogonal polynomials. New York, NY: American Mathematical Society.

Google Scholar

Tahmasbi, A., Saki, F., and Shokouhi, S. B. (2011). Classification of benign and malignant masses based on Zernike moments. Comput. Biol. Med. 41, 726–735. doi: 10.1016/j.compbiomed.2011.06.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Teague, M. R. (1980). Image analysis via the general theory of moments. J. Opt. Soc. Am. 70, 920–930.

Google Scholar

Terriberry, T. B., Joshi, S. C., and Gerig, G. (2005). “Hypothesis testing with nonlinear shape models,” in Proceedings of the Biennial international conference on information processing in medical imaging, (Berlin: Springer), 15–26.

Google Scholar

Tibshirani, R., Walther, G., and Hastie, T. (2001). Estimating the number of clusters in a data set via the gap statistic. J. R. Stat. Soc. Series B Stat. Methodol. 63, 411–423.

Google Scholar

Toriwaki, J., and Yoshida, H. (2009). Fundamentals of three-dimensional digital image processing. Berlin: Springer Science & Business Media.

Google Scholar

Tseng, B. Y., Gundapuneedi, T., Khan, M., Diaz-Arrastia, R., Levine, B., Lu, H., et al. (2013). White matter integrity in physically fit older adults. Neuroimage 82, 510–516.

Google Scholar

Venkatraman, V., Sael, L., and Kihara, D. (2009a). Potential for protein surface shape analysis using spherical harmonics and 3D Zernike descriptors. Cell Biochem. Biophys. 54, 23–32. doi: 10.1007/s12013-009-9051-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Venkatraman, V., Yang, Y. D., Sael, L., and Kihara, D. (2009b). Protein-protein docking using region-based 3D Zernike descriptors. BMC Bioinformatics 10:407. doi: 10.1186/1471-2105-10-407

PubMed Abstract | CrossRef Full Text | Google Scholar

Vranic, D., and Saupe, D. (2001). 3D shape descriptor based on 3D Fourier transform. EURASIP J. Adv. Signal Process. 271–274.

Google Scholar

Wachinger, C., Golland, P., Kremen, W., Fischl, B., and Reuter, M., and Alzheimer’s Disease Neuroimaging Initiative. (2015). Initiative, brainprint: A discriminative characterization of brain morphology. Neuroimage 109, 232–248. doi: 10.1016/j.neuroimage.2015.01.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Wachinger, C., Salat, D. H., Weiner, M., and Reuter, M., and Alzheimer’s Disease Neuroimaging Initiative. (2016). Initiative, Whole-brain analysis reveals increased neuroanatomical asymmetries in dementia for hippocampus and amygdala. Brain 139, 3253–3266. doi: 10.1093/brain/aww243

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, K., Zhu, T., and Wang, J. (2018a). Real-time terrain matching based on 3D Zernike moments. J. Navig. 71, 1441–1459.

Google Scholar

Wang, K., Zhu, T., Gao, Y., and Wang, J. (2018b). Efficient terrain matching with 3-D Zernike moments. IEEE Trans. Aerosp. Electron. Syst. 55, 226–235.

Google Scholar

Wardlaw, J., Smith, E., Biessels, G., Cordonnier, C., Fazekas, F., Frayne, R., et al. (2013). Neuroimaging standards for research into small vessel disease and its contribution to ageing and neurodegeneration. Lancet Neurol. 12, 822–838.

Google Scholar

Wee, C.-Y., and Paramesran, R. (2007). On the computational aspects of Zernike moments. Image Vision Comput. 25, 967–980.

Google Scholar

Wen, W., and Sachdev, P. (2004). The topography of white matter hyperintensities on brain MRI in healthy 60-to 64-year-old individuals. Neuroimage 22, 144–154. doi: 10.1016/j.neuroimage.2003.12.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye, B., and Chen, H. F. (2012). “3D terrain matching algorithm based on 3D Zernike moments,” in Proceedings of the 2012 symposium on photonics and optoelectronics, (Piscataway, NJ: IEEE), 1–4.

Google Scholar

Zeng, N., Li, H., and Peng, Y. (2021). A new deep belief network-based multi-task learning for diagnosis of Alzheimer’s disease. Neural Comput. Appl. 1–12.

Google Scholar

Zhang, L., da Fonseca, M. J., and Ferreira, A. (2007). Survey on 3D shape descriptors. Technical Report: DecorAR (FCT POSC/EIA/59938/2004). Lisboa: Fundação para a Ciência e a Tecnologia.

Google Scholar