CN107844804A - 一种肌纹理磁共振图像的计算方法及基于其的数据建模方法 - Google Patents

Abstract

本发明是一种肌纹理磁共振图像的计算方法及基于其的数据建模方法，属于图像处理技术领域；通过计算机内的MATLAB工具进行小波变换，提取步骤1所获得磁共振图像的12个特征参数；采用支持向量机（Support Vector Machine,SVM）机器学习方法和人工神经网络（Artificial Neural Network,ANN）两种机器学习方法进行数据建模，对肌纹理磁共振图像进行分类预测。

Description

一种肌纹理磁共振图像的计算方法及基于其的数据建模方法

技术领域

本发明是一种肌纹理磁共振图像的计算方法及基于其的数据建模方法，属于图像处理技术领域。

技术背景

假性肥大型进行性肌营养不良症包括两种临床类型，即杜兴氏肌营养不良（Duchenne Muscular Dystrophy，DMD）和贝克氏肌营养不良（Becker MuscularDystrophy，BMD），二者是等位基因的异质性疾病。其中，DMD是一种严重的、X染色体基因遗传的、进行性骨骼肌营养不良的罕见疾病。DMD多发于新生或幼年男性，该疾病以骨骼肌进行性坏死为主要病理特征，是致死性骨骼肌遗传性疾病之一。DMD患者一般在3~5岁起病，13~16岁时丧失独立行走能力，平均死亡年龄为19岁。

目前DMD疾病的有效治疗方法尚未大规模普及，故高效准确地检测该疾病显得非常重要。DMD的传统检测流程一般为：临床疑似->血清肌酸激酶增高->肌电图检查呈肌源性损害->肌肉病理检查或基因检查证实。

现有的DMD临床检测方法包括：观测受试者的临床症状或体征、血清肌酶检测、肌电图、肌活检、基因检测、超声检查及磁共振成像检查等。

1. 血清肌酶检测（Creatine Kinase，CK）

血清肌酶检测法简便实用，对DMD的筛查和检测具有重要价值。但是，肌酶的变化会随患儿年龄的增长呈现出先升高后下降的趋势，故单纯采用血清肌酶检测法不足以判定受试者是否罹患DMD疾病，一般须结合肌电图及肌活检病理检查进行综合判断。

2. 肌电图（Electromyogram，EMG）

20世纪50年代开始，肌电图技术逐步应用于临床，这是检测肌肉肌源性损害和神经源性损害的重要辅助手段。该方法对进行性肌营养不良的检测具有重要意义，为进一步的病理或基因检测提供依据。

作为肌源性肌肉损害的代表性疾病之一，DMD患者的肌电图检查呈现出较为典型的肌源性损害。虽然其异常表现不具有病因学特异性，但可结合病史、体征和肌酸激酶检查，对DMD做出较为正确的判断。

然而，肌电图检查也有不足之处，主要表现在：

（1）肌电图为创伤性检查，会带给病人剧烈的疼痛感，难以被病人（尤其是儿童患者）接受，甚至无法完成检查，并因此影响到检查结果的可靠性。

（2）该检查的灵敏度不高。在进行肌电图检查时，针电极所能探测到的肌肉范围有较为限，故只有当病变达到一定程度范围时，肌电图才有所表现，因此该检查不易对DMD的早期症状做出精确的检测。

（3）该检查对操作者的要求较高。不论是进针的方式，还是对仪器参数的选择，或是对数据的采集及分析，都受到经验因素的制约，并带有一定主观性。因此，并非所有的肌电图检查都可得到正确结论。

3. 肌活检（Muscle biopsy）

一直以来，肌活检（肌肉病理检测）都被认为是许多肌病检测的金标准，免疫病理检查尤为重要，可明确多种肌病或神经源性疾病。但肌活检方法也有其不足之处，比如会受到取样误差的限制，并导致检测准确性下降。特别是当用穿刺术取代手术活检时，该误差对结果的影响会更明显。

4. 超声检查（Ultrasonography）

