基于三元泛函主成分分析和聚类分析的MRI图像分类研究

张嘉茗，林 楠，张梦翰，张 亮，李淑元，熊墨淼，王笑峰，金 力

(1. 复旦大学 生命科学学院，上海 200438；2. 德克萨斯大学 公共卫生学院，休斯顿TX 77030；. 第二军医大学附属长海医院 心血管内科，上海 200433)

基于三元泛函主成分分析和聚类分析的MRI图像分类研究

张嘉茗1，林 楠2，张梦翰1，张 亮3，李淑元1，熊墨淼2，王笑峰1，金 力1

(1. 复旦大学 生命科学学院，上海 200438；2. 德克萨斯大学 公共卫生学院，休斯顿TX 77030；. 第二军医大学附属长海医院 心血管内科，上海 200433)

一元泛函主成分分析(FPCA)已经在fMRI上成功进行了应用，但是目前很少有研究运用多元FPCA对MRI进行探索.本研究将一元FPCA推广到三元并应用于MRI的图像特征提取，并对提取的特征进行了后续研究，提出了一整套MRI病理及正常图像的分析方法.该方法的主要流程是先对MRI图像进行预处理(图像配准和图像分割)，得到脑脊液图像，然后运用三元FPCA对脑脊液进行特征提取，再对提取的特征进行选择，随后利用k- means聚类算法对特征进行聚类，来判断图像所属的类别(正常或异常)，从而达到颅脑MRI图像病变筛查的目的.将该方法应用于颅脑MRI快速自旋回波T2加权像中，结果发现，相比于传统PCA，三元FPCA展现出更好的特征提取能力，可以有效提高图像分类的准确率.

泛函主成分分析； 主成分分析； 颅脑MRI； 聚类分析

随着CT(Computed Tomography, 计算机断层扫描技术)、MRI(Magnetic Resonance Image, 磁共振成像)等影像学技术的发展，医学影像的分辨率在极大提高的同时，MRI等医学影像的数据量也在呈指数级增长.然而，在疾病筛查工作中，图像阅读和诊断的较大工作量是放射科工作人员难以承受的.此外，大量的图像数据和临床信息还会导致MRI图像中存在较多的结构噪音和复杂的疾病状态表现，这极大影响了医生对病症的诊断效率和准确度.因此，如何从大量的MRI图像数据中提取有效的特征，辅助医生进行关于结构或组织病理的诊断，减轻诊断医师的工作负荷并提高诊断的准确性，是目前医学图像分析领域中一个亟待解决的问题.

目前关于MRI图像特征提取的主要方法有主成分分析(Principal Component Analysis, PCA)[1]、独立成分分析(Independent Component Analysis, ICA)[2]等.PCA能从数据的协方差阵中提取不相关的主要活动信号；而ICA不仅能抑制高斯白噪声，还能分解出独立的非高斯信号.但是，ICA方法获得的各个独立信号的幅值缺乏实际意义，导致该方法无法对多个独立信号间的关系进行后续的量化分析.

通常，MRI体数据多为规则结构化体数据集，该数据的特点是在空间中存在明确的拓扑相邻关系，因而可以定义为一个三维空间上的采样函数.但是ICA方法仅将MRI体数据作为一种离散数据进行分析，未能考虑到不同切层图像之间的相邻性以及图像信号随空间变化的潜在规律.这使得上述方法无法有效扩展到MRI三维体数据的分析工作中.针对这类局限性，Ramsay和Dalzell[3]正式提出了泛函数据分析(Functional Data Analysis, FDA)的概念，并使用泛函主成分分析(Functional Principal Component Analysis, FPCA)方法对加拿大温度与降水量的关系进行了实证研究.与传统数据分析相比，FDA将连续变化的曲线作为研究的单元，能有效解决数据采集时间间隔不同或者数据缺失的问题.从动态的角度考虑问题，FDA的优势是没有传统统计方法的参数限制，从而使其具有更强的适用性.另外，FDA将大量数据描述为函数形式，这有助于在处理海量数据的过程中对单一曲线求导或者对多个曲线进行降维，从而为分析曲线之间的差异和曲线内部的动态变化模式提供更多的信息.

