# ARTIFICIAL INTELLIGENCE FOR MEDICAL IMAGE ANALYSIS OF NEUROIMAGING DATA

EDITED BY : Nianyin Zeng, Siyang Zuo, Guoyan Zheng, Yangming Ou and Tong Tong PUBLISHED IN : Frontiers in Neuroscience, Frontiers in Computational Neuroscience and Frontiers in Psychiatry

#### Frontiers eBook Copyright Statement

The copyright in the text of individual articles in this eBook is the property of their respective authors or their respective institutions or funders. The copyright in graphics and images within each article may be subject to copyright of other parties. In both cases this is subject to a license granted to Frontiers. The compilation of articles constituting this eBook is the property of Frontiers.

Each article within this eBook, and the eBook itself, are published under the most recent version of the Creative Commons CC-BY licence. The version current at the date of publication of this eBook is CC-BY 4.0. If the CC-BY licence is updated, the licence granted by Frontiers is automatically updated to the new version.

When exercising any right under the CC-BY licence, Frontiers must be attributed as the original publisher of the article or eBook, as applicable.

Authors have the responsibility of ensuring that any graphics or other materials which are the property of others may be included in the CC-BY licence, but this should be checked before relying on the CC-BY licence to reproduce those materials. Any copyright notices relating to those materials must be complied with.

Copyright and source acknowledgement notices may not be removed and must be displayed in any copy, derivative work or partial copy which includes the elements in question.

All copyright, and all rights therein, are protected by national and international copyright laws. The above represents a summary only. For further information please read Frontiers' Conditions for Website Use and Copyright Statement, and the applicable CC-BY licence.

ISSN 1664-8714 ISBN 978-2-88963-826-0 DOI 10.3389/978-2-88963-826-0

### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office<EMAIL_ADDRESS>
# ARTIFICIAL INTELLIGENCE FOR MEDICAL IMAGE ANALYSIS OF NEUROIMAGING DATA

Topic Editors: Nianyin Zeng, Xiamen University, China Siyang Zuo, Tianjin University, China Guoyan Zheng, Shanghai Jiao Tong University, China Yangming Ou, Harvard Medical School, United States Tong Tong, Independent researcher, China

Citation: Zeng, N., Zuo, S., Zheng, G., Ou, Y., Tong, T., eds. (2020). Artificial Intelligence for Medical Image Analysis of NeuroImaging Data. Lausanne: Frontiers Media SA. doi: 10.3389/978-2-88963-826-0

# Table of Contents

*05 Editorial: Artificial Intelligence for Medical Image Analysis of Neuroimaging Data*

Nianyin Zeng, Siyang Zuo, Guoyan Zheng, Yangming Ou and Tong Tong

*07 Cerebral Micro-Bleeding Detection Based on Densely Connected Neural Network*

Shuihua Wang, Chaosheng Tang, Junding Sun and Yudong Zhang

Jiahang Xu, Fangyang Jiao, Yechong Huang, Xinzhe Luo, Qian Xu, Ling Li, Xueling Liu, Chuantao Zuo, Ping Wu and Xiahai Zhuang

*40 Brain White Matter Hyperintensity Lesion Characterization in T2 Fluid-Attenuated Inversion Recovery Magnetic Resonance Images: Shape, Texture, and Potential Growth*

Chih-Ying Gwo, David C. Zhu and Rong Zhang

*55 Nested Dilation Networks for Brain Tumor Segmentation Based on Magnetic Resonance Imaging*

Liansheng Wang, Shuxin Wang, Rongzhen Chen, Xiaobo Qu, Yiping Chen, Shaohui Huang and Changhua Liu

*69 Diagnosis of Alzheimer's Disease via Multi-Modality 3D Convolutional Neural Network*

Yechong Huang, Jiahang Xu, Yuncheng Zhou, Tong Tong, Xiahai Zhuang and the Alzheimer's Disease Neuroimaging Initiative (ADNI)

*81 Supervised Brain Tumor Segmentation Based on Gradient and Context-Sensitive Features*

Junting Zhao, Zhaopeng Meng, Leyi Wei, Changming Sun, Quan Zou and Ran Su

*92 Prediction and Classification of Alzheimer's Disease Based on Combined Features From Apolipoprotein-E Genotype, Cerebrospinal Fluid, MR, and FDG-PET Imaging Biomarkers*

Yubraj Gupta, Ramesh Kumar Lama, Goo-Rak Kwon and the Alzheimer's Disease Neuroimaging Initiative

Chuanlu Lin, Yi Wang, Tianfu Wang and Dong Ni

# Editorial: Artificial Intelligence for Medical Image Analysis of Neuroimaging Data

#### Nianyin Zeng1 \*, Siyang Zuo2 , Guoyan Zheng3 , Yangming Ou4 and Tong Tong5

*1 Department of Instrumental and Electrical Engineering, Xiamen University, Xiamen, China, 2 Key Laboratory of Mechanism Theory and Equipment Design, Ministry of Education, Tianjin University, Tianjin, China, 3 School of Biomedical Engineering, Institute of Medical Robotics, Shanghai Jiao Tong University, Shanghai, China, 4 Department of Radiology, and Computational Health Informatics Program, Boston Children's Hospital, Boston, MA, United States, 5 College of Physics and Information Engineering, Fuzhou University, Fuzhou, China*

Keywords: medical image analysis, artificial intelligence, machine learning, pattern recognition, computational intelligence

### **Editorial on the Research Topic**

### **Artificial Intelligence for Medical Image Analysis of Neuroimaging Data**

With the development of advanced medical imaging techniques, a huge amount of medical images have been produced in various healthcare institutes and hospitals. Especially, there is a growing research interest in a more multidisciplinary approach for investigating brain structure and function in living humans and animals. In order to better interpret brain images, there is an increasing demand to introduce artificial intelligence methods such as machine learning, expert systems, robotics and perception, and evolutionary computation to automatically exploit useful information besides visual features. It should be pointed out that brain images themselves exhibit several distinguishing features that add to the difficulties in their analysis. In recent years, there have been many new research achievements in each aspect of artificial intelligence for brain image analysis. This Research Topic sought original contributions that address the challenges of artificial intelligence for brain image analysis and welcomed researchers in this field to share their experiences and new research achievements.

We were pleased to receive many submissions from authors of their latest research results on artificial intelligence methods for medical image analysis. Nineteen papers are finally accepted from a total of 29 submissions after rigorous reviews. They were contributed from different countries and regions, including China, the United Kingdom, the United States, Germany, South Korea, Denmark, Canada, and more.

Here, a brief introduction of the 19 accepted papers is given. We refer the readers to the papers in this topic and the references therein for more details. Lin W. et al. established a deep learning approach based on convolutional neural networks (CNN) to accurately predict MCI-to-AD conversion with magnetic resonance imaging (MRI) data. Kazeminejad and Sotero introduced a new biomarker extraction pipeline for Autism Spectrum Disorder that relies on the use of graph-theoretical metrics of fMRI-based functional connectivity to inform a support vector machine. Bi et al. proposed an advanced method, namely an evolutionary weighted random support vector machine cluster, for analysis of Alzheimer's disease. Ladefoged et al. focused on the problem of attenuation correction of PET/MRI in pediatric brain tumor patients based on a deep learning method. Livne et al. established a U-Net deep learning framework for high-performance vessel segmentation in patients with cerebrovascular disease. Wang, Sun et al. proposed a 14-layer convolutional neural network for the identification of multiple sclerosis. Huang C. et al. developed

Edited and reviewed by: *Vince D. Calhoun, Georgia State University, United States*

> \*Correspondence: *Nianyin Zeng<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Brain Imaging Methods, a section of the journal Frontiers in Neuroscience*

Received: *24 November 2019* Accepted: *17 April 2020* Published: *21 May 2020*

#### Citation:

*Zeng N, Zuo S, Zheng G, Ou Y and Tong T (2020) Editorial: Artificial Intelligence for Medical Image Analysis of Neuroimaging Data. Front. Neurosci. 14:480. doi: 10.3389/fnins.2020.00480*

**5**

a new fusion method based on the combination of the shuffled frog leaping algorithm and a pulse coupled neural network for the fusion of SPECT images and CT images to improve the quality of fused brain images. Xin et al. utilized a deep learning method to find differences between the brains of men and women. Zhang Y. et al. proposed an improved wavelet threshold for image de-noising. Lin C. et al. proposed a novel low-rank method for the simultaneous recovery and segmentation of pathological MR brain images. Zhang Z. et al. developed a multiscale time-series model for the diagnosis of brain diseases. Gupta et al. proposed a novel machine learning-based framework to discriminate subjects with AD or MCI, utilizing a combination of four different biomarkers. Zhao et al. proposed a supervised brain tumor segmentation method based on gradient and contextsensitive features. Huang Y. et al. developed a multi-modality 3D convolutional neural network for the diagnosis of Alzheimer's disease. Wang L. et al. presented the use of Nested Dilation Networks for brain tumor segmentation. Gwo et al. developed a method to characterize and quantify the shape, texture, and potential growth of white matter hyperintensity lesions. Xu et al. introduced a fully automatic framework for Parkinson's disease diagnosis. Wang, Xie et al. proposed an AlexNet transfer learning model for alcoholism identification. Wang, Tang et al. developed a densely connected neural network for analysis of cerebral micro-bleeding.

In the end, we strongly hope that this Research Topic will attract more research attention to artificial intelligence methods for medical image analysis. We thank the reviewers for their efforts to guarantee the high quality of this collection. We also thank all of the authors who have contributed.

## AUTHOR CONTRIBUTIONS

NZ wrote the editorial. SZ, GZ, YO, and TT edited the editorial.

### FUNDING

This work was supported in part by the International Science and Technology Cooperation Project of Fujian Province of China under Grant 2019I0003, in part by the UK-China Industry Academia Partnership Programme under Grant UK-CIAPP-276, in part by the Korea Foundation for Advanced Studies, in part by the Fundamental Research Funds for the Central Universities under Grant<PHONE_NUMBER>9, in part by The Open Fund of Provincial Key Laboratory of Eco-Industrial Green Technology-Wuyi University, in part by the Open Fund of Engineering Research Center of Big Data Application in Private Health Medicine of Fujian Province University under Grant KF2020002, and in part by the Harvard Medical School/Boston Children's Hospital Faculty Development Award (YO), and the St Baldrick Research Scholar Award (YO).

**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.

Copyright © 2020 Zeng, Zuo, Zheng, Ou and Tong. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Cerebral Micro-Bleeding Detection Based on Densely Connected Neural Network

#### Shuihua Wang1 \* † , Chaosheng Tang1† , Junding Sun1 \* and Yudong Zhang1,2 \* †

1 School of Computer Science and Technology, Henan Polytechnic University, Jiaozuo, China, 2 Department of Informatics, University of Leicester, Leicester, United Kingdom

#### Edited by:

Nianyin Zeng, Xiamen University, China

#### Reviewed by:

Haiou Liu, Friedrich-Alexander-Universität Erlangen-Nürnberg, Germany Zhi Yong Zeng, Fujian Normal University, China

#### \*Correspondence:

Shuihua Wang<EMAIL_ADDRESS>Junding Sun<EMAIL_ADDRESS>Yudong Zhang<EMAIL_ADDRESS>†These authors have contributed

equally to this work

### Specialty section:

This article was submitted to Brain Imaging Methods, a section of the journal Frontiers in Neuroscience

Received: 14 February 2019 Accepted: 12 April 2019 Published: 17 May 2019

### Citation:

Wang S, Tang C, Sun J and Zhang Y (2019) Cerebral Micro-Bleeding Detection Based on Densely Connected Neural Network. Front.  For example, CMBs can result in long-term disability, neurologic dysfunction, cognitive impairment and side effects from other medications and treatment. Therefore, it is important and essential to detect CMBs timely and in an early stage for prompt treatment. In this research, because of the limited labeled samples, it is hard to train a classifier to achieve high accuracy. Therefore, we proposed employing Densely connected neural network (DenseNet) as the basic algorithm for transfer learning to detect CMBs. To generate the subsamples for training and test, we used a sliding window to cover the whole original images from left to right and from top to bottom. Based on the central pixel of the subsamples, we could decide the target value. Considering the data imbalance, the cost matrix was also employed. Then, based on the new model, we tested the classification accuracy, and it achieved 97.71%, which provided better performance than the state of art methods.

Keywords: DenseNet, CMB detection, transfer learning, cost matrix, deep learning

## INTRODUCTION

Cerebral micro-bleeding (CMB) are small chronic brain hemorrhages that can be caused by structural abnormalities of the small vessels of the brain. According to the recent research reports, the causes of CMBs also can be some other common reasons, including high blood pressure, head trauma, aneurysm, blood vessel abnormalities, liver disease, blood or bleeding disorders and brain tumors (Martinez-Ramirez et al., 2014). It also can be caused by some unusual etiologies, such as cocaine abuse, posterior reversible encephalopathy, brain radiation therapy, intravascular lymphomatosis, thrombotic thrombocytopenic purpura, moyamoya disease, infective endocarditis, sickle cell anemia/β-thalassemia, proliferating angio-endotheliomatosis, cerebral autosomal dominant arteriopathy subcortical infarcts, leukoencephalopathy (CADASIL), genetic syndromes, or obstructive sleep apnea (Noorbakhsh-Sabet et al., 2017). The patients suffering from CMBs can have symptoms where the corresponding area that is controlled by the bleeding area malfunctions, resulting in a rise in intracranial pressure due to the large mass putting pressure on the brain and so on. CMBs could be easily ignored as the similar symptoms and signs of the subarachnoid hemorrhages, unless the patients have more obvious symptoms, such as a headache followed by vomiting. Those symptoms can eventually become worse or occur suddenly, based on the distribution and intensity of the CMBs. Patients suffering from CMBs can result in cognitive impairment, neurologic dysfunction and long-term disability. CMBs could also induce side effects