近年来，有国外学者尝试对肌肉疾病患者的发病部位进行超声成像，并借助计算机辅助手段，通过计算儿童神经肌肉疾病（Neuromuscular Disorders in childhood, NMD）超声图像的灰度特征值，来判断受试者是否患有NMD。结果显示，利用该方法得出的试验灵敏度和特异度分别是87%和67%。

国外学者曾对DMD患者及健康对照组的超声图像进行了研究。该项研究中，DMD患者组及健康对照组均为20人且年龄匹配。研究比较了两组超声组织特征（UltrasonicTissue Characterization，UTC）的心肌参数，发现两组之间存在着显著差异性，可据此作为受试者是否罹患DMD的辅助参考标准之一。

然而，尽管分析肌肉疾病患者的超声图像并得出预测结论是一种有效可行的方法，但总体来说，超声图像的质量令人不甚满意，图像对比度较差，信噪比不佳，且成像质量在一定程度上依赖于操作人员。

因此在DMD疾病的早期临床检测中，迫切需要一种创新的、无创的、灵敏的、有效的手段。如磁共振图像（MRI）检测方法，利用计算机对磁共振图像的自动化识别及数据建模作为辅助检测手段，其原理是通过DMD患者与正常人之间的磁共振图像比对，标注出其差异和该疾病的MRI纹理特征，方便医生进行疾病诊断。

然而，由于DMD疾病是一种罕见病，尚没有针对性的肌纹理磁共振图像计算方法，也缺少肌纹理磁共振图像的对比数据库。

发明内容

本发明的目的在于：提供一种针对DMD疾病等罕见肌肉疾病的肌纹理磁共振图像的计算方法。

本发明的目的通过下述方案实现：一种肌纹理磁共振图像的计算方法，包括以下步骤：

步骤1，通过磁共振成像设备获得磁共振图像输入计算机，或将现有磁共振图像输入计算机，所述磁共振图像为T1加权图像或T2加权图像；

步骤2，通过计算机内的MATLAB工具进行小波变换，提取步骤1所获得磁共振图像的12个特征参数，即水平均值、水平方差和水平能量，垂直均值、垂直方差和垂直能量，对角均值、对角方差和对角能量，全部均值、全部方差和全部能量，作为该磁共振图像的12个纹理特征参数。

在上述方案的基础上，优选的步骤1所述磁共振图像为T1加权图像。

在上述方案的基础上，步骤2包括以下步骤：

步骤2a，使用MATLAB读取磁共振图像，获得图像矩阵X；

步骤2b，进行小波分解，小波分解函数如下：

[c,s] = wavedec2(X,N,wavename)，

其中，c是分解系数矩阵，s是相应分解系数的长度矢量矩阵，X是图像矩阵，N是小波分解层数，wavename是小波基函数，小波分解层数N是2或3或4，小波基函数wavename分别是db4和sym4，即消失矩均取4；

步骤2c，提取第N层的水平小波系数H1D、垂直小波系数V1D、对角小波系数D1D和全部系数HVD，其中，

H1D=c(1+s(1,1)*s(1,2):s(1,1)*s(1,2)+s(2,1)*s(2,2))，

V1D=c(1+s(1,1)*s(1,2)+s(2,1)*s(2,2):s(1,1)*s(1,2)+2*s(2,1)*s(2,2))，

D1D=c(1+s(1,1)*s(1,2)+2*s(2,1)*s(2,2):s(1,1)*s(1,2)+3*s(2,1)*s(2,2))，

HVD=[H1D,V1D,D1D]；

步骤2d，分别对水平小波系数、垂直小波系数、对角小波系数和全部系数求取均值、方差和能量，获得水平均值、水平方差和水平能量，垂直均值、垂直方差和垂直能量，对角均值、对角方差和对角能量，全部均值、全部方差和全部能量，作为该磁共振图像的12个纹理特征参数。

在上述方案的基础上，优选的步骤2b所述小波分解层数N是2。

本发明还提供一种基于权利要求1至4任一所述肌纹理磁共振图像计算方法的数据建模方法，包括：采用支持向量机（Support Vector Machine, SVM）机器学习方法进行数据建模，使用留一法（leave-one-out cross validation）对单一对象或全部对象的全部肌纹理磁共振图像进行分类预测，分为DMD患者一类或健康者一类。