近年来，泛函数据分析方法已经开始受到学者的关注，并在功能性磁共振成像(functional Magnetic Resonance Imaging, fMRI)领域展开了应用.Viviani等[4]将FPCA应用于单个被试的fMRI，通过对比FPCA和传统PCA的结果，发现FPCA可以更有效地提取fMRI中血氧水平的变化特征；Long等[5]运用FPCA估计多个被试fMRI中的非平稳噪声.然而上述研究运用的FPCA仅仅局限于一元情形，针对多元FPCA在MRI图像中的应用，目前极少有文献进行讨论.针对MRI三维数据的特点，本文将一元FPCA推广到三元FPCA，并给予相应的数学证明.

本文以颅脑MRI病变筛查为研究目标，采用三元FPCA的图像特征提取方法，对MRI病理及正常图像的分类进行了研究.主要研究步骤如下： 1) 对MRI图像进行预处理(图像配准和图像分割)，得到脑脊液图像；2) 运用三元FPCA对脑脊液进行特征提取，再对提取的特征进行特征选择；3) 利用k- means聚类算法对已选择的特征进行聚类，来判断图像所属的类别(正常或异常)，从而达到颅脑MRI图像病变筛查的目的.为了验证方法的可靠性和准确性，我们将实验结果与传统PCA方法进行对比，结果证明三元FPCA展现出更好的特征提取能力，可以有效提高图像分类的准确率.

1 三元泛函主成分分析原理

一般来说，多元统计分析处理的是每一个观测样本的多个变量在同一时期或时点上呈现出的数据.为了对多变量之间存在的信息冗余进行约减，通常采用主成分分析方法将原来较多的变量转化为少数几个综合变量.

泛函数据分析的基本思想是将每个观测对象带有泛函性质的数据看作一个整体，而不仅仅是个体观测值的一个序列.泛函分析把传统的有限维空间中的数据扩展到了无限维空间，以一个全新的视角对数据进行分析，然而泛函数据内在的无限维特征也给相关问题的统计推断带来了挑战.对此，可以将多元主成分分析的方法引入到泛函数据中，称为泛函主成分分析(FPCA).根据观测对象测度变量的个数，泛函主成分分析可分为一元FPCA和多元FPCA.本文对多元FPCA进行了研究，在一元FPCA的基础上，推导出三元FPCA的计算方法及其求解过程.下面以三维MRI体数据为例，对三元FPCA进行介绍.

1.1三元FPCA

MRI体数据作为三维离散图像，它可以表示成三维函数f(s,t,u),其中：

f在三维空间坐标点(s,t,u)的值为MR成像时表现的灰度值.

假设观测数据是已经过零均值化处理的MRI体数据集xi(s,t,u),i=1,2,…，N，其中(s,t,u)表示三维空间中的位置坐标.设Ω是空间区域，(s,t,u)∈Ω,则三元泛函变量的主成分f定义如下：

(1)

其中β(s,t,u)是权重函数.

进而，主成分f的方差为

(2)

其中R(s1,t1,u1,s2,t2,u2)=cov(x(s1,t1,u1),x(s2,t2,u2)),是x(s,t,u)的协方差函数.

为了避免主成分的权重函数β(s,t,u)过度粗糙，导致主成分分析结果的可解释性降低，本文在泛函主成分分析过程中引入了粗糙惩罚项来控制函数的粗糙程度[6]，从而得到较平滑的权重函数.

(3)

第k个主成分满足：

(4)

式(4)中μ为平滑参数，约束条件等号左边第二项就是粗糙惩罚项.μ越大，得到的权重函数越平滑，相应的主成分函数的方差受β(s,t,u)粗糙度影响也越大，反之，得到的权重函数波动越激烈.当μ→0时，问题变为未施加惩罚的泛函主成分分析问题.这里平滑参数μ的选择采用交叉验证的方法[7,8].

应用拉格朗日乘数法，上述条件极值问题可转换为如下最大化问题：

(5)

其中ρ为参数.

定义泛函如下：

(6)

则对式(5)的求解转变为对式(6)泛函极值的求解，而泛函极值可以通过变分法[9]求得.

设

(7)

其中δβ(s,t,u)为β(s,t,u)的变分，ε为任意小参数，h(s,t,u)为任意函数.则式(6)的一阶变分为