from medication or treatments. The worse thing is that the death is possible and can happen quickly. Therefore, the early and prompt diagnosis of CMBs is essential and helpful in timely medical treatment.

Due to the paramagnetic susceptibility of the hemosiderin (Allen et al., 2000), CMBs can be visualized by T2∗ -gradient recalled echo (GRE) imaging or susceptibility weighted imaging (SWI). Traditionally, CMBs are manually interpreted based on criteria including shapes, diameters and signal characteristics after imaging. However, the criteria were varied as reported in different studies (Cordonnier et al., 2007), until 2009 when Greenberg et al. (2009) published the consensus on standard criteria for CMB identification. However, manual detection methods involve the human interventions, which can bring biases. Meanwhile, the manual detection is labor intensive, hard to reproduce and difficult to exclude the mimics, which can lead to misdiagnosis.

Therefore, the development of automatic CMB detection is important and essential for the accurate detection and early diagnosis of CMBs. Due to the benefits of advanced imaging technologies, massive computer vision aided systems have been developed for automatic CMB detection. For example, Ateeq et al. (2018) proposed a system based on an ensemble classifier. Their system consisted of three steps: first the brain was extracted, then the initial candidates were detected based on the filter and threshold, and finally, feature extraction and classification model were built to remove the false alarms. Fazlollahi et al. (2015) proposed using a multi-scale Laplacian of Gaussian (msLoG) technique to detect the potential CMB candidates, followed by extracting a set of 3-dimensional Radon and Hessian based shape descriptors within each bounding box to train a cascade of binary random forests. Barnes et al. (2011) proposed a statistical thresholding algorithm to recognize the potential hypointensities. Then, a supervised classification model based on the support vector machine was employed to distinguish true CMBs from other marked hypo-intensities. van den Heuvel et al. (2016) proposed an automatic detection system for microbleeds in MRIs of patients with trauma based on twelve characteristics related with the dark and spherical characteristics of CMBs and the random forest classifier. Bian et al. (2013) proposed a 2D fast radial symmetry transform (RST) based method to roughly detect the possible CMBs. Then the 3D region growing on the possible CMBs was utilized to exclude the falsely identified CMBs. Ghafaryasl et al. (2012) proposed a computer aided system based on following three steps: skull-stripping, initial candidate selection and reduction of false-positives (FPs) by a two-layer classifier. Zhang et al. (2017) proposed voxel-vise detection based on a single hidden layer feed-forward neural network with scaled conjugate gradient. Chen (2017) proposed a seven-layer deep neural network based on the sparse autoencoder for voxel detection of CMBs. Seghier et al. (2011) proposed a system named MIDAS for automatic CMB detection.

All above methods have reached great progress in CMB detection. However, their detection accuracy and robustness are still in need of improvement.

Therefore, in this paper, we employed the SWI for CMB imaging, which was because SWI could provide high resolution as reported in Haacke et al. (2009) and work as the most sensitive techniques to visualize CMBs. Considering the limited amounts of labeled images, and knowledge to recognize representative characters about the medical images, we considered utilizing the DenseNet as the basic algorithm for transfer learning. The reason for this is because the amount of labeled CMB images is typically very limited, and it is hard to effectively train a classifier to get high detection accuracy. In summary, we proposed using transfer learning of DenseNet for CMB detection based on the collected images, which means we use the knowledges obtained from training the related tasks by DenseNet for CMB detection.

The remainder of this paper is organized in a structure as follows: "Materials and Methods" section describes the method used in this research, "Transfer Learning" section explains why we employed the transfer learning, "CMB Detection Based on the Transfer Learning and DenseNet" section describes the research materials used in this paper, including the training set and test set, and also offers the experiment results, and finally, "Discussion" section provides the conclusion and discussion.

 Traditionally, DL consists of multiple layers of nonlinear processing units to obtain the features. The cascaded layers take the output from their previous layer as input. In order to explore the potential of DL, many researchers tried to make the network deeper and wider. However, it suffers from either exploding or the vanishing gradient problem (VGP). Therefore, multiple different structures of DL were proposed. For example, AlexNet, the winner of ImageNet Large Scale Visual Recognition Competition (ILSVRC) 2012, was proposed by Krizhevsky et al. (2012) and has the same structure as LeNet but has max pooling and ReLU non-linearity. VGGNet, proposed by Karen Simonyan (2015), won the second place in ILSVRC 2014 and consisted of deeper networks (19 layers) compared to AlexNet. GoogLeNet, the winner of ILSVRC 2014, provides a deeper and wider network that incorporates 1 × 1 convolutions to reduce the dimensionality of the feature maps before expensive convolutions, together with parallel paths with different receptive field sizes to capture sparse patterns of correlations in the feature map stacks. ResNet, the winner of ILSVRC 2015, offers a 152-layer network that introduces a skip at least two layers or shortcut connections (He et al., 2016). Huang et al. (2017) proposed DenseNet where each layer takes the output from all previous layers in a feed-forward fashion and offers L(L+1)/2 connections for L layers, while traditional convolution networks with L layers provide L connections. According to the report in Huang et al. (2017), DenseNet

can beat the state-of-the-art ResNet structure on ImageNet classification task.

Considering the outstanding performance of DenseNet, we proposed employing DenseNet for cerebral microbleed detection in this paper. The detail of DenseNet was introduced as follows. However, before providing the illustration of the DenseNet, we would first introduce the traditional convolution neural network (CNN) and figure out the difference between CNN and DenseNet later.

### Traditional Convolution Neural Network

The traditional CNN usually includes convolution layer, ReLU Layer, pooling layer, fully connected layer and softmax layer (Zeng et al., 2014, 2016a,c, 2017a,b). The functions of different layers are introduced as follows:

Convolution layer works as the core session of a CNN. The feature maps are generated via the convolution of the input with different kernels. Mathematically, it can be expressed as **Figure 1**, which shows a toy example of convolution operation.

Then, following the convolution layer, we have non-linear activation function, named ReLU, which works to obtain the nonlinear features. The purpose of the ReLU layer is to introduce non-linearity into the network. The mathematic expression of ReLU is shown as Eq. 1:

$$f(\mathbf{x}) = \mathbf{x}^+ = \max(0, \mathbf{x}) \tag{1}$$

The pooling layer works by resizing the feature maps spatially to decrease the number of parameters, memory footprint and to make the computation less intensive in the network. The pooling function works on each feature map, the main approaches used for pooling are max pooling as Eq. 2, average pooling as Eq. 3:

$$a\_j = \max\_{i \in R\_f} (M\_i) \tag{2}$$

$$a\_j = \frac{1}{|R\_j|} \sum\_{i \in R\_j} M\_i \tag{3}$$

In which M stands for the pooling region and Rj represents for the number of elements within the pooling region.

Fully connected layers will calculate the confidential scores, which are stored in a volume of size 1 × 1 × n. Here, n means the number of categories, and each element stands for class scores.

Every neuron of the fully connected layer is connected to all the neurons in the earlier layers.

### Structure Revision of the CNN

In the traditional CNN, all layers are connected gradually as in Eq. 4:

$$\mathbf{x}\_{l} = H\_{l}(\mathbf{x}\_{l-1}) \tag{4}$$

However, as the network becomes deeper and wider, the networks may suffer from either exploding or gradient vanishing. Therefore, researchers proposed different network structures to overcome this problem. For example, ResNet revised this behavior by short connection, and the equation is reformulated as (5).

$$\mathbf{x}\_{l} = H\_{l}(\mathbf{x}\_{l-1}) + \mathbf{x}\_{l-1} \tag{5}$$

Instead of making the sum of the output feature maps of the layer with the incoming feature maps, DenseNet concatenates them sequentially. The expression is reformulated into Eq. 6:

$$\mathbf{x}\_{l} = H\_{l}(\left[\mathbf{x}\_{0}, \boldsymbol{\varkappa}\_{1}, \boldsymbol{\varkappa}\_{2}, \dots, \boldsymbol{\varkappa}\_{l-1}\right])\tag{6}$$

In which l means the index of the layer number, H stands for a non-linear operation and xl stands for the output of the lth layer.

### DenseNet

As expressed in Eq. 6, DenseNet introduces straight forward connections from any layers to all following layers. In other words, the lth layer receives feature-maps from all previous l – 1 layers. However, if the feature maps' size changes, the concatenation operation is not feasible. Therefore, downsampling to change the size of the feature maps are introduced. In order to make the down-sampling in the structure of DenseNet possible, multiple densely connected dense blocks are introduced to divide the network. The layers between the blocks are named as transition layers that have batch normalization, convolution and pooling operations, as shown in **Figure 2**. **Figure 2** describes a case of DenseBlock, in which the layer number is 5 and the growth rate is set as k. Each layer receives feature maps from all earlier layers.

For each operation Hl , it generates k feature maps, which is defined as growth rate. Therefore, the l th layer will have k0 + k(l − 1) feature maps, and k0 is the number of channels in the

FIGURE 5 | CMB samples.

input layer. As the network typically has a large number of inputs, a 1 × 1 convolution is employed as the bottleneck layer before the 3 × 3 convolution layer to reduce the feature maps and improve the computation efficiency.

To further compress the model to improve the model compactness, the feature maps are further reduced by the transition layer. For example, if a dense block generates m feature maps and the compression factor is set as θ ∈ (0, 1], then the feature maps will be reduced to bθmc via the followed transition layer. If θ = 1, the number of feature maps will be the same. **Figure 3** shows the structure of DenseNet, which is composed of three DenseBlocks, an input layer and transition layers. The cropped samples are used as the input, the final layer will tell us whether it is CMB or Non-CMB in this research.

### TRANSFER LEARNING

DenseNet has been widely applied in the medical research. For example, Gottapu and Dagli (2018) proposed using DenseNet for Anatomical Brain Segmentation. Khened et al. (2019) proposed cardiac segmentation based on fully convolutional multi-scale residual DenseNets. Wang H. et al. (2019) offered a system for recognition of mild cognitive impairment (MCI) and Alzheimer's disease (AD), based on the ensemble of 3D densely connected convolution network. Considering the limited amounts of labeled training samples, it is far way from enough to retrain the whole network of DenseNet from scratch to get a high classification accuracy. Therefore, in this paper, we proposed transfer learning, which means frozen the earlier layers and retrain the later layers of DenseNet for CMB detection task. The structure of DenseNet used here is DenseNet 201.

In order to make the pretrained DenseNet 201 for CMB detection feasible, which was a binary classification of CMB or non-CMB, the fully connected (FC) layer with

TABLE 1 | Dividing of the dataset for training and testing.

1000 neuron was replaced by a new FC layer with 2 neurons. The structure of the remaining part of DenseNet 201 was unchanged.

### CMB DETECTION BASED ON THE TRANSFER LEARNING AND DENSENET

### Materials

The subjects used in this research are ten healthy controls and ten patients of CADASIL. Twenty 3D volumetric images were obtained from the 20 patients. Then, Software Sygno MR B17 was utilized to rebuild the 3D volumetric image. Each 3D volumetric image's size is uniformly set as 364∗ 448∗ 48.

In order to mark the CMBs from the subjects manually, we employed three neuron-radiologists with more than twentyyears' experience. The rules were set as follows: (1) via tracking the neighbor slices, blood vessels were first excluded, (2) lesions should be smaller than 10 mm in diameter. The potential CMBs were labeled as either "possible" or "Definite," Otherwise, regarded as non-CMB voxels. In case of the conflictions, we proposed to obey the rule that the minority should be subordinate to the majority.

The sample images were generated from the original image. We applied the sliding window whose size is set as 61 by 61 to the original image.  All the pixels located within a sliding window were used as one input, and the point located in the center of the sliding window was used as the target value. It means that if the central pixel is true or 1, then the target value is 1, otherwise, the

**12**

target label is set as 0. It is expressed in the Eqs 7 and 8:

$$I = \mathcal{W}(\mathfrak{p})\tag{7}$$

$$\text{O}u = \begin{cases} 1, & \text{Central pixel } p \text{ is true (CMB)}\\ 0, & \text{Central pixel } p \text{ is false (non-CMB)} \end{cases} \tag{8}$$

Where I stands for the cropped sample images generated via the sliding window, p represents for the central pixel, W(p) means the pixels which centered on pixel p and were located inside the sliding window, and Ou means the label value. **Figures 4**, **5** show the sample of CMB and non-CMB centered images.

The sliding window was supposed to cover the image from left to right and top to bottom with the stride size as 1. Therefore, we got the total CMB voxels as 68, 847 and non-CMB voxels as 56, 582, 536. The training and test set was divided as **Table 1**. We randomly selected 10000 images for each category of the test, and the remaining images were used for training.