本发明还提供一种基于权利要求1至4任一所述肌纹理磁共振图像计算方法的数据建模方法，包括：采用支持向量机（Support Vector Machine, SVM）机器学习方法和人工神经网络（Artificial Neural Network, ANN）两种机器学习方法进行数据建模，使用十折交叉验证法（10-fold cross-validation）对单一对象或全部对象的全部肌纹理磁共振图像进行分类预测，分为DMD患者一类或健康者一类。

本发明的优点有：提供了针对DMD疾病等罕见肌肉疾病的肌纹理磁共振图像的计算方法；提供了可用于构建DMD疾病等罕见肌肉疾病对比数据库的数据建模方法。

附图说明

图1为本发明的图像计算及采用留一法的数据建模流程示意图；

图2为本发明的图像计算及采用十折交叉验证法的数据建模流程示意图。

具体实施方式

下面结合附图对本发明的实施例作详细说明。

如图1为本发明的图像计算及采用留一法的数据建模流程示意图，图2为本发明的图像计算及采用十折交叉验证法的数据建模流程示意图：

研究对象

本研究涉及的所有MRI均来自上海交通大学附属新华医院小儿神经科门诊，DMD患者组与健康儿童对照组的人数均为15例，均为男性，其年龄范围分别是8.4±1.7岁与8.8±1.9岁。

DMD患者组均为经肌肉活检后，病理结果确诊为DMD的患者；健康儿童对照组均为运动发育正常、无慢性疾病史、无遗传性疾病家族史的志愿者。本项研究已获得新华医院伦理委员会同意，患者家属均已知情同意。

图像计算

步骤1，根据DMD病变位置选定分析图像，获取有效的MRI共为485幅。其中T1图像244幅，包括127幅患者图像和117幅健康儿童图像；T2图像241幅，包括122幅患者图像和119幅健康儿童图像。截取对图像感兴趣的区域（ROI），大小约为180像素*180像素。

步骤2，通过计算机内的MATLAB工具进行小波变换，提取步骤1所获得磁共振图像的12个特征参数，即水平均值、水平方差和水平能量，垂直均值、垂直方差和垂直能量，对角均值、对角方差和对角能量，全部均值、全部方差和全部能量，作为该磁共振图像的12个纹理特征参数。

再分别利用支持向量机和人工神经网络，基于这些参数进行建模分析，以实现对图像的分类预测。

关于小波变换，所涉及的步骤如下：

1. 小波基选取

选择db和sym作为小波基函数，其消失矩均选定为4。即，将每一幅图像都按照db4和sym4函数进行分解，得到图像的水平、垂直和对角这3组原始小波系数。再对这些小波系数进行初步统计，共得到12个纹理特征统计量。最后根据结果进行分析比对，确定选用哪种小波基进行图像的小波分解更适合本课题的研究对象。

2. 分解层数的选择

小波分解的层数L也称尺度数（L即为小波分解层数N）。为进行对比，课题拟对每一幅图像进行L=2、L=3及L=4三层分解，比较在哪种小波分解层数下，结果更具说服力。

3. MRI T1加权图像和MRI T2加权图像

由于成像方法不同，故每一位受试者都有T1和T2两类加权图像。分别对两类图像进行处理，比较哪种加权图像更适合进行DMD疾病的计算机辅助检测。

4. 利用MATLAB工具提取纹理特征

提取对每一幅图像的基本小波系数，再对这些系数进行统计处理，共得到每一幅图像的12个纹理特征统计参数。对图像进行批处理，获取所有图像的全部纹理特征参数。

关键函数解释如下：

（1）二维图像分解

函数：[c,s] = wavedec2(X,level,wavename);

这是利用小波基函数wavename对灰度图像x进行层数为level的分解（level即为小波分解层数N或L）。其中，wavedec2( )是多层二维离散小波分解函数，其基本格式之一为[c,s] = wavedec2(X,N,’wname’)。该函数利用小波基函数’wname’对图像矩阵X在第N层上进行二维离散小波分解，其中N取值为正整数，返回结果为分解系数矩阵C和相应分解系数的长度矢量矩阵S。

分解系数矩阵C的结构是：