(8)

由极值条件

(9)

得

(10)

由h(s,t,u)的任意性，得

(11)

其中λ是特征值，即三元泛函主成分求解问题可以转化为求解式(11)的积分方程问题.

1.2三元FPCA的求解

由上述可知，平滑的泛函主成分的求解可以转化为式(11)特征值和特征函数的求解问题.在实际计算过程中，常用的求解方法是对函数用适当的基展开.本文采用的是傅里叶标准正交基.

假设经过标准化处理后的MRI体数据样本为xi(s,t,u),i=1,2,…，N，{φj(t)}是傅里叶基函数，对于每个j,定义ω2j-1=ω2j=2πj.选择一组基函数{φk(t)}(k=1,2,…，K)，对样本数据进行傅里叶展开有

(12)

(13)

其中⊗表示矩阵的Kronecker积.

类似地有

(14)

(15)

其中b=[b111,…，b11K,…，bKK1，…，bKKK]T.

定义函数向量X(s,t,u)=[x1(s,t,u),…,xN(s,t,u)]T，则所有N个样本数据的基函数展开可以表示为

X(s,t,u)=C(φ(s)⊗φ(t)⊗φ(u)),

(16)

(17)

则有

(18)

将式(17)和(18)代入到积分方程(11)中，得

(19)

进一步可得

(20)

上式也可以写成

(21)

其中u=S-1b.

解得u之后根据b=Su，可得一系列正交的特征向量bj.将bj代入式(14)即可得到一系列的正交特征函数

βj(s,t,u)=[φT(s)⊗φT(t)⊗φT(u)]bjj=1,2,…，J.

(22)

定义扩展的内积为：

(23)

每组体数据集xi(s,t,u)可由这些特征函数表示如下：

(24)

2 基于三元FPCA和k- means的MRI图像分类

为了探究三元FPCA对颅脑MRI图像病变的特征提取能力，本文对同样一组数据，首先进行图像预处理，得到脑脊液图像，然后分别用三元FPCA和传统PCA对脑脊液图像进行特征提取，再对提取的特征进行特征选择，最后利用k- means进行聚类，将分类结果与实际诊断信息对比，从而得到分类准确率，用于比较和分析实验结果，具体步骤如下.

2.1磁共振成像数据

本文使用的38例颅脑MRI由江苏省泰州市人民医院提供，其中正常13例，25例异常均为非急性期检查者.在25例异常中，脑卒中24例(缺血性脑梗死15例，腔隙性脑梗死9例)，脑肿瘤1例.脑卒中诊断参照2007年《中国脑血管病防治指南》，排除其他非血管性原因(如硬膜下血肿、癫痫后麻痹、脑外伤等)造成的脑功能障碍.正常个体指无脑血管病、自身免疫病、血液系统疾病、肿瘤、肝肾功能不全等疾病的个体.所有被试者体内均无金属植入物或没有MRI检查禁忌症，自愿参加本研究并签署知情同意书.本研究已获得复旦大学生命科学学院伦理委员会批准.

使用Siemens公司生产的Veior 3.0T磁共振成像系统对被试者脑部进行MRI扫描.常规三平面定位后，采用快速自旋回波(Turbo Spin Echo, TSE)脉冲序列，获得T2加权像(T2 Weighted Image, T2 WI)，其扫描参数如下： TR=3220ms，TE=93ms，FA=150°，层数： 56～60层，层厚=3mm，FOV=22cm×22cm，矩阵=256×256，Nex=1.

2.2MRI数据格式转换及预处理

采用MRIcroN软件(http:∥www.cabiatl.com/mricro/mricron/dcm2nii.html)将原始DICOM图像转换成NIfTI文件格式.转换格式后，使用SPM8[10](Statistical Parametric Mapping, Wellcome Trust Centre for Neuroimaging, London, UK)软件对T2 WI进行预处理，首先采用仿射变换(affine transformation)和非线性形变(non- linear deformation)进行空间标准化，将实验数据标准化到SPM8自带的T2模板上，体素大小重采样至1.5mm×1.5mm×1.5mm；然后对标准化后的脑结构图像进行解剖分割，得到灰质、白质和脑脊液3种图像；最后对分割后的图像进行高斯平滑(半高全宽参数=8mm)，以提高信噪比，使数据更接近高斯随机场模型.