To make the images suitable for DenseNet, which should be resized as 224 × 224 × 3, we padded the images with zero. The preprocessed image sample is shown as **Figure 6**. Then, **Figure 7** shows the flowchart of the DenseNet, including number of feature maps generated by each layer.

From **Table 1**, we can find that the Non-CMB training data dominates the majority type CMB, which will cause the classifier more bias toward to the Non-CMB. Therefore, it may cause difficulties in controlling false positives and false negatives, which means the model is hard to find the CMB samples. Therefore, in order to overcome this challenge, we introduced cost matrix (Wu and Kim, 2010). The cost ratio ct was set as 961 via Eq. 9:

$$ct = \mathcal{N}\_{non-CMB} / \mathcal{N}\_{CMB} \tag{9}$$

In which Nnon−CMB means the number of non-CMB training samples and NCMB stands for number of CMB training samples. The reason for why we employ the cost matrix instead of over sampling or down sampling is mainly because we have more concerns about the false positives and false negatives, therefore it is better to highlight the imbalanced learning problem by using cost matrices instead of creating a balanced data distribution forcefully.

### Experiment Design

The goal of this research is to identify the input image as either CMB or Non-CMB. In order to achieve this goal, we proposed using DenseNet 201 as the basic algorithm for transfer learning, based on the excellent performance of DenseNet on ImageNet classification task. Section "Materials" stated the materials used in this research. Based on the original images, we created 68, 847 CMB subsamples and 56, 582, 536 Non-CMB subsamples. 10000 samples were randomly selected as test samples. The remaining sub-samples were used for training. In order to overcome the problem of data imbalance, we proposed cost matrix to show the more concerns in false positive and false negatives. The experiment is carried on by Matlab on the Windows 10 Operation System with 2.88 GHz processor and 8 GB memory. The following experiments were carried out: (1) CMB detection based on DenseNet. The measurements used here include accuracy,

TABLE 2 | Confusion matrix of detected CMB and Non-CMB.

sensitivity and specificity. The definition of the measurements can be found in Zhang et al. (2018b). (2) Different cutting points of transfer learning. (3) In order to show the performance of proposed methods, we compared with other state of art work. Considering the measurements provided in other research, we only used sensitivity, specificity and accuracy.

In order to provide better illustration of DenseNet, we added a flowchart with feature map size, learnable weights of each layer. As we only noted the size of the width, the length should be same with the width.

### CMB Detection Result Based on DenseNet

The rebuilt network was composed of four DenseBlocks, one input layer, three transition layers, one fully connected layer with two neurons, a softmax layer and a classification layer, as described in **Figure 8**.

**Table 2** provides the detection result. The correctly detected CMBs were 9777, and for Non-CMB they were 9764. 236 non-CMBs were incorrectly detected as CMBs, and 223 CMBs were wrongly detected as non-CMBs. The sensitivity was achieved as 97.78%, the specificity was 97.64%, the accuracy was 97.71% and the precision was 97.65%. Above measurements were obtained based on the average of 10 runs as shown in **Table 3**.

### Comparison to the Different Cases of Transfer Learning

In order to achieve the best performance of transfer learning, different cutting points for transfer learning were designed as shown in **Figure 8**. Due to the limited subjects, we mainly focused on retraining the later layers of DenseNet. Therefore, in case A, the DenseNet 201 except for the last three layers, was used as the feature extractor for this research, and we retrained the newly added three layers.

In case B, C, and D, we included extra layers for retraining. For example, case B retrained the DenseBlock 4, Batch normalization, Average pooling, Fully connected (FC) layer, softmax layer and

TABLE 3 | Measurements value CMB detection based on transfer learning of DenseNet (Units: %).

#### TABLE 4 | Comparison of different cases of transfer learning (Unit: %).

classification layer. It was implemented via setting the learning rate to 0 for earlier layers and setting learning rate factor to 10 for layers to be retrained. **Table 4** illustrates the comparison results.

From **Table 4**, we can find that Case A performed slightly better than the other three cases in the terms of sensitivity and

TABLE 5 | Comparison to the state of art methods.

accuracy. Considering that in medical research, we focus more on the sensitivity and accuracy than on the other two terms, we thought Case A provided the best performance among all the cases. **Figure 9** shows the error bar of the measurement values. From the point of storage consuming, all four cases take about the same RAM as we did when not using the precomputation method.

### Comparison to the State of Art Work

In order to validate our proposed method, we compared different state of the art methods, including traditional machine learning methods and DL methods.

From **Table 5**, we compared our method with single-hidden layer feed-forward neural-network (SLFN)+ (leaky) rectified linear unit, 4-layer sparse encoder, 7-layer sparse encoder, different layers of CNN, Naive Bayesian Classifier and so on. We can find that our proposed method offers the best performance. DenseNet works as a logical extension of ResNet but provides more compact models and fully uses the features.

**Figure 10** shows the bar chart of the comparison of the state of the state of art methods. It shows that our proposed method performs slightly better than the current best method, but largely improved compared to the traditional method naïve Bayes classifier (NBC).

### DISCUSSION

In this paper, we proposed to employ DenseNet to detect CMBs in patients with CADASIL. DenseNet was proposed by Huang et al. (2017) and competed the other DL methods for ImageNet classification task because of its model compactness and fully used features. DenseNet are quite similar with ResNet, however, instead of the summation, DenseNet proposed the concatenation of all feature maps from previous layers, which encourages the feature reuse, the VGP alleviation, and the decreased number of parameters.

Therefore, in this paper, we proposed using DenseNet for CMB detection by supposing CMB detection has similarity with ImageNet classification. However, because of the data imbalance, we utilized cost matrix to avoid the model bias toward non-CMB, which means the model would be hard to find CMBs if trained under the imbalanced dataset. As there are some other

### REFERENCES

methods for data imbalance, such as over sampling and down sampling, we have more concerns about the false negatives or false positives. Therefore, instead of enforcing the data into balanced distribution, we employed the cost matrix. In order to check the best cutting point, we test different cases of transfer learning and the results are shown in **Table 4**, however, the difference is not so obvious. On the other hand, training less layers can help us save time and decrease the computation cost if we import the strategy of precomputation.

In the future, we will try to collect more samples and test more different structures for CMB detection. Meanwhile, the training cost long term is very high (Zeng et al., 2018), therefore it is necessary to optimize the algorithm to make the training fast (Zeng et al., 2016b,d). We will consider other precomputation and some optimization methods (Zhang et al., 2018a).

### DATA AVAILABILITY

The datasets for this manuscript are not publicly available because due to the privacy of the subjects. Requests to access the datasets should be directed to<EMAIL_ADDRESS>
### ETHICS STATEMENT

This research was approved by Institutional Review Board (IRB) of the First Affiliated Hospital of Nanjing Medical University. We obtained written informed consent from the participants of this study.

### AUTHOR CONTRIBUTIONS

SW proposed the study and wrote the draft. CT and JS designed the model and interpreted the results. CT and YZ analyzed the data. SW and YZ acquired the preprocessed data. All authors gave critical revision and consent for this submission.

### FUNDING

This project was financially supported by Natural Science Foundation of Jiangsu Province (No. BK20180727).

resonance images. Mag. Reson. Imag. 29, 844–852. doi: 10.1016/j.mri.2011. 02.028

**Conflict of Interest Statement:** 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.

Copyright © 2019 Wang, Tang, Sun and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Alcoholism Identification Based on an AlexNet Transfer Learning Model

Shui-Hua Wang1,2,3†, Shipeng Xie4†, Xianqing Chen5†, David S. Guttery 3† , Chaosheng Tang1 \*, Junding Sun1 \* and Yu-Dong Zhang1,3,6 \* †

*1 School of Computer Science and Technology, Henan Polytechnic University, Jiaozuo, China, 2 School of Architecture Building and Civil Engineering, Loughborough University, Loughborough, United Kingdom, 3 Department of Informatics, University of Leicester, Leicester, United Kingdom, 4 College of Telecommunications and Information Engineering, Nanjing University of Posts and Telecommunications, Nanjing, China, 5 Department of Electrical Engineering, College of Engineering, Zhejiang Normal University, Jinhua, China, 6 Guangxi Key Laboratory of Manufacturing System and Advanced Manufacturing Technology, Guilin, China*

Aim: This paper proposes a novel alcoholism identification approach that can assist radiologists in patient diagnosis.

Edited by: *Nianyin Zeng, Xiamen University, China*

### Reviewed by:

*Liansheng Wang, Xiamen University, China Chen Yang, Hangzhou Dianzi University, China*

#### \*Correspondence:

*Chaosheng Tang<EMAIL_ADDRESS>Junding Sun<EMAIL_ADDRESS>Yu-Dong Zhang<EMAIL_ADDRESS>
*†These authors have contributed equally to this work*

#### Specialty section:

*This article was submitted to Neurodegeneration, a section of the journal Frontiers in Psychiatry*

Received: *14 February 2019* Accepted: *21 March 2019* Published: *11 April 2019*

#### Citation:

*Wang S-H, Xie S, Chen X, Guttery DS, Tang C, Sun J and Zhang Y-D (2019) Alcoholism Identification Based on an AlexNet Transfer Learning Model. Front. Psychiatry 10:205. doi: 10.3389/fpsyt.2019.00205* Method: AlexNet was used as the basic transfer learning model. The global learning rate was small, at 10−4 , and the iteration epoch number was at 10. The learning rate factor of replaced layers was 10 times larger than that of the transferred layers. We tested five different replacement configurations of transfer learning.

Results: The experiment shows that the best performance was achieved by replacing the final fully connected layer. Our method yielded a sensitivity of 97.44%± 1.15%, a specificity of 97.41 ± 1.51%, a precision of 97.34 ± 1.49%, an accuracy of 97.42 ± 0.95%, and an F1 score of 97.37 ± 0.97% on the test set.

Conclusion: This method can assist radiologists in their routine alcoholism screening of brain magnetic resonance images.

Keywords: alcoholism, transfer learning, AlexNet, data augmentation, convolutional neural network, dropout, local response normalization, magnetic resonance imaging

### INTRODUCTION

Alcoholism (1) was previously composed of two types: alcohol abuse and alcohol dependence. According to current terminology, alcoholism differs from "harmful drinking" (2), which is an occasional pattern of drinking that contributes to increasing levels of alcohol-related ill-health. Today, it is defined depending on more than one of the following conditions: alcohol is strongly desired, usage results in social problems, drinking large amounts over a long time period, difficulty in reducing alcohol consumption, and usage resulting in non-fulfillment of everyday responsibilities.

Alcoholism affects all parts of the body, but it particularly affects the brain. The size of gray matter and white matter of alcoholism subjects are less than age-matched controls (3), and this shrinkage can be observed using magnetic resonance imaging (MRI). However, neuroradiological diagnosis using MR images is a laborious process, and it is difficult to detect minor alterations in the brain of alcoholic patient. Therefore, development of a computer vision-based automatic smart alcoholism identification program is highly desirable to assist doctors in making a diagnosis.

Within the last decade, studies have developed several promising alcoholism detection methods. Hou (4) put forward a novel algorithm called predator-prey adaptive-inertia chaotic particle swarm optimization (PAC-PSO), and applied it to identify alcoholism in MR brain images. Lima (5) proposed to use Haar wavelet transform (HWT) to extract features from brain images, and the authors used HWT to detect alcoholic patients. Macdonald (6) developed a logistic regression (LR) system. Qian (7) employed the cat swarm optimization (CSO) and obtained excellent results in the diagnosis of alcoholism. Han (8) used wavelet Renyi entropy (WRE) to generate a new biomarker; whereas Chen (9) used a support vector machine, which was trained using a genetic algorithm (SVM-GA) approach. Jenitta and Ravindran (10) proposed a local mesh vector co-occurrence pattern (LMCoP) feature for assisting diagnosis.

Recently, deep learning has attracted attention in many computer vision fields, e.g., synthesizing visual speech (11), liver cancer detection (12), brain abnormality detection (13), etc. As a result, studies are now focused on using deep learning techniques for alcoholism detection. Compared to manual feature extraction methods (14–18), deep learning can "learn" the features of alcoholism. For example, Lv (19) established a deep convolutional neural network (CNN) containing seven layers. Their experiments found that their model obtained promising results, and the stochastic pooling provided better performance than max pooling and average pooling. Moreover, Sangaiah (20) developed a ten-layer deep artificial neural network (i.e., three fully-connected layers and seven conv layers), which integrated advanced techniques, such as dropout and batch normalization, into their neural network.

Transfer learning (TL) is a new pattern recognition problemsolver (21–23). TL attempts to transfer knowledge learned using one or more source tasks (e.g., ImageNet dataset) and uses it to improve learning in a related target task (24). In perspective of realistic implementation, the advantages of TL compared to plain deep learning are: (i) TL uses a pretrained model as a starting point; (ii) fine-tuning a pretrained model is usually easier and faster than training a randomly-initialized deep neural network.

The contribution of this paper is that we may be the first to apply transfer learning in this field of alcoholism identification. We used AlexNet as the basic transfer learning model and tested different transfer configurations. Further, the experiments showed that the performance (sensitivity, specificity, precision, accuracy, and F1 score) of our model is >97%, which is superior to state-of-the-art approaches. We also validated the effectiveness of using data augmentation which further improves the performance of our model.

### DATA PREPROCESSING

### Datasets