其中，A是行向量对应图像矩阵X的近似系数，H是细节系数矩阵的水平分量，V是细节系数矩阵的垂直分量，D是细节系数矩阵的对角分量。A(N)代表第N层的低频系数，H(N)|V(N)|D(N)代表第N层的高频系数，分别是水平、垂直及对角高频，依次类推到H(1)|V(1)|D(1)。

相应的分解系数长度矩阵S的结构是：

S(1,:)是A(N)的长度（即A(N)的原矩阵的行数和列数）；

S(2,:)是H(N)|V(N)|D(N)的长度；

S(3,:)是H(N-1)|V(N-1)|D(N-1)的长度；

……

S(N+1,:)是倒数第二行H(1)|V(1)|D(1)的长度；

S(N+2,:)是倒数第一行X的长度（大小）。

（2）提取小波系数

H1D=c(1+s(1,1)*s(1,2):s(1,1)*s(1,2)+s(2,1)*s(2,2));

V1D=c(1+s(1,1)*s(1,2)+s(2,1)*s(2,2):s(1,1)*s(1,2)+2*s(2,1)*s(2,2));

D1D=c(1+s(1,1)*s(1,2)+2*s(2,1)*s(2,2):s(1,1)*s(1,2)+3*s(2,1)*s(2,2));

HVD=[H1D,V1D,D1D];

提取第Level层的水平小波系数H1D、垂直小波系数V1D、对角小波系数D1D和全部系数HVD。

（3）对原始小波系数进行统计，以水平方向小波系数为例：

temparr(1)=mean(H1D);

temparr(2)=std(H1D);

HH=H1D.*H1D;

temparr(3)=mean(HH);

分别对所得的小波系数进行统计，mean为求取均值的函数，std为求取标准差的函数，对系数的平方求均值则可求取能量。

选择db和sym作为小波基函数，其消失矩阵的阶数选定4，分别对MRI T1和T2图像进行Level=2/3/4三种尺度的分解。由小波变换过程可知，对一幅图像每做一次（层）小波变换，可获得水平（H）、垂直（V）和对角（D）共3组小波系数，这是最基本的小波变换系数。

进一步地，实验者根据以往研究经验，分别对水平、垂直、对角以及全部小波系数求取均值、方差和能量，可获得每幅图像在第i层上的12个特征参数，即水平均值、水平方差和水平能量；垂直均值、垂直方差和垂直能量；对角均值、对角方差和对角能量；全部均值、全部方差和全部能量。这些数值将作为本实验暂定的纹理特征参数。

数据建模

为探索何种图像数据更适用于本研究的建模预测，分别对T1和T2两类图像进行三种尺度的小波分解，从每幅图像中均提取出12个纹理特征参数。之后依次以受试者个体和所有磁共振图像为研究对象，分别采用留一法和十折交叉验证法来组织样本数据，利用ANN和SVM两种机器学习算法对样本数据进行建模，实现对未知属性（患者或非患者）图像的分类预测。

对图1的展开阐述，以受试者个体为研究对象：

采用db4作为小波基函数，分解层数Level=2/3/4，图像包括T1和T2两种，对每位受试者所有图像（需要区分T1和T2）的12个参数求取均值，故每位受试者得到一条记录。例如，某受试者共有20张MRI，包括T1和T2图像各10张，这两种类型的图像需分开处理。以T1图像为例，对10张T1图像的12个纹理特征参数分别求取均值，最终得到一条记录。

利用SVM进行机器学习和分类预测时，采用RBF径向基核函数和多项式核函数，取（C,γ）的默认值，即C=1，γ=1/K=0.0833（其中K=12，即每条记录中的参数数目）。由于以受试者个体为研究对象时，样本数量较少，故适合采用“留一法”进行交叉验证。

1）样本类别号

本研究针对的是二值分类的问题，即受试者的每一幅图像或归属于DMD患者一类，或归属于健康者一类。故此，对训练样本集中的DMD患者赋类别号1，对健康者赋类别号-1。经SVM建模预测后，若测试样本集中的对象类别号为1，则判断该受试者为DMD患者；若为-1，则判断该受试者为健康者。