2.3预处理后的脑脊液图像

本研究采用的25例患者均在非急性期接受MRI检查，除1例脑肿瘤外，其余皆为缺血性脑梗死(15例)和腔隙性脑梗死(9例).而缺血性脑梗死在MRI的T2 WI上的表现取决于梗死的时间，过了急性期后大的梗死灶最终形成囊性软化灶，表现近于脑脊液的信号[11]；腔隙性脑梗死在慢性期，MRI的T2 WI呈高信号，与脑脊液信号类似[12].此外，对样本图像进行预处理后得到的脑脊液图像，也显示其中保留了异常病灶信号(图1～图3(见第46页)所示).

因此，针对样本MRI图像的特点以及要筛查出的病变，我们考虑利用分割后的脑脊液图像做后续的特征提取.

图1 正常人原始图像及分割后的组织Fig.1 Original image of a normal person and its segmented tissues

图2 腔隙性脑梗死原始图像及分割后的组织Fig.2 Original image of a patient with lacunar infarction and its segmented tissues

图3 脑肿瘤原始图像及分割后的组织Fig.3 Original image of apatient with brain tumor and its segmented tissues

2.4MRI图像特征提取及聚类

对每例MRI体数据，分别用三元FPCA和传统PCA提取主成分，再通过确定性贪心选择算法(deterministic greedy selection algorithm for column subset selection)[13]对提取的主成分进行特征选择，从而得到最终的泛函主成分系数矩阵SFPCA和传统主成分系数矩阵SPCA.

采用k- means聚类算法[14]对主成分系数矩阵进行多次聚类(100次)，从中选择最优的结果，得到分类标签向量.对比分类得到的标签向量和实际诊断结果，计算准确率，用于比较和评价两种特征提取方法对分类准确率的影响.

3 结果与分析

实验结果(表1)表明，利用三元FPCA进行k- means分类的结果良好，并在准确率和特异性上优于PCA方法.实验最后使用了96个泛函主成分进行聚类分析，这些主成分的累积方差贡献率达到95%.图4比较了两种方法前1～5个主成分对应的累积方差贡献率.图4显示，用三元FPCA得到的第一个主成分方差贡献率就已经达到0.49，而PCA仅为0.09，远不及三元FPCA的特征提取能力.另外，PCA提取到第5个主成分的累积贡献率也才为0.29，而FPCA为0.61.

表1 三元FPCA和PCA的分类结果比较

以上结果表明，较之传统PCA，三元FPCA在对图像进行数据降维和特征提取方面更有优势.原因在于该方法将每例MRI体数据集看作一个连续的整体，考虑到了每例被试不同切层图像之间的相邻关系，对所有切层图像数据进行了整体的降维，从而使图像中差异最大的特征得到了更大程度的保留.而传统PCA对每例被试的每层图像数据进行“独立”地降维，并不能发现数据集中的非线性结构.此外，上述结果也验证了基于三元FPCA进行颅脑MRI图像分类识别的可行性和有效性.

图4 FPCA、PCA前1～5个主成分的累积方差贡献率比较Fig.4 Comparison between the cumulative contribution of variances of the first 5 principal components obtained by FPCA and that obtained by PCA

4 结 语

本文以颅脑MRI疾病筛查为研究目标，在图像特征提取方面将一元FPCA推广到了三元FPCA，再结合已有的MRI图像预处理软件，建立了一整套MRI病理/正常图像分析方法.该方法首先对样本图像进行图像配准和分割，然后运用三元FPCA对脑脊液进行特征提取，紧接着再对特征进行选择，随后用k- means对已选择的特征进行聚类，从而达到颅脑MRI图像病变筛查的目的.将该方法应用于38例颅脑MRI数据集上，结果证实三元FPCA较传统PCA有更好的数据降维和特征提取能力，有效地提高了MRI图像分类的准确率及效率.本文所提出的一整套MRI病理/正常图像聚类分析方法，可以实现对颅脑MRI从图像配准、图像分割、图像特征提取、分类识别的全过程控制，为未来大规模应用到MRI图像中提供了可行性.