This study was approved by the ethical committee of Henan Polytechnic University.   The division is shown in **Table 1**.

TABLE 1 | Dataset division into training, validation, and test sets.

TABLE 2 | Data augmentation.

### Data Augmentation

To improve the performance of deep learning, data augmentation (DA) (25) was introduced. This was done because our deep neural network model has many parameters, so we needed to show that our model contains a proportional amount of sample images to achieve optimal performance. For each original image, we generated a horizontally flipped image. Then, for both original and horizontal-flipped images, we applied the following five DA techniques: (i) noise injection, (ii) scaling, (iii) random translation, (iv) image rotation, and (v) gamma correction. Each of those methods produced 30 new images.

Gaussian noise with zero-mean and variance of 0.01 was applied to every image. Scaling was used with a scaling factor of 0.7–1.3, with an increase of 0.02. Random translation was utilized with a random shift within [−40 40] pixels. Image rotation with rotation angle varies from −30o to 30o and a step of 2o was employed. Gamma correction with gamma value varies from 0.4 to 1.6 with a step of 0.04 was utilized.

The DA result is shown in **Table 2**. Each image generated (1+30∗ 5)∗ 2 = 302 images including itself. After DA, the training set had 24,160 alcoholism brain images and 24,160 healthy brain images. Altogether, the new training set consisted of a balanced 160∗ 320 = 48,320 samples.

### METHODOLOGY

### Fundamentals of Transfer Learning

The core knowledge of transfer learning (TL) is shown in **Figure 1**. The core is to use a relatively complex and successful pre-trained model, trained from a large data source, e.g.,

ImageNet, which is the large visual database developed for visual object recognition research (26). It contains more than 14,000,000 hand-annotated images and at least one million images are provided with bounding boxes. ImageNet contains more than 20,000 categories (27). Usually, pretrained models are trained on a subset of ImageNet with 1,000 categories. Then we "transferred" the learnt knowledge to the relatively simplified tasks (e.g., classifying alcoholism and non-alcoholism in this study) with a small amount of private data.

Two attributes are important to help the transfer (28): (i) The success of the pretrained model can promote the exclusion of user intervention with the boring hyper-parameter tuning of new tasks; (ii) The early layers in pretrained models can be determined as feature extractors that help to extract low-level features, such as edges, tints, shades, and textures.

Traditional TL only retrains the new layers (29). In this study, we initially used the pretrained model, and then re-trained the whole structure of the neural network. Importantly, the global learning rate is fixed, and the transferred layers will have a low factor, while newly-added layers will have a high factor.

### AlexNet

AlexNet competed in the ImageNet challenge (30) in 2012, achieved a top-5 error of only 15.3%, more than 10.8% better than the result of the runner-up that used the shallow neural network. Original AlexNet was performed on two graphical processing units (GPUs). Nowadays, researchers tend to use only one GPU to implement AlexNet. **Figure 2** illustrates the structure of AlexNet. This study only counts layers associated with learnable weights. Hence, AlexNet contains five conv layers (CL) and three fullyconnected layers (FCL), totaling eight layers.

The details of learnable weights and biases of AlexNet are shown in **Table 3**. The total weights and biases of AlexNet are 60,954,656 + 10,568 = 60,965,224. In Matlab, the variable is stored in single-float type, taking four bytes for each variable. Hence, in total we needed to allocate 233 MB.

### Common Layers in AlexNet

Compared to traditional neural networks, there are several advanced techniques used in AlexNet. First, CLs contain a set of learnable filters. For example, the user has a 3D input with a size of PW×PH×D, a 3D filter with a size of QW×QH×D. As a consequence, the size of the output activation map is SW×SH. The value of SW and SH can be obtained by

$$S\_W = 1 + \frac{P\_W - Q\_W + 2\beta}{\mu} \tag{1}$$

$$S\_H = 1 + \frac{P\_H - Q\_H + 2\beta}{\mu} \tag{2}$$

where µ is the stride size and β represents the margin. Commonly, there may be T filters. One filter will generate one 2D feature map, and T filters will yield an activation map with a size of SW×SH×T. An illustration of convolutional procedure is shown in **Figure 3**. The "feature learning" in the filters here, can be regarded as a replacement of the "feature extraction" in traditional machine learning.

Second, the rectified linear unit (ReLU) function was employed to replace the traditional sigmoid function S(x) in terms of the activation function (31). The reason is because the sigmoid function may come across a gradient vanishing problem in deep neural network models.

$$S(\mathbf{x}) = \frac{1}{1 + \exp(-\mathbf{x})} \tag{3}$$

Therefore, the ReLU was proposed and defined as follows:

$$\text{ReLU}(\mathfrak{x}) = \max(0, \mathfrak{x}) \tag{4}$$

The gradient of ReLU is one at all times, when the input is larger than or equal to zero. Scholars have proven that the convergence speed of deep neural networks, with ReLU as the activation function, is 6x quicker than traditional activation functions. Therefore, the new ReLU function greatly accelerates the training procedure.

Third, a pooling operation is implemented with two advantages: (i) It can reduce the size of the feature map, and thus reduce the computation burden. (ii) It ensures that the representation becomes invariant to the small translation of the input. Map pooling (MP) is a common technique that chooses the maximum value among a 2 × 2 region of interest. **Figure 4** shows a toy example of MP, with a stride of 2 and kernel size of 2.

The fourth improvement is the "local response normalization (LRN)." Krizhevsky et al. (26) proposed the LRNs in order to aid generalization. Suppose that ai represents a neuron computed by applying kernel i and ReLU non-linearity, then the responsenormalized neuron bi will be expressed as:

$$b\_i = \frac{a\_i}{\left(m + \alpha \sum\_{s=\max(0, i-z/2)}^{\min(Z-1, i+z/2)} a\_s^2\right)^\beta} \tag{5}$$

where z is the window channel size, controlling the number of channels used for normalization of each element, and Z is the gross number of kernels in that layer. Hyperparameters are set as: β = 0.75, α = 10−4 , m = 1, and z = 5.

TABLE 3 | Learnable layers in AlexNet.

Fifth, the fully connected layers (FCLs) have connections to all activations in the previous layer, so they can be modeled as multiplying the input by a weight matrix and then adds a bias vector. The last fully-connected layer includes the equal number of artificial neurons as the number of total classes C. Therefore, each neuron in the last FCL represents the score of that cognate class, as shown in **Figure 5**.

Sixth, the softmax layer (SL), utilizes the multiclass generalization of logistic regression (32), also known as softmax function. SL is commonly connected after the final FCL. From the perspective of the activation function, the sigmoid/ReLU function works on a single input single output mode, while the SL serves as a multiple input multiple output mode, as shown in **Figure 6**. A toy example can be imagined. Suppose we have a four input at the final SL layer with values of (1–4), then after a softmax layer, we have an output of [0.032, 0.087, 0.236, 0.643].

Suppose that T(f) symbolizes the prior class probability of class f, and T(h|f) means the conditional probability of sample h given class f. Then we can conclude that the likelihood of sample h belonging to class f is

$$T(f \mid h) = \frac{T(h \mid f) \times T(f)}{\sum\_{i=1}^{F} T(h \mid i) \times T(i)} \tag{6}$$

Here F stands for the whole number of classes. Let Ωf equals

$$
\Omega\_f = \ln \left[ T(h, f) \times T(f) \right] \tag{7}
$$

Afterwards, we get

$$T(f \mid h) = \frac{\exp\left(\Omega\_f(h)\right)}{\sum\_{i=1}^{F} \exp\left(\Omega\_i(h)\right)}\tag{8}$$

Finally, a dropout technique is used, since training a big neural network is too expensive. Dropout freezes neurons at random with a dropout probability (PD) of 0.5. During training phase, those dropped out neurons are not engaged in both a forward and backward pass. During the test phase, all neurons are used but with outputs multiplied by PD of 0.5 (33).

It can be regarded as taking a geometric mean of predictive distributions, generated by exponentially-many smallsize dropout neural networks. **Figure 7A** shows a plain neural network with numbers of neurons at each layer as (2, 4, 8, 10),

and **Figure 7B** shows the corresponding dropout neural network with PD of 0.5, where only (1, 2, 4, 5) neurons remain active at each layer.

### Transfer AlexNet to Alcoholism Identification

First, we needed to modify the structure. The last FCL was revised, since the original FCLs were developed to classify 1,000 categories. Twenty randomly selected classes were listed as: scale, barber chair, lorikeet, miniature poodle, Maltese dog, tabby, beer bottle, desktop computer, bow tie, trombone, crash helmet, cucumber, mailbox, pomegranate, Appenzeller, muzzle, snow leopard, mountain bike, padlock, diamondback. We observed that none of them are related to the brain image. Hence, we could not directly apply AlexNet as the feature extractor. Therefore, fine-tuning was necessary.

Since the length of output neurons in orthodox AlexNet (1000) is not equal to the number of classes in our task (2), we needed to revise the corresponding softmax layer and classification layer. The revision is shown in **Table 4**. In our transfer learning scheme, we used a new randomly-initialized fully connected layer with two neurons, a softmax layer, and a new classification layer with only two classes (alcoholism and non-alcoholism).

Next, we set the training options. Three subtleties were checked before training. First, the whole training epoch should be small for a transfer learning. In this study, we set the number of training epochs to 10. Second, the global learning rate was set to a small value of 10−4 to slow learning down, since the early parts of this neural network were pre-trained. Third, the learning rate of new layers were 10 times that of the transferred layer, since the transferred layers with pre-trained weights/biases and new layers were with random-initialized weights/biases.

Third, we varied the numbers of transferred layers and tested different settings. The AlexNet consists of five conv layers (CL1, CL2, CL3, CL4, and CL5) and three fully-connected layers (FCL6, FL7, FL8). As a result, we tested five different transfer learning settings as shown in **Figure 8** in total, in all experiments. For example, here Setting A means that the layers from the first layer to layer A are transferred directly with learning rate as 10−4 × 1 = 10−4 . The late layers from layer A to the last layer are randomly initialized with a learning rate of 10−4 × 10 = 10−3 .

### Implementation and Measure

We ran the experiment many times. Each time, the trainingvalidation-test division was set at random again. The training procedure stopped when either the algorithm reached maximum

epoch, or the performance of validation decreased over a preset training epoch. We repeatedly tuned the hyperparameters and found those optimal hyper-parameters based on a validation set. After the hyperparameters were fixed, we ran the final model on the test set for 10 runs. The test set confusion matrix across all runs was recorded, and the following five measures were calculated: sensitivity (SEN), specificity (SPC), precision (PRC), accuracy (ACC), and F1 score. Assume TP, TN, FP, and FN stands for true positive, true negative, false positive, and false negative, respectively, all five measures were defined as:

$$\text{SEN} = \frac{TP}{TP + FN} \tag{9}$$

$$\text{SPC} = \frac{TN}{TN + FP} \tag{10}$$

$$\text{PRC} = \frac{\text{TP}}{\text{TP} + \text{FP}} \tag{11}$$

$$ACC = \frac{TP + TN}{TP + TN + FP + FN} \tag{12}$$

F1 considers both the precision and the sensitivity to computer the score (34). That means, the measure of the "F1 score" is

the harmonic mean of the previous two measures: precision and sensitivity.

$$F\_1 = \left(\frac{SEN^{-1} + PRC^{-1}}{2}\right)^{-1} \tag{13}$$

Using simple mathematical knowledge, we can obtain:

$$\begin{array}{l} F\_1 = 2 / \left( \frac{T P + FN}{T P} + \frac{T P + FP}{T P} \right) \\ = 2 / \left( \frac{2 T P + FP + FN}{T P} \right) \\ = \frac{2 \times T P}{2 \times T P + FP + FN} \end{array} \tag{14}$$

Then, the average and standard deviation (SD) of all five measures of 10 runs of the test set were calculated and used for comparison. For ease of understanding, a pseudocode of our experiment is listed below in **Table 5**. The first block is to split the dataset into non-test and test sets. In the second block, the non-test set was split into training and validation randomly. The performance of the retrained AlexNet model was recorded and used to select the optimal transfer learning setting S∗ . In the final block, the performance on the test set via the retrained AlexNet using setting S∗ was recorded and outputted.

### RESULTS

### Data Augmentation Results

**Figure 9** shows the horizontally flipped image. Here, vertical flipping was not carried out because it can be seen as a combination of horizontal flipping with 180-degree rotation. TABLE 5 | Pseudocode of our experiment.

```
for S = [A, B, C, D, E]
```
for i = 1:10 [train(i), valid(i)] = split(NonTest), Model(S, i) = TrainNetwork(AlexNet, train(i), valid(i), Setting = S),

PerfValid(S, i) = Predict(Model(S, i), valid(i)),

end

PerfValid(S) = mean(PerfValid(S, i)),

End

S\* = argmax[Performance(S)],

for i = 1:10

[train(i), valid(i)] = split(NonTest),

Model(S\*, i) = TrainNetwork(AlexNet, train(i), valid(i), Setting = S\*),

PerfTest(S\*, i) = predict(Model(S\*, i), Test),

```
End
```
PerfTest(S\*) = mean(PerfTest(S\*, i)),

Output PerfTest(S\*),

FIGURE 9 | Data augmentation by horizontal flipping. (A) Original image. (B) Flipped image.