2）样本数据的格式

使用LIBSVM平台进行分类预测，必须将样本数据转化为相应的格式。实验中，每位受试者对象都有12个纹理特征参数（的均值），在每个特征参数前分别加上1~12，可得到合适的数据格式。即：

“1：参数1” “2：参数2”……“11：参数11” “12：参数12”

3）样本数据的归一化（标准化）

样本数据的归一化处理是数据挖掘的基础工作之一。由于不同的评价指标往往具有不同的量纲或数量级，从而会对数据分析的结果产生不良影响。为解决此类问题，需进行样本数据的归一化处理。原始数据经过归一化后，各指标将会处在同一数量级，更适合进行综合比对和评价。

在LIBSVM软件中，应用可执行文件SVM-scale即能完成数据归一化工作。SVM-scale的用法是：SVM-scale [-l lower] [-u upper] [-y y_lower y_upper] [-s save_ filename] [-r restore_filename] filename

几个关键参数为：

-l：数据的下限标记；lower：归一化后数据的下限。

-u：数据的上限标记；upper：归一化后数据的上限。

-s save_filename：将归一化的规则进行保存。

-r restore_filename：载入归一化规则文件并按此进行归一化。

filename：待归一化的数据文件。

本实验中，我们将所有样本数据都归一化在[-1,+1]范围内。值得注意的是，应该将训练样本集与测试样本集置于同一文件中并进行归一化，完成后再将归一化后的样本数据重新划分为训练样本集和测试样本集。

4）样本分组

我们以受试者个体为研究对象，采用“留一法”交叉验证来完成SVM的分类预测。本实验包括DMD患者正常人各15名，每次抽取1名受试者的记录作为预测样本，而将另外29名受试者的记录作为训练样本。依次重复30次，确保样本集中的所有样本均做过一次预测。

5）分析建模

在LIBSVM软件中，应用可执行文件SVM-train即能完成数据归一化工作。SVM-train的用法是：SVM-train [options] training_set_file [model_file]

几个关键参数（包含在options中）：

-t：设置核函数类型。0为线性核，1为多项式核，2为RBF径向基核。t值默认为2，即RBF多项式核。本研究试验了多项式核和RBF核，即对应的参数t分别取1和2。

-g：设置核函数中的γ。默认值为1/ K（K为每条记录中的数据项数，K=12。）

-c：设置惩罚系数C。默认值为1。

6）分类预测

在LIBSVM软件包中，应用可执行文件SVM-predict即能完成数据归一化工作。SVM- predict的用法：SVM-predict [options] test_file model_file output_file options

几个关键参数（包含在options中）：

-b probability_estimates：是否需要概率估计预测，可选0或1，默认值为0。

model_file 是SVM-train 产生的建模文件。

test_file是将要进行分类预测的数据文件。

output_file是由SVM-predict预测后的输出文件，内有预测结果。输出结果包括均方误差（Mean squared error）和相关系数（Squared correlation coefficient），如支持向量的个数等。

若需要调整核函数、惩罚参数C和核函数γ的值，可利用批处理过程对其进行选择，根据分类结果准确率筛选最优（C,γ）组合。本实验选用RBF和多项式作为SVM的核函数，同时设置多种(C,γ)组合。根据以往研究经验，分别选择C=1、10、100、1000；γ=0.083、5、80、1000。训练模型并实现对测试样本的分类预测，根据各项结果指标来评价16种不同(C,γ)参数组合下的分类情况，寻求推广能力最强的SVM模型。

对图2的展开阐述，以全体图像为研究对象：

采用db4作为小波基函数，分解层数Level=2/3/4，图像包括T1和T2两种，再以所有受试者的所有图像为研究对象（T1和T2需分开处理）。由于图片数量较多，故采用十折交叉验证法处理。将所有图像的数据记录随机分组，平均（或基本平均）分成10组，每次取其中一组作为测试样本集，而其余9组作为训练样本集。重复循环10次，直到所有样本集均做过一次测试样本集方才结束。其技术路线与留一法类似。

测试数据

表1 ANN模型的预测结果（个体为对象）

其中“T1_L2”表示对T1图像进行LEVEL=2的小波分解。其余以此类推。