[1] JOLLIFFE I T. Principal component analysis [M]. New York: Springer, 2002.

[2] HYVARINEN A, OJA E. Independent component analysis: Algorithms and applications[J].NeuralNetworks, 2000,13(4- 5): 411- 430.

[3] RAMSAY J O, DALZELL C J. Some tools for functional data analysis[J].JournaloftheRoyalStatisticalSociety, 1991,53(3): 539- 572.

[4] VIVIANI R, GRON G, SPITZER M. Functional principal component analysis of fMRI data [J].HumanBrainMapping, 2005,24(2): 109- 129.

[5] LONG C J, BROWN E N, TRIANTAFYLLOU C,etal. Nonstationary noise estimation in functional MRI [J].Neuroimage, 2005,28(4): 890- 903.

[6] RAMSAY J O, SILVERMAN B W. Functional data analysis [M]. New York: Springer, 1997.

[7] 靳刘蕊.函数性数据分析方法及应用研究[D].厦门: 厦门大学，2008.

[8] 李红广.基于B样条基粗糙惩罚的某些约束函数型数据光滑方法研究[D].上海: 华东师范大学，2008.

[9] SAGAN H. Introduction to the calculus of variations [M]. New York: Dover Publications, Inc, 1992.

[10] FRISTON K J, HOLMES A P, WORSLEY K J,etal. Statistical parametric maps in functional imaging: A general linear approach [J].HumanBrainMapping, 1994,2(4): 189- 210.

[11] 王维治.CT及MRI在脑血管疾病诊断中的应用及评价[J].中国实用内科杂志,2001,21(8): 496- 498.

[12] 金征宇.医学影像学[M].北京: 人民卫生出版社，2010: 69.

[13] CIVRIL A, MAGDON- ISMAIL M. Deterministic sparse column based matrix reconstruction via greedy approximation of SVD[C]∥Proc 19th International Symposium on Algorithms and Computation(ISAAC 2008). New York: Springer, 2008: 414- 423.

[14] MACQUEEN J B. Some methods for classification and analysis of multivariate observations[C]∥Proceedings of 5- th Berkeley Symposium on Mathematical Statistics and Probability[M]. Berkeley: University of California Press, 1967.

Abstract： One dimensional functional principal component analysis(FPCA) has been successfully applied in functional magnetic resonance imaging(fMRI) data. But there is scarce study focused on the application of high dimensional FPCA in magnetic resonance imaging(MRI) data. Three dimensional FPCA extended from one dimensional FPCA was provided and its application in MRI data to extract features was studied. A full set of techniques used to discriminate between pathological images and normal images was designed. It mainly consisted of the following three steps. Firstly, MRI data were preprocessed. The preprocessing steps included image registration and segmentation. The cerebrospinal fluid was obtained after image preprocessing. Secondly, functional PCA was employed to extract features from the segmented images. After that a feature selection method was performed for the extracted features. Thirdly, binary classification results were derived byk- means clustering method using features selected, which reached the goal to screen diseases. The application of this new method to brain MR turbo spin echo T2 weighted image data showed that FPCA outperformed multivariate PCA in feature extraction and classification accuracy.

Keywords： functional principal component analysis; principal component analysis; brain MRI; cluster analysis

ResearchonAlgorithmsBasedonThreeDimensionalFunctionalPrincipalComponentAnalysisandClusterAnalysisforMRIImageClassification

ZHANG Jiaming1, LIN Nan2, ZHANG Menghan1, ZHANG Liang3, LI Shuyuan1, XIONG Momiao2, WANG Xiaofeng1, JIN Li1

(1.SchoolofLifeSciences,FudanUniversity,Shanghai200438,China; 2.SchoolofPublicHealth,UniversityofTexas,HoustonTX77030,USA; 3.DepartmentofCardiology,ChanghaiHospital,SecondaryMilitaryMedicalUniversity,Shanghai200433,China)

R4；O1；TP3

A

0427- 7104(2017)01- 0040- 08

2016- 04- 05

科技部国际科技合作专项(2014DFA32830)；国家科技支撑项目(2011BAI09B02)

张嘉茗(1991—)，女，硕士研究生；金 力，男，教授，博士生导师，通信联系人，E- mail: lijin@fudan.edu.cn.