**Figure 10** shows the data augmentation results of five different techniques: (a) noise injection; (b) scaling; (c) random translation; (d) image rotation; (e) Gamma correction. Due to the page limit, the data augmentation results on the flipped image are not shown.

### Comparison of Setting of TL

In this experiment, we compared five different TL settings on the validation set. The results of Setting A are shown in **Table 6**, where the last row shows the mean and standard deviation value. The results of Setting E are shown in **Table 7**. Due to page limit, we only show the final results of Setting B, C, and D in **Table 8**.

Here, it can be seen from **Table 8** that Setting E, i.e., replacing the FCL8, achieves the greatest performance among all five settings with respect to all measures. The reason may be (i) we expanded a relatively small dataset to a large training set using data augmentation; and (ii) the dissimilarity of our data and the original 1,000-category dataset. The first fact ensures that retraining avoids overfitting; and the latter fact

TABLE 7 | Ten runs of validation performance of transfer learning using Setting E.

TABLE 8 | Comparison of different setting.

*Bold means the best.*

suggests that it is more practical to put most of the layers initialized with weights from a pretrained model, than freezing those layers. For clarity, we plotted the error bar and show it in **Figure 11**.

### Analysis of Optimized TL Setting

The structure of the optimal transfer learning model (Setting E) is listed in **Table 9**. Compared to the traditional AlexNet model, the weights and biases of FCL8 were reduced from 4,096,000 to 8,192, and from 1,000 to 2, respectively. The main reason is that we only had two categories in our classification task. Thus, the whole weight of the deep neural network reduced slightly from 60,954,656 to 56,866,848.

Nevertheless, we can observe that FCL6 and FCL7 still constitutes too many weights and biases. For example, FCL6 occupied 37,748,736/56,866,848 = 66.38% of the total weights in this optimal model, and FCL7 occupied 16,777,216/56,866,848 = 29.50% of the total weights. Additionally, the FCL subtotal comprised 95.90% of the total weights. This is the main limitation of our method. To solve it, we need to replace the fully connected layers with 1 × 1 conv layers. Another solution is to choose smallsize transfer learning models, such as SqueezeNet, ResNet, GoogleNet, etc.

### Effect of Data Augmentation

This experiment compared the performance of runs with data augmentation against runs without data augmentation (DA). Configuration of transfer learning was set to Setting E. All the other parameters and network structures were the same as the previous experiments. The performance of the 10 runs without using DA are shown in **Table 10**. The results in terms of all measures are equal to or slightly above 95%.

The comparison of using DA against not using DA is shown in **Table 11**. We can discern that DA indeed enhances the classification performance. The reason is that having a large dataset is crucial for good performance. The alcoholism image dataset is commonly of small size, and its size can be augmented to the order of tens of thousands (48,320 in this study). AlexNet can make full use of all its parameters with a big dataset. Without using DA, overfitting is likely to occur in the transferred model.

### Results of Proposed Method

In this experiment, we chose Setting E (replace the final block) as shown in **Figure 8**. Here, the retrained neural network was tested on the test set. The results over all 10 runs on

Name Weights Weights (%) Biases Biases (%) CL1 (Ours) 11\*11\*3\*96 = 34,848 0.06 1\*1\*96 = 96 1.00 CL2 (Ours) 5\*5\*48\*256 = 307,200 0.54 1\*1\*256 = 256 2.68 CL3 (Ours) 3\*3\*256\*384 = 884,736 1.56 1\*1\*384 = 384 4.01 CL4 (Ours) 3\*3\*192\*384 = 663,552 1.17 1\*1\*384 = 384 4.01 CL5 (Ours) 3\*3\*192\*256 = 442,368 0.78 1\*1\*256 = 256 2.68 FCL6 (Ours) 4096\*9216 = 37,748,736 66.38 4096\*1 = 4,096 42.80 FCL7 (Ours) 4096\*4096 = 16,777,216 29.50 4096\*1 = 4,096 42.80 FCL8 (AlexNet) 1000\*4096 = 4,096,000 1000\*1 = 1,000 FCL8 (Ours) 2\*4096 = 8,192 0.01 2\*1 = 2 0.02 CL Subtotal (AlexNet) 2,332,704 1,376 CL Subtotal (Ours) 2,332,704 4.10 1,376 14.38 FCL Subtotal (AlexNet) 58,621,952 9,192 FCL Subtotal (Ours) 54,534,144 95.90 8,194 85.62 Total (AlexNet) 60,954,656 10,568 Total (Ours) 56,866,848 100 9,570 100

TABLE 9 | Learnable layers in optimal transfer learning model.

the test set are listed in **Table 12** with details of sensitivity, specificity, precision, accuracy, and the F1 score of each run. Setting E yielded a sensitivity of 97.44 ± 1.15%, a specificity of 97.41 ± 1.51%, a precision of 97.34 ± 1.49%, an accuracy of 97.42 ± 0.95%, and an F1 score of 97.37% ± 0.97%. Comparing **Table 12** with **Table 7**, we can see that the mean value of test performance is slightly worse than that of the validation performance, but the standard deviation of the test performance is much smaller than that of the validation performance.

### Comparison to Alcoholism Identification Approaches

This proposed transfer learning approach was compared with seven state-of-the-art approaches: PAC-PSO (4), HWT (5), LR (6), CSO (7), WRE (8), SVM-GA (9), and LMCoP (10). The comparison results are shown in **Table 13**. The cognate bar plot is shown in **Figure 12**. We can observe that our AlexNet transfer learning model has more than 3% improvement compared to the next best approach.

The reason is that this proposed model did not need to find features manually; nevertheless, it only used a learned feature from a pretrained model as initialization, and utilized the enhanced training set to fine-tune those learned features. This has two advantages: First, the development is quite fast, which can be reduced to <1 day. Second, the features can be fine-tuned to be more appropriate to this alcoholism classification task than other manually-designated features.

The bioinspired-algorithm may help retraining our AlexNet model. Particle swarm optimization (PSO) (35–37) and other methods will be tested. Cloud computing (38) in particular can be integrated into our method, to help diagnosis of remote patients.

#### TABLE 10 | Ten runs without using data augmentation (Setting E).

#### TABLE 11 | Effect of using data augmentation technique.

TABLE 12 | Ten runs of proposed method on the test set (Setting E).

TABLE 13 | Comparison with state-of-the-art approaches.

### CONCLUSIONS

In this study, we proposed an AlexNet-based transfer learning method and applied it to the alcoholism identification task. This paper may be the first paper using transfer learning in the field of alcoholism identification. The results showed that this proposed approach achieved promising results with a sensitivity of 97.44 ± 1.15%, a specificity of 97.41 ± 1.51%, a precision of 97.34 ± 1.49%, an accuracy of 97.42 ± 0.95%, and an F1 score of 97.37 ± 0.97.

Future studies may include the following points: (i) other deeper transfer learning models, such as ResNet, DenseNet, GoogleNet, SqueezeNet, etc. should be tested; (ii) other data augmentation techniques should be attempted. Currently our dataset is small, so data augmentation may have a distinct effect on improving the performance; (iii) how to set the learning rate factor of each individual layer in the whole neural network, remains a challenge and needs to be solved; (iv) this method is ready to run on a larger dataset and can assist radiologists in their routine screening of brain MR images.

### DATA AVAILABILITY

The datasets for this manuscript are not publicly available because we need approval from our affiliations. Requests to access the datasets should be directed to<EMAIL_ADDRESS>
### AUTHOR CONTRIBUTIONS

S-HW and Y-DZ conceived the study. SX and XC designed the model. CT, JS, and DG analyzed the data. S-HW, XC, and Y-DZ acquired the preprocessed data. SX and CT wrote the draft. S-HW, CT, JS, and Y-DZ interpreted the results. DG provided

### REFERENCES

English revision of this paper. All authors provided critical revision and consent for this submission.

### FUNDING

The authors are grateful for the financial support of the Zhejiang Provincial Natural Science Foundation of China (LY17F010003, Y18F010018), the National key research and development plan (2017YFB1103202), the Open Fund of Guangxi Key Laboratory of Manufacturing System and Advanced Manufacturing Technology (17-259-05-011K), the Natural Science Foundation of China (61602250, U1711263, U1811264), and the Henan Key Research and Development Project<PHONE_NUMBER>29).

**Conflict of Interest Statement:** 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.

Copyright © 2019 Wang, Xie, Chen, Guttery, Tang, Sun and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# A Fully Automatic Framework for Parkinson's Disease Diagnosis by Multi-Modality Images

Jiahang Xu1,2, Fangyang Jiao3 , Yechong Huang1 , Xinzhe Luo1 , Qian Xu4 , Ling Li5 , Xueling Liu6 , Chuantao Zuo5 , Ping Wu5 \* and Xiahai Zhuang1,2 \*

1 School of Data Science, Fudan University, Shanghai, China, 2 Fudan-Xinzailing Joint Research Center for Big Data, Fudan University, Shanghai, China, 3 Department of Nuclear Medicine, Daping Hospital, Army Medical University, Chongqing, China, 4 Department of Nuclear Medicine, North Huashan Hospital, Fudan University, Shanghai, China, 5 PET Center, Huashan Hospital, Fudan University, Shanghai, China, 6 Department of Radiology, Huashan Hospital, Fudan University, Shanghai, China

Background: Parkinson's disease (PD) is a prevalent long-term neurodegenerative disease. Though the criteria of PD diagnosis are relatively well defined, current diagnostic procedures using medical images are labor-intensive and expertise-demanding. Hence, highly integrated automatic diagnostic algorithms are desirable.

### Edited by:

Siyang Zuo, Tianjin University, China

#### Reviewed by:

Delia Cabrera DeBuc, University of Miami, United States Danny J. J. Wang, University of Southern California, United States

#### \*Correspondence:

Xiahai Zhuang<EMAIL_ADDRESS>Ping Wu<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Brain Imaging Methods, a section of the journal Frontiers in Neuroscience

Received: 14 February 2019 Accepted: 05 August 2019 Published: 23 August 2019

#### Citation:

Xu J, Jiao F, Huang Y, Luo X, Xu Q, Li L, Liu X, Zuo C, Wu P and Zhuang X (2019) A Fully Automatic Framework for Parkinson's Disease Diagnosis by Multi-Modality Images. Front. Neurosci. 13:874. doi: 10.3389/fnins.2019.00874 Methods: In this work, we propose an end-to-end multi-modality diagnostic framework, including segmentation, registration, feature extraction and machine learning, to analyze the features of striatum for PD diagnosis. Multi-modality images, including T1-weighted MRI and 11C-CFT PET, are integrated into the proposed framework. The reliability of this method is validated on a dataset with the paired images from 49 PD subjects and 18 Normal (NL) subjects.

Results: We obtained a promising diagnostic accuracy in the PD/NL classification task. Meanwhile, several comparative experiments were conducted to validate the performance of the proposed framework.

Conclusion: We demonstrated that (1) the automatic segmentation provides accurate results for the diagnostic framework, (2) the method combining multi-modality images generates a better prediction accuracy than the method with single-modality PET images, and (3) the volume of the striatum is proved to be irrelevant to PD diagnosis.

#### Keywords: Parkinson's disease, multi-modality, image classification, U-Net, striatum

### INTRODUCTION

Parkinson's disease (PD) is the second-most prevalent long-term neurodegenerative disease characterized by bradykinesia, rigidity and rest tremor (Postuma et al., 2015). At present, PD is responsible for about 346,000 deaths per year and is thus one of the major concerns of neuroscience community (Roth et al., 2018). The diagnosis of PD mainly refers to the Movement Disorder Society Clinical Diagnostic Criteria for Parkinson's disease (MDS-PD Criteria) (Postuma et al., 2015). According to the MDS-PD criteria, the motor symptoms of PD are linked with the loss of dopaminergic neurons, which mainly affects the anatomical regions of the striatum (SARs). Therefore, SARs, which include the caudate nucleus, the putamen and the pallidum, are commonly explored (Strafella et al., 2017).

Functional neuroimaging of the presynaptic dopaminergic system is highlighted according to the MDS-PD criteria (Liu et al., 2018). Positron-emission tomography (PET) is one of the neuroimaging modalities that indicate the regional activity of the tissues. Accordingly, PET tracers are developed to observe the activity of the dopamine transporter (DAT) in early stage of PD, such as 11C-CFT, which is a biomarker of the presynaptic dopaminergic system with high sensitivity (Kazumata et al., 1998; Ilgin et al., 1999; Wang et al., 2013). However, due to the low resolution of the PET images, the anatomical and structural information related to the brain that PET can provide is limited. Therefore, the structural neuroimaging methods, such as T1 weighted magnetic resonance imaging (T1-MRI), are introduced to assist the multi-modality diagnosis of PD (Ibarretxe-Bilbao et al., 2011). Bu et al. (2018) worked on the subtypes of multiple system atrophy (MSA) utilizing T1-MRI and 11C-CFT PET. Huang et al. (2019) combined these two modalities with 18F-FDG PET and analyzed Rapid Eye Movement (REM) Sleep Behavior Disorder research. In both of their studies, T1-MR images were registered to PET images to identify the region of interest (ROIs) in the PET images.