若以全部图像为研究对象，结论有：

可以看出，“T1_L2”是最佳情况，且（1,0.083）是最佳参数组合。同样，当以“全部图像”为对象时，也有相同的结论。

在SVM方法下，以受试者个体为研究对象，结论有：

可以看出，“T1_L2”是最佳情况，且（1,0.083）是最佳参数组合。同样，当以“全部图像”为对象时，也有相同的结论。

综上所述，在对本研究项目中的特定DMD磁共振图像数据的SVM和ANN机器学习方法的分类方法进行了建模和计算，其结果表明，通过比较筛选，在这两类模型下，完全可以确定较好的参数组合，使得这些模型产生有效的分类识别效果，有把握让这些模型的分类的灵敏度、特异度及准确率均可达95%以上，从而获取较为满意的分类一致性和实用价值。

Claims (6)

1.一种肌纹理磁共振图像的计算方法，其特征在于，包括以下步骤：

步骤1，通过磁共振成像设备获得肌纹理磁共振图像输入计算机，或将现有肌纹理磁共振图像输入计算机，所述磁共振图像为T1加权图像或T2加权图像；

步骤2，通过计算机内的MATLAB工具进行小波变换，提取步骤1所获得磁共振图像的12个特征参数，即水平均值、水平方差和水平能量，垂直均值、垂直方差和垂直能量，对角均值、对角方差和对角能量，全部均值、全部方差和全部能量，作为该磁共振图像的12个纹理特征参数。

2.根据权利要求1所述的一种肌纹理磁共振图像的计算方法，其特征在于：步骤1所述磁共振图像为T1加权图像。

3.根据权利要求1所述的一种肌纹理磁共振图像的计算方法，其特征在于，步骤2包括以下步骤：

步骤2a，使用MATLAB读取磁共振图像，获得图像矩阵X；

步骤2b，进行小波分解，小波分解函数如下：

[c,s] = wavedec2(X,N,wavename)，

其中，c是分解系数矩阵，s是相应分解系数的长度矢量矩阵，X是图像矩阵，N是小波分解层数，wavename是小波基函数，小波分解层数N是2或3或4，小波基函数wavename分别是db4和sym4，即消失矩均取4；

步骤2c，提取第N层的水平小波系数H1D、垂直小波系数V1D、对角小波系数D1D和全部系数HVD，其中，

H1D=c(1+s(1,1)*s(1,2):s(1,1)*s(1,2)+s(2,1)*s(2,2))，

V1D=c(1+s(1,1)*s(1,2)+s(2,1)*s(2,2):s(1,1)*s(1,2)+2*s(2,1)*s(2,2))，

D1D=c(1+s(1,1)*s(1,2)+2*s(2,1)*s(2,2):s(1,1)*s(1,2)+3*s(2,1)*s(2,2))，

HVD=[H1D,V1D,D1D]；

步骤2d，分别对水平小波系数、垂直小波系数、对角小波系数和全部系数求取均值、方差和能量，获得水平均值、水平方差和水平能量，垂直均值、垂直方差和垂直能量，对角均值、对角方差和对角能量，全部均值、全部方差和全部能量，作为该磁共振图像的12个纹理特征参数。

4.根据权利要求3所述的一种肌纹理磁共振图像的计算方法，其特征在于：步骤2b所述小波分解层数N是2。

5.一种基于权利要求1至4任一所述肌纹理磁共振图像计算方法的数据建模方法，其特征在于：采用支持向量机（Support Vector Machine, SVM）机器学习方法进行数据建模，使用留一法（leave-one-out cross validation）对单一对象或全部对象的全部肌纹理磁共振图像进行分类预测，分为患者一类或健康者一类。

6.一种基于权利要求1至4任一所述肌纹理磁共振图像计算方法的数据建模方法，其特征在于：采用支持向量机（Support Vector Machine, SVM）机器学习方法和人工神经网络（Artificial Neural Network, ANN）两种机器学习方法进行数据建模，使用十折交叉验证法（10-fold cross-validation）对单一对象或全部对象的全部肌纹理磁共振图像进行分类预测，分为患者一类或健康者一类。