Recently, researchers attempt to improve the accuracy in diagnostic methods with the help of machine learning algorithms, for example, the support vector machine (SVM) has been widely used. Long et al. (2012) used SVM to distinguish early PD patients from NL subjects utilizing resting-state functional MRI, and obtained an accuracy of 86.96%. Haller et al. (2012) used SVM and reached an accuracy of 97% when classifying PD from other atypical forms of Parkinsonism by combining Diffusion Tensor Imaging (DTI) and 123I ioflupane Single-Photon Emission Computed Tomography (SPECT). These works combining multi-modality imaging have proved the reliability of artificial intelligence (AI)-assisted PD diagnosis, while few works are reported including 11C-CFT PET, to the best of our knowledge.

In this work, we proposed an end-to-end multi-modality diagnostic framework for PD combining T1-MR and 11C-CFT PET images.  The resulting segmentation was then used to locate the SARs of the PET images by registration. Finally, features extracted from these SARs were used to diagnose PD. Our main contributions include:

### METHODOLOGY

The proposed framework is shown in **Figure 1**. It contains three major steps: (1) segmentation, (2) registration, and (3) feature extraction and prediction. In the first two steps, MRI-assisted PET segmentation is performed by MRI segmentation and MRI-PET registration, and in the subsequent step, only information of PET images is considered for PD diagnosis.

### Striatum Segmentation via Deep Neural Network

To obtain the fine structure of the brain tissues, a 3D deep neural network, i.e., U-Net (Ronneberger et al., 2015; Wong et al., 2018), is implemented to segment the striatum in the MR images. The obtained segmentation is used as a reference for SAR localization and extraction in the subsequent procedures.

**Figure 2** shows the network architecture for the segmentation, which outputs a mask indicating the segmented labels of the input image. The network further incorporates the idea of deep supervision introduced by Mehta and Sivaswamy (2017) for faster training convergence. Specifically, the network comprises encoding and decoding paths. The encoding path captures contextual information by residual blocks and max-pooling operations at different resolutions, while the decoding path sequentially recovers the spatial resolution and object boundaries. Besides, skip connections between the upsampled feature maps in the decoder and the corresponding feature maps in the encoder are employed for the combination of local and contextual information. Moreover, the deep supervision scheme is adopted to allow more direct backpropagation to the hidden layers for faster convergence. A final 1 × 1 × 1 convolution layer with a softmax function produces the segmentation probabilities. Gaussian blurring and dropout operations are adopted to avoid overfitting. A loss function is defined to handle the relatively small anatomical structures of labels for accurate segmentation, i.e.,

$$L = \left. \boldsymbol{w}\_{\rm D} \boldsymbol{L}\_{\rm Dict} + \boldsymbol{w}\_{\rm C} \boldsymbol{L}\_{\rm Cross},\tag{1}$$

where, wD and wC denote the weights of LDice and LCross, respectively; LDice denotes the Dice-related loss, and LCross denotes the cross-entropy. They are respectively, given by

$$L\_{Dice} = E\_i \left[ \left( -\ln Dice\_i \right)^{\mathcal{V}} \right],\tag{2}$$

with

$$Dice\_i = \frac{2\left(\sum\_{\mathbf{x}} \delta\_{il}\left(\mathbf{x}\right) \cdot p\_i\left(\mathbf{x}\right)\right)}{\sum\_{\mathbf{x}} \left(\delta\_{il}\left(\mathbf{x}\right) + p\_i\left(\mathbf{x}\right)\right)},\tag{3}$$

and

$$L\_{Cross} = E\_{\mathbf{x}} \left[ -\ln p\_l \left( \mathbf{x} \right) \right]. \tag{4}$$

In Eq. (3), δil (x) is the Kronecker delta, which equals to 1 if the segmentation label i(x) equals to the ground-truth label l(x) at the voxel position x, and 0 otherwise; pi (x) is the probability of voxel x being labeled asi. In our implementations, we chose wD = 0.8, wC = 0.2 and γ = 0.3 for the loss function and pretrained the model using an Adam optimizer with a learning rate of 1 × 10−3 for 10 epochs (Kingma and Ba, 2014). Due to the computational limitations, an ROI of MR images with a size of 96 × 96 × 96 voxels was cropped, which contains the whole structure of SARs.

We employ T1-MR brain images from the Alzheimer's Disease Neuroimaging Initiative (ADNI1 ) database for pretraining, for the size of the clinical data used in this work is far from being enough for the U-Net training. Note that the information related to AD or other modality data is not used in this study, namely we solely employ the 1859 brain T1-MR images to assist the U-Net training. The ADNI MR images are segmented by the multi-atlas label propagation with the expectation-maximization (MALP-EM2 ) framework (Ledig et al., 2015). The manual segmentation of the caudate nucleus, the putamen and the pallidum are chosen to be the gold standards in the pretraining stage.

We propose to combine two-modality images for the automatic diagnosis of PD, where T1-MRI provides the morphological information of SARs, and 11C-CFT offers pathological information related to PD. The extraction of the SAR information from the MR images is achieved by the DNN segmentation method, as described in Section "Striatum Segmentation via Deep Neural Network." With this information, one can extract the shape or substructure features from each of anatomical regions. For the combination of the two-modality images, we propose to use image registration, which propagates the anatomical and structural information of SARs in the MRI to the PET.

Combining Two-Modality Images via Image Registration

1 adni.loni.usc.edu

2https://biomedia.doc.ic.ac.uk/software/malp-em/

The registration in the multi-modality diagnostic framework is achieved via the zxhproj3 platform (Zhuang et al., 2011). Firstly, the image with prior label information is registered to the target PET image. The resulting transformation is then used to propagate the prior label information to the PET, which results in the automatic localization of the SARs for the target PET image. Since the MRI and PET images are from the same subject at the same acquisition session, we propose to use a rigid registration. By registration, the caudate nucleus, the putamen and the pallidum, as well as the parieto-occipital regions are labeled.

For comparisons, we propose a single-modality diagnostic framework using solely PET images. To achieve the fully automated diagnosis, we propose to achieve the anatomical information in the PET images via the same registration method used for the multi-modality scheme. In this scenario, the image with prior label information is defined using a pre-labeled PET template, and the registration between the template and the target PET is achieve via an affine registration following a prerigid registration.

### Feature Extraction and Prediction

To extract adequate features from the SARs, the caudate nucleus and the putamen are further divided into three substructures using a k-means algorithm (Ng et al., 2006). After clustering, statistics of image intensity are calculated to represent the feature information in each region, including maximum, minimum, median, 1st and 3rd quantile, and mean of PET intensity. Several studies characterize radioactive uptake by the striatal-to-occipital ratio (SOR), as the parieto-occipital region is widely considered of lacking CFT uptake (Ma et al., 2002; Carbon et al., 2004; Huang et al., 2007). In this work, the SOR, which is defined as (striatumoccipital)/occipital, is calculated with each kind of intensity value. Meanwhile, the volumes of the six anatomical SARs are included into the feature set. In all, 90 features are generated (for a list of specific features, see **Supplementary Table S1**).

After feature extraction, a t-test is performed to analyze the significance of each feature. Setting significance level α = 0.01, features with a p-value less than 0.01 are considered as being statistically significant. Only significant features would be regarded as the arguments of the machine learning models.

Consequently, the SVM classifier is trained to classify the subjects (Haller et al., 2012; Long et al., 2012). Furthermore, to estimate the generalization ability and stability of the method, the leave-n-out cross-validation strategy is employed to evaluate the performance of the models. In addition, we implement the random forest algorithm to calculate the importance of the features (Gray et al., 2013).

### EXPERIMENTS

The following parts in this section are organized as follows. Section "Data Acquisition" describes the data used in this work; Section "Evaluation of Automatic Striatum Segmentation" validates the reliability of the automatic segmentation method; Section "Advantages of Multi-Modality Images" investigates the advantages of combining multi-modality images; and Section "Efficacy of Volume Features" explores the efficacy of the volume features of SARs for the diagnostic of PD.

### Data Acquisition

Data used in this study was collected from the Department of Neurology, Huashan Hospital, Fudan University. It contains paired 11C-CFT PET and T1-MR images of PD patients and healthy participants.  Each MR image was visually inspected to rule out motion artifacts (Bu et al., 2018; Huang et al., 2019). PET images were acquired by a Siemens Biograph 64 PET/CT scanner (Siemens, Munich, Germany) in three-dimensional (3D) mode. A CT transmission scan was first performed for attenuation correction. Static emission data were acquired 60 min after the intravenous injection of 370 MBq of 11C-CFT and lasted for 15 min. All subjects were scanned in a supine position with a dimly lighted and noise-free surrounding (Bu et al., 2018; Huang et al., 2019). The synchronous MRI data were acquired using a T1-weighted 3D inversion recovery spoiled gradient recalled acquisition (IR-SPGR) with the following parameters: TE/TR = 2.8/6.6 ms, inversion time = 400 ms, flip angle = 15◦ , matrix = 256 × 256 × 170, field-of-view = 24 cm, and slice thickness = 1 mm. MR and PET images acquisition for each subject had a time interval of no more than 3 months.

Forty-nine patients with PD and 18 age-matched normal control (NL) subjects were recruited. All subjects were screened and clinically examined by a senior investigator of movement disorders before entering the study and were followed up for at least 1 year. The diagnosis of PD was made referring to the MDS-PD Criteria. The Unified Parkinson's Disease Rating Scale (UPDRS) and Hoehn and Yahr scale (HY) were assessed after the cessation of oral anti-parkinsonian medications (if used) for at least 12 h. The following exclusion criteria were used for the NL subjects' recruitment: (1) being tested positive by the REM Sleep Behavior Disorder Single-Question Screen (Postuma et al., 2012), (2) a history of neurological or psychiatric illness, (3) a prior exposure to neuroleptic agents or drugs, (4) an abnormal neurological examination. The data are summarized in **Table 1**. In this study, gender proportion differences between groups could be ignored, as previous studies have shown no significant difference in DAT bindings between genders (Eshuis et al., 2009). The research was approved by the Ethics Committee of Huashan

TABLE 1 | Summary for the studied dataset.

For gender, the expression means Male/Female, and for age and Unified Parkinson's Disease Rating Scale, the expression means mean ± standard deviation.

3http://www.sdspeople.fudan.edu.cn/zhuangxiahai/0/zxhproj/

Hospital. All subjects or legally responsible relatives signed written informed consent in accordance with the Declaration of Helsinki before the study.

After data acquisition, both sides of the caudate nucleus, the putamen and the pallidum of each MR image were manually labeled by an experienced clinician from the Department of Neurology, Huashan Hospital. To ensure the qualities of the segmentation results, boundaries of these anatomical structures were double-checked by another clinician from the same department.

### Evaluation of Automatic Striatum Segmentation

To test the performance of the segmentation network, three-fold cross-validation was performed. The whole dataset was split into three disjoint parts, and the model was fine-tuned for 5 epochs on the union of every two disjoint subsets. **Table 2** illustrates the average Dice Similarity Coefficient (DSC) of each anatomical region, and **Figure 3** provides a visualization of the segmentation results of five example cases. One can find that the left pallidum (colored goldenrod in **Figure 3**) is worst segmented with the maximal standard deviation while the right putamen (colored olive drab in **Figure 3**) is best segmented with the minimal standard deviation.

**Figure 4** shows the average accuracy (ACC) and the number of wrong predictions with leave-n-out cross-validation of the different segmentation methods, i.e., automatically and manually. Both accuracies reached 100% when n = 1, and the accuracies and the numbers of wrong predictions of the two experiments result in no significant difference in a pairwise t-test (p-value = 0.1017). Furthermore, when training classifiers using features of manually segmented images and testing it using features of automatically segmented images, we still obtained 100% accuracy. All results indicate that the automatic segmentation provides accurate results for the proposed diagnostic framework of PD.

### Advantages of Multi-Modality Images

To evaluate the influence of multi-modality images, the singlemodality method using solely PET images was compared. In the multi-modality scheme, the MR images provides accurate anatomical and structural SAR information of the subject. By contrast, in the single-modality method this information is achieved by registering the PET images to a pre-labeled Automated Anatomical Labeling (AAL) PET template. We conducted the rest of the pipeline in the same way for the two methods.

**Figure 5** shows the results of the comparative experiments with leave-n-out cross-validation. The results demonstrate that with the assistance of MR images, the performance of the multimodality group is better than the single-modality PET group in the PD/NL task. When n = 1, the accuracy of the multi-modality

FIGURE 4 | The results for the leave-n-out cross-validation of the classification with automatic segmentation and manual segmentation. Panel (A) presents the average ACC, and panel (B) presents the average number of wrong predictions. The horizontal axes in the two panels represent subjects numbers of the test set, i.e., the n in the leave-n-out cross-validation.

FIGURE 5 | The results for the leave-n-out cross-validation of the classification by the multi-modality diagnostic method and the single-modality method. Panel (A) presents the average ACC, and panel (B) presents the average number of wrong predictions. The horizontal axes in the two panels represent subjects numbers of the test set, i.e., the n in the leave-n-out cross-validation.

The importance values were calculated by the random forest algorithm.

group reached 100% in the PD/NL task, while the accuracy of the single-modality PET group was 98.51%.

To test the uniformity of the classifiers based on the different groups, we also trained the classifier using features of multi-modality images and tested it using features of singlemodality PET, and the accuracy dropped to 88.05%, with 8 subjects misclassified.

TABLE 4 | Feature importance of groups with manual segmentation results and automatic segmentation results.

The importance values were calculated by the random forest algorithm.

### Efficacy of Volume Features

In the feature extraction step, t-tests were performed to evaluate the significances of all feature, and results indicated that the features of the volume are not statistically significant with α = 0.01 (see **Supplementary Table S1** for more details). To further evaluate the effects of the volume of SARs, we compared the importance of different features based on groups with

fnins-13-00874 August 22, 2019 Time: 17:44 # 6

and without volume, as **Table 3** shows. One can see that the importance values of the two groups are similar. Hence, the effect of volume to the model is negligible. Note that the volume is calculated based on the original MR images without downsampling.

### DISCUSSION AND CONCLUSION

In this work, we proposed a fully automatic framework for PD diagnosis. This method utilized two modalities, i.e., 11C-CFT PET and T1-MR imaging, performed MRI-assisted PET segmentation, selected features and employed SVM to give the predictions. To validate the performance of the framework, we applied the proposed method on the clinical data from Huashan Hospital.

One of the major differences between the proposed method and the traditional methods is that the SARs are located according to the labels of the automatic segmentation by U-Net. To evaluate the performance of the U-Net, we calculated the DSCs between automatic and manual segmentation. In addition, we compared the proposed pipeline, whose SARs were located according to the automatic segmentation, to the method whose SARs were manually segmented. The leave-nout experiment shows the two methods performed comparably,

FIGURE 6 | The comparison of gold standard and wrongly placed SARs of the wrongly predicted subject. Panel (A) shows the segmentation result in gold standard, and panel (B) shows the segmentation in the wrongly predicted subject. Images are in sagittal plane and have the same cursor position.

fnins-13-00874 August 22, 2019 Time: 17:44 # 7

indicating that the automatic segmentation could provide accurate results for the proposed diagnostic framework of PD. Further investigation of the feature importance of the two groups is illustrated in **Table 4**. It indicates that the minimum has lower importance than the first five features. Given that the striatum region has a higher uptake value compared with its adjacent areas, voxels with minimal intensity value are more likely to appear on the edge of the SARs. Therefore, the inaccurate delineation of the anatomical boundary as a potential result of the automatic segmentation could not cause a significant decline in the performance of the overall diagnostic framework.

An alternative way to locate SARs for subsequent feature extraction is to apply a pre-labeled PET template by registration. In Section "Advantages of Multi-Modality Images," AAL PET template was used as the PET template, and was registered to PET images for the localization of SARs. Experiments show that the diagnostic capability of this singlemodality PET group is worse than the proposed multi-modality framework. Though the single-modality PET approach gives a favorable prediction, the multi-modality approach performs better. This is because the localization of the SARs occupies a significant place in the diagnostic framework, and the additional structural information from MR images can better locate SARs. **Figure 6** demonstrates that the single-modality PET approach might be affected by the erroneous delineation of the SARs. The error could be attributed to the ignored inter-subject variations in brain structures when defining SARs from a PET template.

To test the uniformity of the classifiers based on different segmentation approaches, we trained classifiers using features of manually segmented multi-modality images, and tested it using features of other methods. When testing with features of multi-modality automatic segmentation method, we still obtained 100% accuracy, indicating that the features of manual and automatic segmentation are highly consistent. However, testing with the singlemodality method resulted in an accuracy of 88.05%. The lower accuracy might be explained by the lack of adequate extracted features due to the falsely located SARs. Hence, compared with the multi-modality group, single-modality PET group naturally needs more feature engineering and better-designed algorithms.

In the feature extraction, features of volume were rejected according to the t-test. This could be the reason why the volume of SARs does not change significantly with the progression of PD, as concluded from the literature (Ibarretxe-Bilbao et al., 2011). **Figure 7** shows the heatmaps of feature distribution on the SARs, displaying the influence of each subregion for the classification in the PD/NL task. The difference of influence is expressed by the color scale. One can find that the most relevant region influencing the separation of PD/NL is localized in the middle and rear of the putamen, then the pallidum, and the caudate nucleus reveals the least significance on this task.

Several future studies could be completed explored based on the current pipeline. Firstly, the classifiers can be trained with Parkinsonian disorders (PDS) dataset to classify PD and atypical PDS, such as MSA and Progressive Supranuclear Palsy (PSP), which has important clinical values but is with great challenges. Secondly, this framework only contains medical imaging information currently, while other aspects of information, such as age, gender, motor ratings and other biomarkers are not included, which may further improve the diagnostic accuracy. Future research could be undertaken to incorporate additional multimodal data for better disease prediction. Finally, the sample size of subjects in this work is relatively small, and a bigger dataset is expected to validate our experiment results and improve the performance of the framework.

To conclude, we proposed a fully automatic framework combining the two modalities for PD diagnosis. This framework obtained a promising diagnostic accuracy in the PD/NL task. In addition, this work also emphasized the high value of the 11C-CFT PET in the PD diagnosis.

### DATA AVAILABILITY

The datasets generated for this study are available on request to the corresponding author.

## ETHICS STATEMENT

All data were collected from Department of Neurology, Huashan Hospital, Fudan University, and the study was approved by the Ethics Committee of Huashan Hospital. All subjects or a legally responsible relative were given written informed consent before the study.

## AUTHOR CONTRIBUTIONS

XZ is the principle investigator of this work, designed the experiments, and supervised and revised the manuscript. PW co-investigated the research and revised the manuscript. CZ co-investigated the research. JX led the implementation and experiments and wrote the manuscript. YH coled the work and wrote the manuscript. XiL provided support to the work of coding, experiments, and manuscript writing. FJ, QX, and LL collected the data. XuL collected the data and segmented the striatum of the subjects.

### FUNDING

This work was funded by the Science and Technology Commission of Shanghai Municipality Grant (17JC1401600), and the National Natural Science Foundation of China (NSFC) Grant (61971142). Collection of the PD data was supported by the National Natural Science Foundation of China (81671239 and<PHONE_NUMBER>3).

### ACKNOWLEDGMENTS

fnins-13-00874 August 22, 2019 Time: 17:44 # 9

We thank Lei Li and Yuncheng Zhou, whose dedicated comments helped us to modify the organization of the manuscript and improve its quality.

### REFERENCES

### SUPPLEMENTARY MATERIAL

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

**Conflict of Interest Statement:** 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.

Copyright © 2019 Xu, Jiao, Huang, Luo, Xu, Li, Liu, Zuo, Wu and Zhuang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Brain White Matter Hyperintensity Lesion Characterization in T2 Fluid-Attenuated Inversion Recovery Magnetic Resonance Images: Shape, Texture, and Potential Growth

Chih-Ying Gwo1 \* † , David C. Zhu2† and Rong Zhang3†

*1 Department of Information Management, Chien Hsin University of Science and Technology, Zhongli District, Taiwan, 2 Department of Radiology and Psychology, and Cognitive Imaging Research Center, Michigan State University, East Lansing, MI, United States, 3 Department of Neurology and Neurotherapeutics, Department of Internal Medicine, University of Texas Southwestern Medical Center and Institute for Exercise and Environmental Medicine, Texas Health Presbyterian Hospital Dallas, Dallas, TX, United States*

### Edited by:

*Nianyin Zeng, Xiamen University, China*

### Reviewed by:

*Xia-an Bi, Hunan Normal University, China Yuan Yang, Northwestern University, United States*

> \*Correspondence: *Chih-Ying Gwo<EMAIL_ADDRESS>
*†These authors have contributed equally to this work as co-authors*

#### Specialty section:

*This article was submitted to Brain Imaging Methods, a section of the journal Frontiers in Neuroscience*

Received: *31 January 2019* Accepted: *27 March 2019* Published: *16 April 2019*

#### Citation:

*Gwo C-Y, Zhu DC and Zhang R (2019) Brain White Matter Hyperintensity Lesion Characterization in T*2 *Fluid-Attenuated Inversion Recovery Magnetic Resonance Images: Shape, Texture, and Potential Growth. Front. Neurosci. 13:353. doi: 10.3389/fnins.2019.00353* Prior methods in characterizing age-related white matter hyperintensity (WMH) lesions on T2 fluid-attenuated inversion recovery (FLAIR) magnetic resonance images (MRI) have mainly been limited to understanding the sizes of, and occasionally the locations of WMH lesions. Systematic morphological characterization has been missing. In this work, we proposed innovative methods to fill this knowledge gap. We developed an innovative and proof-of-concept method to characterize and quantify the shape (based on Zernike transformation) and texture (based on fuzzy logic) of WMH lesions. We have also developed a multi-dimension feature vector approach to cluster WMH lesions into distinctive groups based on their shape and then texture features. We then developed an approach to calculate the potential growth index (PGI) of WMH lesions based on the image intensity distributions at the edge of the WMH lesions using a region-growing algorithm. High-quality T2 FLAIR images containing clearly identifiable WMH lesions with various sizes from six cognitively normal older adults were used in our method development Analyses of Variance (ANOVAs) showed significant differences in PGI among WMH group clusters in terms of either the shape (*P* = 1.06 × 10−2 ) or the texture (*P* < 1 × 10−20) features. In conclusion, we propose a systematic framework on which the shape and texture features of WMH lesions can be quantified and may be used to predict lesion growth in older adults.

Keywords: brain T2 FLAIR hyperintensity, shape, texture, potential growth, morphology

## INTRODUCTION

The presence of white matter hyperintensities (WMH) on T2 fluid-attenuated inversion recovery (FLAIR) magnetic resonance images (MRI) is common 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

**40**

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 cerebrovascular and neurodegenerative disease diagnosis and to assess treatment effects (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 most commonly used methods for WMH quantification in brain aging, vascular, and AD type of dementia 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).

WMH shape is a basic morphological feature which can be derived from T2 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). Current shape classification methods include mainly the following: (1) one-dimensional function shape representation (Kauppinen et al., 1995; Yadav et al., 2007; Zhang and Lu), (2) polygonal approximation (ShuiHua and ShuangYuan), (3) spatial interrelation feature (Sebastian et al., 2004; Guru and Nagendraswamy, 2007; Bauckhage and Tsotsos), (4) moments (Mukundan, 2004; Celebi and Aslandogan; Taubin and Cooper), (5) scale-space methods (Zhang and Lu, 2003; Kpalma and Ronsin, 2006), and (6) shape transform domains (Chen and Bui, 1999; Zhang and Lu). These methods for shape classification may be suitable for specific applications in various fields but have major limitations for shape characterization of brain lesions. For example, method (1) is highly sensitive to noise, and inaccurate boundary definition can cause large errors; method (2) can only represent the object appearance but not all shape features; method (3) may be used to describe the general appearance of an object, but is limited by the orientation and size of the object; method (4) contains redundant information in the image feature vectors and thus unique images cannot be reconstructed back; method (5) is limited to the shapes which have shallow concavities/convexities; and method (6) requires the definition of the shape contour starting point derived from other methods. In order to characterize the complex shapes of brain lesions, a technique needs to be invariant to the orientation of a lesion, be resistant to image noise and be able to define a one-to-one relationship between feature vector and shape. In this regard, Zernike transformation can satisfy these criteria (Khotanzad and Hong, 1990). Similar to Fourier analysis, shape features of an object captured on MRI can be represented by the coefficients of shape function of the Zernike polynomial expansion (i.e., Zernike transform), referred to as Zernike moments (ZMs) (Zernike, 1934). In this study, we applied Zernike transformation to extract WMH shape features for pattern recognition and classification in cognitively normal older adults.

Image texture is another morphological feature which can be categorized through modeling (Chen et al., 1989), structure (Chow and Rahman, 2007), transformation (Tsai and Hsiao, 2001), and statistics based methods (Haralick et al., 1973; Iivarinen et al., 1996). Model-based and structure-based methods work best for repeating texture patterns but are not suitable for irregular texture patterns such as those in brain lesion images. The transformation based method works best in identifying subregions with known characteristics, but does not work well on unknown and potentially complicated patterns such as those in brain lesion images. Statistical-based methods describe the texture in the distribution of and relationships between gray-level values in an image. These statistics-based methods can normally describe objects better than the structure and transformation based methods because they are invariant to the orientation, the size of an object, and also robust to the noise inside the object (Castellano et al., 2004). Since WMH lesions often have various sizes, orientations, and locations, and are manifested across multiple image slices, a statistics-based method is likely to be the best choice to accommodate these complexities. Therefore, we adopted a statistical method based on fuzzy logic to construct the image intensity histogram of WMH lesions for texture feature extraction.

Finally, we have thought that as a potential imaging biomarker of brain aging and CSVD, the size, shape, and image texture of WMH lesion may change with time (Sachdev et al., 2007; Godin et al., 2011) which may reflect the progression of the underlying pathological process. In this regard, recent studies have shown that the immediate surrounding areas of clearly defined WMH lesions may be at risk for further tissue damage and conversion to lesions (Maillard et al., 2014; Promjunyakul et al., 2016). These areas are classified as WMH penumbras (Maillard et al., 2014). To characterize WMH lesions as well as their penumbras, we developed a seed-based region-growing algorithm to characterize WMH boundaries to explore the potential growth of WMH lesions. We defined this specific WMH boundary characteristic as potential growth index (PGI). To explore whether the shape and texture characterization techniques can potentially be used to predict lesion growth, we assessed whether different shape and texture patterns are related to PGI.

### METHODS AND RESULTS

### MRI Acquisition

Full-brain 2D T2 FLAIR images were collected on a Philips Achieva 3T scanner (Philips Healthcare, Best, the Netherlands) with the following parameters: axial, time of echo (TE) = 125 ms, time of repetition (TR) = 11 s, time of inversion (TI) = 2,800 ms, field of view (FOV) = 23 cm × 23 cm, slice thickness = 2.5 mm, number of slices = 64 with no gaps, acquisition matrix size = 352 × 212, and reconstructed matrix size = 512 × 512. All subjects signed informed consent approved by the Institutional Review Boards of the UT Southwestern Medical Center and Texas Health Presbyterian Hospital of Dallas. Six T2 FLAIR brain image datasets (two male, four female, 75 ± 4 years old and normal cognition), which contained clearly identifiable white matter hyperintensity (WMH) lesions with various sizes, were selected from an healthy aging study we published previously (Tarumi et al., 2014).

## T2 FLAIR Image Segmentation

T2 FLAIR WMH regions were segmented on each 2D image 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 T2 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 (Maldjian et al., 2013). The fitness of a new T2 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.

### Lesion Size Distribution

WMH binary masks generated from 2D T2 FLAIR images (**Figure 1**) were used to obtain WMH size distribution. To minimize artifacts, only those masks with more than 10 connected WMH voxels (voxel size: 0.45 mm × 0.45 mm) on an image were considered probable lesions and were used for further characterization, which resulted in a total of 993 WMH lesions. Fitting each of these lesions within a square, these lesions had a size range of from 6 × 6 to 176 × 176 voxels. The lesion size distributions of six subjects are shown in **Figure 2**. Of note, more than 93% of these lesions are ≤ 60 × 60, and only about 1.5% are ≥ 120 × 120 voxels.

### WMH Shape Feature Extraction and Classification in 2D WMH Shape Feature Extraction Using Zernike Transformation

Zernike transformation has been used extensively in imaging shape feature extraction and pattern recognition (Papakostas et al., 2007; Wee and Paramesran, 2007). The coefficients of Zernike polynomial expansion of an object are referred to as Zernike moments (ZMs) which are used to represent the shape features of analyzed objects.  The lower-order ZMs describe global contour and gross shape features, and the higher-order ZMs describe regional and fine topological details of a shape (Gwo and Wei, 2016). Of note, the magnitudes of ZMs are not only rotational invariant but also robust to small perturbations on the contour of a shape image (Teh and Chin, 1988).

For a 2D image object (a WMH lesion image segmented from a T2 FLAIR image in this work), using polar coordinates, the complex Zernike moments of order n with repetition m can be represented as the inner product of a shape function f (r, θ) with the basis function of Zernike polynomials, Vnm (r, θ), specifically as

$$Z\_{nm} = \frac{n+1}{\pi} \int\_0^{2\pi} \int\_0^1 f\left(r, \theta\right) V\_{nm}^\*\left(r, \theta\right) r dr d\theta, \qquad |r| \le 1,\tag{1}$$

where V ∗ nm (r, θ) denotes the complex conjugation of Vnm (r, θ). The basis function of Zernike polynomial is given by

$$V\_{nm}\left(r,\theta\right) = R\_{nm}\left(r\right)e^{im\theta}, \quad i = \sqrt{-1} \tag{2}$$

where the radial polynomial, Rnm (r), is defined as follows:

$$R\_{nm}(r) = \sum\_{k}^{\frac{n-|m|}{2}} (-1)^k \frac{(n-k)!}{k! \left(\frac{n+|m|}{2} - k\right)! \left(\frac{n-|m|}{2} - k\right)!} r^{n-2k} \tag{3}$$

where 0 ≤ |m| ≤ n, n − |m| is an even integer, and n ≥ 0.

Since the shape features represented by ZMs at orders higher than six are usually too small (small ZM magnitude) to be detected reliably by human eyes (Charman, 2005), the maximum Zernike transformation order was set to five in this study (**Figure 3**). In Zernike transformation, Vn,+m (r, θ) = Vn,−m (r, θ) , and Zn,+m  = Zn,−m . The number of distinctive ZM magnitudes for an expansion up to order n is computed as follows:

$$\begin{cases} \left(\frac{n+2}{2}\right)^2 & \text{if ordern is even} \\ \frac{(n+3)(n+1)}{4} & \text{if order n is odd} \end{cases} \tag{4}$$

Shape feature extraction procedures based on Zernike transformation are illustrated in **Figure 4**. In this illustration, we chose three WMH masks, two with a similar shape but different sizes, and one with both different shape and size. To simplify computation complexity, these image masks with different sizes were first scaled to the same size of 60 × 60 voxels so that the ZM magnitudes can be compared on a same scale. Each Znmwas calculated using Equation (1–3). The calculation resulted with 21 ZM complex coefficients with maximum order n = 5. Based on the magnitudes of the ZM coefficient, only 12 coefficients were needed to extract shape features since the WMH shape (mask) generated with tissue segmentation is rotational invariant. As shown in the right column of **Figure 4**, the two lesion images (a) and (b) with a similar shape have similar ZMs magnitudes at all 12 coefficients. On the contrary, the ZM magnitudes of WMH lesions with a different shape are different from the other two.

### WMH Shape Classification

We then classified the lesion images to different clusters (groups) based on the similarity on shape features. Current common clustering algorithms, such as the K-means clustering algorithm,

lesion prediction algorithm (LPA) showing in red, and a WMH binary mask after tissue segmentation, which was used in shape feature extraction.

requires data-specific a priori selection on the number of clusters (Zhao, 2012). For instance, if the number of clusters is too small, the WMH lesion images with noticeable different shapes may be grouped inappropriately into a same cluster. On the other hand, if the number of clusters is too large, lesion images with trivial differences may be assigned into different clusters, confounding potential clinical significance. Finding the appropriate number of clusters using model simulation is one way to resolve this dilemma (Zhao, 2012). However, this procedure has to be carried out for all choices of shape feature dimensions. To simplify the procedures, the estimation of cluster characteristic indices based on sum of within-cluster dispersions [Wk in Equation (5)] or its variants were proposed (Ball and Hall, 1965; Calinski and Harabasz, 1974; Xu, 1997; Tibshirani et al., 2001). For a better understanding of the influences of the classifiable number of clusters and the feature dimensions derived from the Zernike transform on Wk , we plotted Wk as a function of cluster numbers and feature dimensions which are equivalent to the numbers of the distinctive magnitudes of the ZMs (**Figure 5**). The ZMs for WMH shapes from one to 10 orders were calculated to generate 2 to 36 dimensional feature vectors (Equation 4). Euclidean distance was then calculated to assess the similarity between the feature vectors. The K-means clustering algorithm was applied for grouping purpose. Wk was calculated based on the 2 to 20 cluster groups at each feature dimensions. Wk , in general, decreases with the increase of the number of clusters but increases with the increase of the number of feature dimensions (**Figure 5**). For a specified feature dimension, a better grouping result is likely achieved at a lower Wk value by finding a local minimum. However, in some feature dimensions, a local minimum cannot be found even after Wk decreases to nearly constant. For example, when two is selected as the feature dimension, the Wk value remains small even at small number of clusters because two-dimensional feature vectors only represents

FIGURE 3 | (A) The 21 basis functions of Zernike polynomials, *Vnm* (*r*, θ), with order *n* ≤ 5, are illustrated. The polynomials have a radial range of [−1, 1] (|*Vnm* (*r*, θ)| ≤ 1), shown by the color bar on the left column; (B) 12 distinctive magnitude images, which are rotational invariance, are shown, corresponding to the polynomials in (A).

The images are normalized to the same size of 60 × 60 voxels shown in the middle column. The magnitudes of 12 Zernike moments coefficients based on ZM orders ≤5 of Zernike polynomial expansion are shown in the right column for comparison.

the gross global contour and thus cannot differentiate shapes with enough details. Therefore, selecting an appropriate feature dimension is also crucial and will be discussed later. Nevertheless, once a feature dimension is selected [which was selected to be 12 in this study based on our exploration of data features (**Figures 4**, **5**)], the optimal number of clusters can be determined based on the estimation of cluster characteristics discussed below (Desgraupes, 2013).

The distance between two points in a share feature vector space can be calculated based on Euclidian distance. The overall distance of all points in a cluster to their mean indicates the compactness of a cluster, or within-cluster dispersion. To determine the optimal number of shape clusters for WMH shape classification, we then employed a gap statistics method proposed by Tibshirani et al. (2001). This method estimates the optimal number of clusters by comparing the logarithm of the sum of within-cluster dispersions of a set of clusters to that from the reference datasets created through sampling uniformly at random from the original dataset. The sum of all within-cluster dispersions decreases gradually with the increase of number of clusters but becomes nearly constant at some points as demonstrated in **Figure 5**. This is so called "elbow" phenomenon, which has been used to find the optimal number of clusters (Tibshirani et al., 2001). The algorithm used to estimate the optimal number of WMH shape clusters based on the gap statistic is presented below:

1. Group the shape vectors by varying the number of shape clusters from k = 1, 2, . . . , N (pre-defined as the maximum to evaluate), and compute the sum of the within-cluster dispersion Wk for each choice k.

$$W\_k = \sum\_{r=1}^k \sum\_{\mathbf{x}\_i \in \mathcal{C}\_r} (\mathbf{x}\_i - \bar{\mathbf{x}}\_r)^2 \tag{5}$$

where xi is a data point, Cr denotes cluster r, and x¯r is the vector mean of Cr .

Generate reference datasets (total number = B) by sampling uniformly at random from the original dataset within its distribution ranges of all dimensions. Although a better statistically randomness is likely achieved with a large B, the choice of B is bounded by computation demand. For each reference dataset b, we can generate k clusters, and we can calculate the sum of the within-cluster dispersion Wkb for each k based on Equation (5) above, where b = 1, 2, . . . , B; k = 1, 2, . . . , N. The gap statistics for each k is calculated as below:

$$\log p\left(k\right) = \frac{1}{B} \sum\_{b=1}^{B} \log \left(W\_{kb}\right) - \log \left(W\_k\right) \tag{6}$$

2. letl = (1/B) P b log (Wkb), compute the standard deviation

$$sd\_k = \left[\frac{1}{B} \sum\_{b=1}^{B} \left(\log\left(W\_{kb}\right) - l\right)^2\right]^{1/2} \tag{7}$$

Let sk = sdk √ (1 + 1/B). Choose the optimal number of shape clusters kopt by Equation (8)

$$k\_{opt} = \text{smallest} \quad k \quad \text{such} \quad \text{that } \text{Gap}\left(k\right) \ge \text{Gap}\left(k+1\right) - s\_{k+1} \text{(8)}$$

In the gap statistic procedure above, N is a pre-selected number of shape clusters such that kopt can be determined in the range of [1, N]. B is selected such that the value of sdk converges. In this study, N and B were set to 20 and 10, respectively.

For the WMH shape datasets in this study, based on the previous discussion and **Figure 5**, the "elbow" phenomenon was sufficiently noticeable when the feature dimension was set to 12. At this feature dimension, only ZM magnitudes corresponding to ZM orders of n = 0 to 5 were used in clustering [cf. Equation (4)]. The Gap values were calculated and displayed in **Figure 6**; the optimal number of shape clusters was selected to be six according to Equation (8).

**Figure 7** shows the WMH shape classification results using the K-means algorithm based on the cluster number of six and feature dimension of 12. Unique shape difference can be visualized between the six clusters.

### WMH Texture Feature and Classification Texture Feature Extraction

Image texture characterizes the voxel signal intensity distribution patterns in a WMH region. Statistics-based methods quantify the distribution and relationships of voxel signal values in an image region. These methods often provide better discrimination indexes than structure and spectral transformation based methods (Castellano et al., 2004).

WMH lesions often have various sizes, orientations and locations, and manifest across multiple image slices. In this study, the distributions of WMH lesion size measured in the number of voxels are presented in **Figure 8**.

Most lesions are small, with 48.64% of lesions ≤40 voxels (8 mm2 on a slice). Therefore, a robust texture analysis method needs to satisfy three requirements: (1) Texture feature should be independent of lesion orientation and location; (2) texture feature should be able to quantify small lesions, and (3) texture characterization needs to go beyond a single image slice. A statisticsbased method for WMH texture feature extraction is described next.

Since WMH lesions manifests across multiple slices, we used the "WMH3D" term to emphasize the 3D perspective. Specifically, if a WMH lesion image in a slice connects directly either above, below, or diagonally to another WMH lesion image in an adjacent slice, we treat these lesion images belonging to the same lesion, called it a "WMH3D" for texture analysis. This treatment also reduces the chance of false positive in lesion identification. To characterize the voxel signal intensity distribution, potential voxel spike noise, which is often seen in images, needs to remove first. This can be accomplished by setting the voxel intensities within the boundaries of above or below three standard deviations of the mean values. To reduce the slice variation in signal intensity, a min-max normalization was applied to a WMH3D to normalize its voxel intensity based on the equation,

$$s\left(\mathbf{x},\boldsymbol{y},\boldsymbol{z}\right) = \frac{f\left(\mathbf{x},\boldsymbol{y},\boldsymbol{z}\right) - \operatorname{gMin}}{\operatorname{gMax} - \operatorname{gMin}}\tag{9}$$