CN103745227A - 一种基于多维信息的肺结节良恶性鉴别方法 - Google Patents

Abstract

一种基于多维信息的肺结节良恶性鉴别方法，包括如下步骤：步骤1、三维结节的二维表示；步骤2、建立特征模型；步骤3、构建模糊分类器；步骤4、评价分类性能。本发明的有益效果：准确的特征建模对肺结节的良恶性鉴别至关重要。本发明采用影像学诊断特征、图像处理常用形状及纹理特征以及患者信息为肺结节良恶性鉴别提供更加客观的依据，特征的提取均在基于螺旋扫描技术生成的二维图像上进行并对特征建模采用了新的方法使得提取的特征更加准确。本发明使用模糊C均值聚类(FCM)算法对疑似结节进行良恶性鉴别，给出其为良性或恶性的概率，这样更符合医生思维模式。

Description

一种基于多维信息的肺结节良恶性鉴别方法

技术领域

本发明涉及一种基于多维信息的肺结节良恶性鉴别方法，尤其涉及特征建模以及基于模糊C均值聚类分类器构建。

技术背景

目前：人们在肺结节良恶性鉴别方面做了大量工作。在国外，Suzuki等人使用6个大规模训练的人工神经网络构成的CAD系统来对肺结节进行良恶性鉴别；Sumiaki等人利用线性判别函数对他们提出的量化收敛指数滤波器(QCI filter)的八个属性进行良恶性鉴别；Dolejei等人提出给漏诊的肺结节和错误分类的组织不同权重的一种自适应增强算法，降低了假阴率和假阳率；Antonelli等人模拟一个医师团队，建立一个多分类器系统，通过投票来完成肺结节的良恶性鉴别；Kumar等人提出采用模糊推理来对肺结节进行良恶性鉴别。在国内，魏颖等人采用基于加权的马氏距离来对肺结节进行良恶性鉴别；陈武凡等人利用支持向量机(SVM)来对肺结节进行良恶性鉴别；张丽秋等人利用反向传播的人工神经网络来进行良恶性鉴别。

影像学诊断特征、图像处理常用形状及纹理特征及患者信息为肺结节良恶性鉴别提供更加客观的依据。现有的肺结节良恶性计算机诊断具有唯一确定结果，即要么为恶性，要么为良性。但在临床诊断中，有时医生不能给出确切的诊断结果，而是给出肺结节为恶性或者良性的概率，这样的结果导致误诊的发生，延误病人的治疗。

发明内容

为了解决现有技术的不足，本发明提供一种基于多维信息的肺结节良恶性鉴别方法，通过对疑似结节进行良恶性鉴别，已达到给出其为良恶性的概率的目的。

本发明的方案是这样实现的：一种基于多维信息的肺结节良恶性鉴别方法，包括如下步骤：

步骤1、三维结节的二维表示；

采用螺旋扫描技术生成在候选结节表面上有序均匀分布的视点，在生成有序、均匀分布视点之后，以候选结节的中心为起始点向视点方向均匀的指定长度的射线，按顺时针方向将这些射线上的像素垂直地排列得到三维候选结节的二维图像。

步骤2、建立特征模型；

步骤2.2建立结节影像学特征模型；

(1)建立分叶特征模型

采用弦长与弦距比、特征点个数、最小凹角、可疑分叶点数、最大凹弧线比、分叶等级、最小凸角、最大凸弧线比来描述分叶征；

首先特征点检测，方法如下：a，根据二维图像分割边界，求各边界点与其相邻后边的边界点的斜率，获取各边界点斜率；b，判断特征点，具体为相邻两边界点斜率符号为正负则前点为凸；为负正则为凹；否则向下检测；将特征点依次放入集合A；c，优化特征点集合，检测相邻凹点之间的距离小于一定阈值且中间凸点也小于一定阈值，满足两条件则将三个特征点中的后两个从集合A中删除；

a表示A中的特征点，a_i表示A中的第i个特征点，特征点集A中三个相邻的特征点分别为a_i-1，a_i，a_i+1；a_i-1、a_i+1具有相同的凹凸性，与a_i相反；a_i-1与a_i+1之间的距离为弦长d_ci；a_i到直线a_i-1a_i+1的距离为弦距d_vi；其中，d表示距离，下角标c表示弦长，下角标v表示弦距，弦长与弦距比为

其中，R表示弦长与弦距比，R_i表示a_i对应的弦长与弦距比，以a_i为顶点所形成的夹角度数为α_i；

●获取平均凸弧线比

R_m=mean{R_i|a_i为凸点}   (1)

其中，R为弦长与弦距比，下角标m表示平均，R_m表示平均凸弧线比，a表示A中的特征点，a_i为特征点集A中第i个特征点；

●统计特征点个数

N_f=card A

其中，A为特征点集合，N表示个数，下角标f标识为特征点，N_f表示特征点个数；

●获取最小凹角与最小凸角

相邻三个特征点a_i-1，a_i，a_i+1对应的坐标分别为(x_i-1,y_i-1)，(x_i,y_i)，(x_i+1,y_i+1)；通过余弦定理可求得a_i为顶点所形成的夹角α_i；

$\cos {\mathrm{\α}}_i=\frac{{(x_{i+1}-x_i)}^2+{(y_{i+1}-y_i)}^2+{(x_{i-1}-x_i)}^2+{(y_{i-1}-y_i)}^2-({(x_{i+1}-x_{i-1})}^2+{(y_{+1}-y_{i-1})}^2)}{2\sqrt{{(x_{i+1}-x_i)}^2+{(y_{i+1}-y_i)}^2}\sqrt{{(x_{i-1}-x_i)}^2+{(y_{i-1}-y_i)}^2}}$

最小凹角大小：α_c=min{α_i|a_i为凹点}

最小凸角大小：α_t=min{α_i|a_i为凸点}

其中，a_i为特征点集A中第i个特征点，x_i，y_i分别表示特征点a_i在二维图像上的横坐标和纵坐标，α表示角度，α_i表示坐a_i为顶点所形成的夹角度数，下角标c表示凹角，t表示凸角，α_c表示最小凹角，α_t表示最小凸角；

●统计可疑分叶点数

N_s=card{R_i>0.5}

其中，R_i表示a_i对应的弦长应弦距比，N表示个数，下角标s标识为可疑的分叶，N_s表示特征点个数；

●获取最大凹弧线比应最大凸弧线比

最大凹弧线比：R_c=max{R_i|a_i为凹点}

最大凸弧线比：R_t=max{R_i|a_i为凸点}

其中，a_i为特征点集A中第i个特征点，R_i表示a_i对应的弦长与弦距比，R为弦长与弦距比，下角标c表示凹角，t表示凸角，R_c表示最大凹弧线比，R_t表示最大凸弧线比；

(2)建立毛刺特征模型

具体步骤如下：

①对二维图像进行伽马(gamma)校正；

②利用卷积函数imfilter求图像竖直边缘Ix，水平边缘Iy，然后根据竖直边缘Ix和水平边缘Iy来求得边缘强度sqrt(Ix.^2＋Iy.^2)和边缘斜率Iy./Ix；

③将图像每m'×n'个像素分到一个单元(cell)中，对于大小为M'×N'的二维灰度图像，可以分成至少两个单元，这里每个单元大小为16*16，图像大小为32*412；

④对于每个单元求其梯度方向直方图，通常取Q个方向，也就是每360/Q分到一个方向，方向大小按像素边缘强度加权，最后归一化直方图；

⑤每2*2个单元合成一个块(block)，所以这里就有(2-1)*(25-1)=24个块；

⑥每个块中都有2*2*9个特征，这里一共有24个块，所以24个36维的特征构成总特征；

(3)建立胸膜凹陷特征模型

胸膜凹陷征表现为肺内结节与邻近胸膜之间的一条或几条线状影，在结节的二维图像中，如果结节存在胸膜凹陷，则相应列上会出现上下两组边界，当出现胸膜凹陷征时，拟使用动态规划算法搜索2D图像的上边界，利用上边界所含区域的大小PC来衡量胸膜凹陷的程度；

其中，d为结节中心到肺边界的距离，d_max为结节的最大直径，如果

则该结节可能为与胸膜粘连的结节；如果d在[d_max/2,2d_max]范围内，则该结节可能存在胸膜凹陷征，用搜索边界上部区域面积PC来度量胸膜凹陷程度；如果d>2d_max，则该结节远离肺壁，不应胸膜粘连；

(4)建立空泡特征模型

空泡征表现为肺结节内的类圆形含气腔，在肺窗上呈低密度影，使用变量空泡体积百分比PE_va来表示

${\mathrm{PE}}_{\mathrm{va}}=\frac{V_{\mathrm{va}}}{V}---\left(3\right)$

其中，V为肺结节的体积，V_va为肺结节内空泡体积，由公式(8)计算得到；

在二维图像中，取结节边界以内的部分为感兴趣区域，对感兴趣区域进行阈值分割得到空泡边界，然后根据空泡边界求得空泡体积；

(5)建立边界清晰度特征模型

利用肺结节边界点的径向梯度和GR来表示结节边界清晰程度，如果径向梯度和越大，表明结节的边界越清晰，

$\mathrm{GR}=\frac{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}(f(j,\mathrm{boundary}\left(j\right)+1)-f(j,\mathrm{boundary}(j))}{N}---\left(4\right)$

其中，f为结节的二维灰度图像，j为边界中第几个边界点，boundary(j)为结节边界第j个边界点对应的纵坐标，N为边界点的数量；

(6)建立磨玻璃样特征模型

根据边界清晰程度与平均密度来度量结节的磨玻璃样，

$G_{\mathrm{go}}={\mathrm{\ω}}_1\frac{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}(f(j,\mathrm{boundary}\left(j\right)+1)-f(j,\mathrm{boundary}\left(j\right))}{N}+{\mathrm{\ω}}_2\frac{\underset{(x,y)\∈\mathrm{RE}}{\mathrm{\Σ}}f(x,y)}{N_{\mathrm{RE}}}---\left(5\right)$

其中，f为结节的二维灰度图像，boundary(j)为结节边界第j个边界点对应的纵坐标，N为边界点的数量，RE为二维图像中结节分割边界以内的区域，N_RE为RE区域像素点个数，ω₁为边界清晰程度调节系数，ω₂为平均密度的调节系数，ω₁,ω₂使二者处于相近的数量级；

(7)建立钙化含量征特征模型

采用钙化含量百分比PR_ca来量化结节的钙化程度，

${\mathrm{PR}}_{\mathrm{ca}}=\frac{V_{\mathrm{ca}}}{V}\×100\%---\left(6\right)$

其中，V为肺结节的体积，V_ca为肺结节内钙化部分体积，由公式(8)计算得到；

在二维图像中，取结节边界以内的部分为感兴趣区域，对感兴趣区域进行动态规划分割得到钙化边界，然后根据钙化边界求得钙化体积。

步骤2.3建立图像形状及纹理的特征模型；包括如下步骤：

(1)获取结节表面积

由于螺旋扫描技术产生的视点近似均匀的分布在结节表面，在结节分割后的图像中，距结节中心越远表面积越大，距离结节中心为r个像素的视点所在区域表面积是距结节中心1个像素视点所在区域表面积的r²倍，

$S=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}s\×{\left(\mathrm{boundary}\left(j\right)\right)}^2---\left(7\right)$

其中，N为结节边界点的数量，s为距结节中心1个像素视点所在表面积

boundary(j)为结节边界第j个边界点对应的纵坐标，由二维图像的生成原理可知，二维图像中结节边界点的纵坐标也表示边界点距离结节中心的距离；

(2)获取结节体积

肺结节的体积常用于估计肺结节的倍增时间，它是医生判断肺结节良、恶性的重要指标，在本专利中，由于螺旋扫描技术产生的视点近似均匀的分布在结节表面，因此，在结节分割后的图像中，距结节中心越远体积越大，距离结节中心为r个像素的视点所在体积是距结节中心1个像素视点所在体积的r³倍，

$V=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}v\×{\left(\mathrm{boundary}\left(j\right)\right)}^3---\left(8\right)$

其中，N为结节边界点的数量，v为距结节中心1个像素视点所在体积

boundary(j)为二维结节边界第j个边界点对应的纵坐标，由二维图像的生成原理可知，二维图像中结节边界点的纵坐标也表示边界点距离结节中心的距离；

(3)获取结节球形度

球形度是描述结节近似球形的程度，在[0，1]范围内，当球形度越接近1，则说明该结节越近似于球形；

$D=\frac{\sqrt[6]{36\mathrm{\π}}\sqrt[3]{V}}{\sqrt{S}}---\left(9\right)$

其中，V为肺结节的体积，由公式(8)计算可得，S为肺结节的表面积，由公式(7)计算可得；

(4)获取结节最大直径

最大直径d_max是描述结节大小的一个重要特征，由于螺旋扫描技术产生的视点近视对称的分布在结节表面，在结节分割后的图像中，过结节中心，且过表面上下、左右对称视点的射线长度为最大直径d_max；

(5)获取结节平均密度

$M_n=\frac{\underset{(x,y)\∈\mathrm{RE}}{\mathrm{\Σ}}f(x,y)}{N_{\mathrm{RE}}}---\left(10\right)$

其中，f(x,y)为二维图像位于(x,y)点的灰度值，RE为二维图像中结节分割边界以内的区域，N_RE为RE区域像素点个数，M_n表示结节平均密度；

(6)获取结节密度方差

$V_n=\sqrt{\frac{1}{N_{\mathrm{RE}}}\underset{(x,y)\∈\mathrm{RE}}{\mathrm{\Σ}}{(f(x,y)-M_n)}^2}---\left(11\right)$

其中，f(x,y)为二维图像位于(x,y)点的灰度值，RE为二维图像中结节分割边界以内的区域，N_RE为RE区域像素点个数，M_n表示结节平均密度，V_n表示结节密度方差；

(7)获取结节傅里叶描述子

结节边界是在二维图像上通过动态规划获得；

$\mathrm{FD}\left(u\right)=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\mathrm{boundary}\left(j\right)e^{-j^{\′}2\mathrm{\πu}(j-1)/N}---\left(12\right)$

其中，FD(u)为结节傅里叶描述子，u=0,1,2,…,N-1，N为结节边界点的数量，boundary(j)为结节边界第j个边界点对应的纵坐标，j'表示复系数，由二维图像的生成原理可知，二维图像中结节边界点的纵坐标也表示边界点距离结节中心的距离；

(8)获取平均标准化半径

提取平均标准化半径之前需要对半径进行归一化，

$d_{\mathrm{mr}}=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{d}_{\mathrm{no}}\left(j\right)---\left(13\right)$

其中，d_mr为平均标准化半径，d_no(j)为结节边界点距离结节中心的归一化距离，N为结节边界点的数量；

(9)获取标准化半径方差

$\mathrm{\σ}=\sqrt{\frac{1}{N-1}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(d_{\mathrm{no}}\left(j\right)-d_{\mathrm{mr}})}^2}---\left(14\right)$

其中，d_mr为平均标准化半径，d_no(j)为结节边界点距离结节中心的归一化距离，N为结节边界点的数量；

(10)获取标准化半径的熵E(entropy of the normalized radial)

$E=-\underset{k=1}{\overset{100}{\mathrm{\Σ}}}{p}_k\left(\log \left(p_k\right)\right)---\left(15\right)$

其中，式中的p_k表示概率值，这个概率值是通过标准直方图来得出的，标准直方图的获取方法是将最大半径与最小半径之间的差值平均划分为100等分，然后将所有的半径分别装入自己所属的区段，最后得到每个区段所占的比例，即p_k，它表示在第k个区段内的半径数所占的比例；

(11)获取面积比率A_r

$A_r=\frac{1}{d_{\mathrm{mr}}N}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}(d_{\mathrm{no}}\left(j\right)-d_{\mathrm{mr}})---\left(16\right)$

其中，d_no(j)-d_mr=0，

d_mr为平均标准化半径，d_no(j)为结节边界点距离结节中心的归一化距离，N为结节边界点的数量；

(12)获取粗糙度R(boundary roughness)

$\mathrm{Rough}=\frac{1}{N}\underset{j=1}{\overset{N-1}{\mathrm{\Σ}}}|d_{\mathrm{no}}\left(j\right)-d_{\mathrm{no}}(j+1)|---\left(17\right)$

其中，d_no(j)为结节边界点距离结节中心的归一化距离，N为结节边界点的数量；

(13)获取灰度共生矩阵统计量

采用灰度共生矩阵统计量对结节的二维图像进行分析，二维图像采用螺旋扫描技术生成，在二维图像中，水平相邻的两像素点距离结节中心越近，在三维结节中两像素点间的实际距离越近，

设f(x,y)为一幅二维数字图象，其大小为M×N，灰度级别为Ng，则满足一定空间关系的灰度共生矩阵为：

P(m,n)=#{(x1,y1),(x2,y2)∈M×N|f(x1,y1)=m,f(x2,y2)=n}    （18）

其中#(x)表示集合X中的元素个数，显然P为Ng×Ng的矩阵，(x1,y1)、(x2,y2)为图像中的两点，对应的灰度级分别为m，n，若(x1,y1)与(x2,y2)间距离为dist，两者与坐标横轴的夹角为θ，则可以得到各种间距及角度的灰度共生矩阵P(m,n,dist,θ)；

在提取灰度共生矩阵的特征参数之前，要作正规化处理，令p(m,n,dist,θ)=P(m,n,dist,θ)/H，其中H是正规化常数，是灰度共生矩阵中全部元素之和；

定义用于纹理分析的灰度共生矩阵特征参数如下：

●获取角二阶矩(Angular second moment)

$\mathrm{ASM}=\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}{p}^2(m,n,\mathrm{dist},\mathrm{\θ})---\left(19\right)$

其中，式中q为选取的灰度级Ng的大小；

角二阶矩也称能量，是灰度共生矩阵元素值的平方和，反映了图像灰度分布均匀程度和纹理粗细度，如果共生矩阵的所有值均相等，则ASM值小；相反，如果其中一些值大而其余值小，则ASM值大，当共生矩阵中元素集中分布时，此时ASM值大，ASM值大表明一种较均一和规则变化的纹理模式；

●获取对比度(Contrast)

$\mathrm{Constrast}=\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}\left[{(m-n)}^2{p}^2(m,n,\mathrm{dist},\mathrm{\θ})\right]---\left(20\right)$

其中，式中q为选取的灰度级Ng的大小，

对比度反映了图像的清晰度和纹理沟纹深浅的程度，纹理沟纹越深，其对比度越大，视觉效果越清晰；反之，对比度小，则沟纹浅，效果模糊，灰度差即对比度大的象素对越多，这个值越大，灰度共生矩阵中远离对角线的元素值越大，Constrast越大；

●获取相关性(Correlation)

$\mathrm{Correlation}=\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}(m\×n\×p(m,n,\mathrm{dist},\mathrm{\θ})-u_1\×u_2)/d_1{\×d}_2---\left(21\right)$

其中，式中q为选取的灰度级Ng的大小；

$\begin{array}{cc}{u}_1=\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}m\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}p(m,n,\mathrm{dist},\mathrm{\θ})& u_2=\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}j\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}p(m,n,\mathrm{dist},\mathrm{\θ})\end{array}$

$\begin{array}{cc}{d_1}^2=\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}{(m-u_1)}^2\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}p(m,n,\mathrm{dist},\mathrm{\θ})& {d_2}^2=\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}{(n-u_2)}^2\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}p(m,n,\mathrm{dist},\mathrm{\θ})\end{array}$

用于度量空间灰度共生矩阵元素在行或列方向上的相似程度，相关值大小反映了图像中局部灰度相关性，当矩阵元素值均匀相等时，相关值就大;相反，如果矩阵像元值相差很大则相关值小，如果图像中有水平方向纹理，则水平方向矩阵的Correlation大于其余矩阵的Correlation值；

●获取熵(Entropy)

$\mathrm{Entropy}=-\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}p(m,n,\mathrm{dist},\mathrm{\θ})\×\mathrm{lgp}(m,n,\mathrm{dist},\mathrm{\θ})---\left(22\right)$

其中，式中q为选取的灰度级Ng的大小；

熵是图像具有的信息量的度量，纹理信息也属于图像的信息，是一个随机性的度量，当共生矩阵中所有元素有最大的随机性、空间共生矩阵中所有值几乎相等时，共生矩阵中元素分散分布时，熵较大，它表示了图像中纹理的非均匀程度或复杂程度；

●获取方差(Variance)

$\mathrm{Variance}=\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}{(m-\mathrm{mean})}^{2}p(m,n,\mathrm{dist},\mathrm{\θ})---\left(23\right)$

其中，式中q为选取的灰度级Ng的大小，mean为p(m,n,dist,θ)的均值；

●获取方差和(Sum of variance)

$\mathrm{SOV}=\underset{l=2}{\overset{2q}{\mathrm{\Σ}}}{(l-\mathrm{SOA})}^2{p}_x\left(l\right)---\left(24\right)$

其中，式中q为选取的灰度级Ng的大小，SOA由公式(25)计算得到，

$p_x\left(l\right)=\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}p{(m,n,\mathrm{dist},\mathrm{\θ})}_{|m+n|=l}l=\mathrm{2,3},\·\·\·,2q;$

●获取均值和(Sum of average)

$\mathrm{SOA}=\underset{l=2}{\overset{2q}{\mathrm{\Σ}}}{l\×p}_x\left(l\right)---\left(25\right)$

其中，式中q为选取的灰度级Ng的大小，

$p_x\left(l\right)=\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}p{(m,n,\mathrm{dist},\mathrm{\θ})}_{|m+n|=l}l=\mathrm{2,3},\·\·\·,2q;$

均值和是图像区域内像素点平均灰度值的度量，反映了图像的明暗深浅，适用于灰度图像；

●获取逆差矩(Inverse difference moment)

$\mathrm{IDM}=\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}p(m,n,\mathrm{dist},\mathrm{\θ})/(1+{(m-2)}^2)---\left(26\right)$

其中，式中q为选取的灰度级Ng的大小；

逆差矩又称为局部平稳，它是图像纹理局部变化的度量，反映了纹理的规则程度，纹理越规则，IDM就越大；反之亦然；

●获取差的方差(Variance of difference)

$\mathrm{VOD}=\underset{l=0}{\overset{q-1}{\mathrm{\Σ}}}{[l-\underset{l=0}{\overset{q-1}{\mathrm{\Σ}}}l\×p_y\left(l\right)]}^2\×p_y\left(l\right)---\left(27\right)$

其中，式中q为选取的灰度级Ng的大小，

$p_y\left(l\right)=\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}p{(m,n,\mathrm{dist},\mathrm{\θ})}_{|m+n|=l}l=\mathrm{0,1},\·\·\·,q-1;$

差的方差是邻近像素对灰度值差异的方差，对比越强烈，VOD越大；反之亦然；

●获取和熵(Sum of entropy)

$\mathrm{SOE}=-\underset{l=2}{\overset{2q}{\mathrm{\Σ}}}{p}_x\left(l\right)\×\log \left[p_x\left(l\right)\right]---\left(28\right)$

其中，式中q为选取的灰度级Ng的大小，

$p_x\left(l\right)=\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}p{(m,n,\mathrm{dist},\mathrm{\θ})}_{|m+n|=l}l=\mathrm{2,3},\·\·\·,2q;$

●获取差熵(Difference of entropy)

$\mathrm{DOE}=-\underset{l=0}{\overset{q-1}{\mathrm{\Σ}}}{p}_y\left(l\right)\×\log \left[p_y\left(l\right)\right]---\left(29\right)$

其中，式中q为选取的灰度级Ng的大小，

$p_y\left(l\right)=\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}p{(m,n,\mathrm{dist},\mathrm{\θ})}_{|m+n|=l}l=\mathrm{0,1},\·\·\·,q-1;$

●获取聚类阴影(Clustering shadow)

$\mathrm{CS}=-\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}{[(m-u_1)+(n-u_2)]}^{3}p(m,n,\mathrm{dist},\mathrm{\θ})---\left(30\right)$

其中，式中q为选取的灰度级Ng的大小，

$\begin{array}{cc}{u}_1=\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}m\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}p(m,n,\mathrm{dist},\mathrm{\θ})& u_2=\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}n\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}p(m,n,\mathrm{dist},\mathrm{\θ})\end{array};$

●获取显著聚类(Significant clustering)

$\mathrm{SC}=-\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}{[(m-u_1)+(n-u_2)]}^{4}p(m,n,\mathrm{dist},\mathrm{\θ})---\left(31\right)$

其中，式中q为选取的灰度级Ng的大小，

$\begin{array}{cc}{u}_1=\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}m\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}p(m,n,\mathrm{dist},\mathrm{\θ})& u_2=\underset{m=0}{\overset{q-1}{\mathrm{\Σ}}}n\underset{n=0}{\overset{q-1}{\mathrm{\Σ}}}p(m,n,\mathrm{dist},\mathrm{\θ})\end{array};$

●获取最大概率(Maximum probability)

$\mathrm{MP}=\underset{m,n}{\mathrm{MAX}}\left(p(m,n,\mathrm{dist},\mathrm{\θ})\right)---\left(32\right).$

步骤2.4结合患者的相关信息。

包括如下信息：

(1)患者年龄；

(2)患者性别；

(3)患者是否患过其他癌症；

患者是否痰中咯血。

步骤3、构建模糊分类器；

具体如下：

通过模糊C-均值聚类分类器，并用马氏距离代替欧氏距离构建；

求下式中目标函数J的最小值：

$J(U,Z)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\underset{i=1}{\overset{C}{\mathrm{\Σ}}}{\left({\mathrm{\μ}}_{\mathrm{ij}}\right)}^m{d}^2(x_j,z_i)---\left(33\right)$

式中，C为类别数目；U={μ_ij}是C×N维模糊分类矩阵，m是一个常数，用于控制分类的隶属度，d²(x_j,v_i)表示样本x_j到第i类中心v_i的欧式距离，d²(x_j,z_i)表示样本x_j到第i类中心z_i的欧式距离，表示为d²(x_j,z_i)=||x_j-z_i||²。Z={z₁,z₂,...,z_c}是C维聚类中心矩阵，表示为 $z_i=\frac{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{\left({\mathrm{\μ}}_{\mathrm{ik}}\right)}^m{x}_k}{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{\left({\mathrm{\μ}}_{\mathrm{ik}}\right)}^m}1\≤i\≤C.$ μ_ij表示样本x_j属于聚类中心z_i的隶属度，表示为 ${\mathrm{\μ}}_{\mathrm{ij}}=\frac{1}{\underset{k=1}{\overset{C}{\mathrm{\Σ}}}{\left[\frac{d_{\mathrm{ij}}}{d_{\mathrm{kj}}}\right]}^{2/(m-1)}};$

在分类过程中，不断迭代来更新隶属度函数和聚类中心。当满足迭代终止条件(1)迭代次数小于100，或(2)目标代价函数值J小于1.2×10⁵，或(3)连续两次迭代产生的聚类中心的差值小于0.1时，停止迭代取得最终的聚类中心z_i；

设s维空间的样本集X={x₁,x₂,…,x_n}，样本x到样本集X的马氏距离如公式(34)所示；

$D_M\left(x\right)=\sqrt{{(x-\mathrm{\μ})}^T{\mathrm{\Σ}}^{-1}(x-\mathrm{\μ})}---\left(34\right)$

式中μ--样本集X的均值向量;

∑--为协方差矩阵， $\mathrm{\Σ}=\frac{1}{l}\underset{j=1}{\overset{l}{\mathrm{\Σ}}}(x_j-p){(x_j-p)}^T;$

在计算马氏距离的过程中，若协方差矩阵为奇异矩阵，则协方差矩阵的逆阵不存在，将马氏距离无法直接求得，根据矩阵理论，任一秩为r的实对称半正定矩阵∑都可按公式(33)分解；

∑=A^TG^-1A    (35)

式中G--为r×r的对角阵，它由∑的非0特征值构成，G是非奇异的；

A--为r×m矩阵，由G中特征值所对应的特征向量构成，它张成了∑的非退化子空间，且AA^T为r×r的单位阵；

通过对公式(35)的分解，就可根据G的逆来求∑的伪逆矩阵，如公式(36)所示；

∑^-1=A^TG^-1A    (36)

用马氏距离来代替欧氏距离应用于模糊C均值聚类算法中；

将公式(34)代入到传统模糊C均值聚类算法的全局代价函数(33)中，得到基于马氏距离的模糊C均值聚类算法的全局代价函数，如公式(37)所示，约束条件如公式(38)、(39)所示：

$J_{\mathrm{FCM}}(U,Z)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\underset{i=1}{\overset{C}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ij}}^m{d}^2(x_j-z_i)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\underset{i=1}{\overset{C}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ij}}^m{(x_j-z_i)}^T{\mathrm{\Σ}}^{-1}(x_j-z_i)---\left(37\right)$

$\begin{array}{cc}0\≤{\mathrm{\μ}}_{\mathrm{ij}}\≤1& \∀i\∈\end{array}\{\mathrm{1,2},\·\·\·,c\},j\∈\{\mathrm{1,2},\·\·\·,n\}---\left(38\right)$

$\begin{array}{cc}\underset{i=1}{\overset{C}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ij}}=1& \∀j\∈\{\mathrm{1,2},\·\·\·,n\}\end{array}---\left(39\right)$

依据GK算法，聚类中心、协方差矩阵和隶属度值的计算公式分别如公式(40)、(41)、(42)所示；

$z_i=\frac{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ij}}^m{x}_j}{\underset{i=1}{\overset{C}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ij}}^m},1\≤i\≤C---\left(40\right)$

${\mathrm{\Σ}}_i=\frac{\underset{j=1}{\overset{l}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ij}}^m(x_j-p_i){(x_j-p_i)}^T}{\underset{j=1}{\overset{l}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ij}}^m},1\≤i\≤C---\left(41\right)$

(42)

式中D_ij--为马氏距离 $D_{\mathrm{ij}}=\sqrt{{(x_j-z_i)}^T{\mathrm{\Σ}}_I^{-1}(x_j-z_i)}.$

所述的模糊C均值聚类算法的算法如下：

(1)设置目标函数精度ε，模糊指数m，最大迭代次数T_m；

(2)初始化模糊聚类中心；

(3)由步骤(3)-(5)更新模糊划分矩阵U={μ_ij}和聚类中心Z={z₁,z₂,...,z_c}；

(4)若|J(t)-J(t-1)|<ε或t>T_m，则结束聚类；t←t+1并转到(3)步；

(5)由U={μ_ij}得到各像素点分类结果。

步骤4、评价分类性能，具体如下：

采用的样本数据来自肺部图像数据库联盟(LIDC)，同时采用敏感性、准确性、特异性和假阳性四个指标来分别评价分类性能，具体为：

(1)根据敏感性(Sensitivity)(又称为真阳性率)，判断对结节样本恶性的检测性能；

$\mathrm{Sensitivity}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}---\left(43\right)$

(2)根据准确性(Accuracy)，判断所有的检测结果中，检测恶性的样本占全部样本的比率；

$\mathrm{Accuracy}=\frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{FP}+\mathrm{TN}+\mathrm{FN}}---\left(44\right)$

(3)根据特异性(Specificity)，判断对结节样本良性的检测性能；

$\mathrm{Specificity}=\frac{\mathrm{TN}}{\mathrm{TN}+\mathrm{FP}}---\left(45\right)$

(4)根据假阳性率，判断检测结果中被错误检测为良性结节样本所占的比率；

$\mathrm{FPF}=\frac{\mathrm{FP}}{\mathrm{TN}+\mathrm{FP}}---\left(46\right)$

上式中：TP为真阳性结节的数量，FP为假阳性结节的数量，TN为假阴性结节的数量，FN为真阴性结节的数量。

本发明的有益效果：准确的特征建模对肺结节的良恶性鉴别至关重要。本发明采用影像学诊断特征、图像处理常用形状及纹理特征以及患者信息为肺结节良恶性鉴别提供更加客观的依据，特征的提取均在基于螺旋扫描技术生成的二维图像上进行并对特征建模采用了新的方法使得提取的特征更加准确。现有的肺结节良恶性计算机诊断具有唯一确定结果。但在临床诊断中，有时医生不能给出确切的诊断结果，而是给出肺结节为良性或恶性的概率。本发明使用模糊C均值聚类(FCM)算法对疑似结节进行良恶性鉴别，给出其为良性或恶性的概率，这样更符合医生思维模式。

附图说明

图1为本发明的总流程图；

图2为有序均匀分布的视点示意图；

图3为三维候选结节的二维图像；

图4为特征点检测流程图；

图5为含胸膜凹陷征结节的2D整体图像；

图6为本发明的特征提取方案的结构图；

图7为本发明的毛刺特征提取的流程图；

图8为本发明的模糊C-均值聚类训练的流程图。

具体实施方式

本发明结合具体实施方式进行说明：

一种基于多维信息的肺结节良恶性鉴别方法，包括如下步骤：如图1所示，

步骤1、三维结节的二维表示；

步骤2、建立特征模型；

步骤3、构建模糊分类器；

步骤4、评价分类性能。

步骤1、三维结节的二维表示

本发明采用螺旋扫描技术生成在候选结节表面上有序均匀分布的视点。视点的产生是由方位角、仰角来确定。方位角的范围为0到2π，仰角的范围为0到π。由于方位角范围是仰角范围为的二倍，初始的方位角及仰角分别被分成2N份和N份，这样就产生2N²个分布在球面的视点。螺旋化就是让每个试点都处在不同的纬度上，而不再是2N个视点集中在一个平面上。在第一阶段生成了2N²个视点，因此2N²个视点应该处在不同的纬度上。对螺旋线上的视点进行重采样，采样间隔为4N²/π，可以获得有序均匀分布的视点，如图2所示，其中(a)具有均分方位角和仰角的视点，(b)使用螺旋扫描技术使的所有视点按照仰角均匀分布，(c)对于螺旋线上的点重新采样，使得视点均匀分布。

在生成有序、均匀分布视点之后，以候选结节的中心为起始点向视点方向均匀的指定长度的射线。按顺时针方向将这些射线上的像素垂直地排列得到三维候选结节的二维图像如图3所示，其中(a)结节的CT断层图像，(b)结节的2D整体图像。在生成二维图像时，采用双线性插值算法来获得图像中像素点的灰度值。

本发明在特征建模所用到的二维图像均采用此种方法获得。

步骤2、建立特征模型

本发明拟使用三类特征构成描述结节的特征向量。(1)结节的影像学特征，例如分叶征，毛刺征，胸膜凹陷征等；(2)图像的形状及纹理特征，例如傅里叶描述子、灰度共生矩阵等；(3)患者的相关信息，例如患者年龄，性别等。具体特征提取方案，如图6所示，首先，二维灰度图像是肺结节的最基本信息，基于结节的二维图像，本发明拟建立灰度共生矩阵统计量特征模型等；其次，本发明拟对二维灰度图形进行gamma校正得到伽马校正图像，建立毛刺特征模型；然后，将二维灰度图像进行动态规划分割得到肺结节边界，根据肺结节的边界，本发明拟建立空泡特征模型等；接着，将肺结节边界与二维灰度图像进行叠加，基于此图像，本发明拟建立胸膜凹陷特征模型等；在肺结节边界上提取特征点，基于提取的特征点，本发明拟建立分叶特征模型；最后，结合患者信息构成总的特征向量。

2.1建立结节影像学特征模型

(1)建立分叶特征模型

分叶征是指肿块的轮廓并非纯粹的圆形或椭圆，表面常呈凹凸不平的多个弧形，在二维图像中结节边缘明显高低不平的弧形。本专利采用弦长与弦距比、特征点个数、最小凹角、可疑分叶点数、最大凹弧线比、分叶等级、最小凸角、最大凸弧线比来描述分叶征。

如图4所示，特征点检测：第一步，计算各边界点斜率。根据二维图像分割边界，求各边界点与其相邻后边的边界点的斜率。第二步，判断特征点。相邻两边界点斜率符号为正负则前点为凸；为负正则为凹；否则向下检测。将特征点依次放入集合A。第三步，优化特征点集合。检测相邻凹点之间的距离小于一定阈值且中间凸点也小于一定阈值，满足两条件则将三个特征点中的后两个从集合A中删除。

由图4过程得到特征点集，即特征点集A，a表示A中的特征点，a_i表示A中的第i个特征点。特征点集A中三个相邻的特征点分别为a_i-1，a_i，a_i+1。a_i-1、a_i+1具有相同的凹凸性，与a_i相反。a_i-1与a_i+1之间的距离为弦长d_ci。a_i到直线a_i-1a_i+1的距离为弦距d_vi。其中，d表示距离，下角标c表示弦长，下角标v表示弦距。弦长与弦距比为

其中，R表示弦长与弦距比，R_i表示a_i对应的弦长与弦距比。以a_i为顶点所形成的夹角度数为α_i。

●计算平均凸弧线比

R_m=mean{R_i|a_i为凸点}    (1)

其中，R为弦长与弦距比，下角标m表示平均，R_m表示平均凸弧线比，a表示A中的特征点，a_i为特征点集A中第i个特征点。

●统计特征点个数

N_f=card A

其中，A为特征点集合，N表示个数，下角标f标识为特征点，N_f表示特征点个数。

●计算最小凹角与最小凸角

相邻三个特征点a_i-1，a_i，a_i+1对应的坐标分别为(x_i-1,y_i-1)，(x_i,y_i)，(x_i+1,y_i+1)。通过余弦定理可求得a_i为顶点所形成的夹角α_i。

$\cos {\mathrm{\&alpha;}}_i=\frac{{(x_{i+1}-x_i)}^2+{(y_{i+1}-y_i)}^2+{(x_{i-1}-x_i)}^2+{(y_{i-1}-y_i)}^2-({(x_{i+1}-x_{i-1})}^2+{(y_{i+1}-y_{i-1})}^2)}{2\sqrt{{(x_{i+1}-x_i)}^2+{(y_{i+1}-y_i)}^2}\sqrt{{(x_{i-1}-x_i)}^2+{(y_{i-1}-y_i)}^2}}$

最小凹角大小：α_c=min{α_i|a_i为凹点}

最小凸角大小：α_t=min{α_i|a_i为凸点}

其中，a_i为特征点集A中第i个特征点，x_i，y_i分别表示特征点a_i在二维图像上的横坐标和纵坐标，α表示角度，α_i表示以a_i为顶点所形成的夹角度数，下角标c表示凹角，t表示凸角，α_c表示最小凹角，α_t表示最小凸角。

●统计可疑分叶点数

N_s=card{R_i>0.5}

其中，R_i表示a_i对应的弦长与弦距比，N表示个数，下角标s标识为可疑的分叶，N_s表示特征点个数。

●计算最大凹弧线比与最大凸弧线比

最大凹弧线比：R_c=max{R_i|a_i为凹点}

最大凸弧线比：R_t=max{R_i|a_i为凸点}

其中，a_i为特征点集A中第i个特征点，R_i表示a_i对应的弦长与弦距比，R为弦长与弦距比，下角标c表示凹角，t表示凸角，R_c表示最大凹弧线比，R_t表示最大凸弧线比。

(2)、建立毛刺特征模型

肺结节的毛刺征表现为肺结节边缘向周围伸展的、放射状的、无分支的、直而有力的细短线条影，本专利采用在二维图像中用方向梯度直方图(Histogram of Oriented Gradient，HOG)来描述毛刺征。

HOG的核心思想是所检测的局部物体外形能够被光强梯度或边缘方向的分布所描述。通过将整幅图像分割成小的连接区域(称为cells)，每个cell生成一个方向梯度直方图或者cell中pixel的边缘方向，这些直方图的组合可表示出(所检测目标的目标)描述子，毛刺特征建模流程如图7所示，具体为：

①对二维图像进行gamma校正；

②求图像竖直边界，水平边界，边缘强度，边缘斜率；

③将图像每16×16个像素分到一个cell中，对于大小为32×412的二维灰度图像，可以分成2×25个cell；

④对于每个cell求其梯度方向直方图。通常取9个方向，也就是每360/9=40度分到一个方向，方向大小按像素边缘强度加权。最后归一化直方图。

⑤每2*2个cell合成一个block，所以这里就有(2-1)*(25-1)=24个block。

⑥每个block中都有2*2*9个特征，一共有24个block，所以总的特征有24*36个。

(3)、建立胸膜凹陷特征模型

胸膜凹陷征表现为肺内结节与邻近胸膜之间的一条或几条线状影。在结节的二维图像中，如果结节存在胸膜凹陷，则相应列上会出现上下两组边界，如图3所示。本项目中，当出现胸膜凹陷征时，拟使用动态规划算法搜索2D图像的上边界如图5所示。利用上边界所含区域的大小PC来衡量胸膜凹陷的程度。

其中，d为结节中心到肺边界的距离，d_max为结节的最大直径。如果

则该结节可能为与胸膜粘连的结节；如果d在[d_max/2,2d_max]范围内，则该结节可能存在胸膜凹陷征，用搜索边界上部区域面积PC来度量胸膜凹陷程度；如果d>2d_max，则该结节远离肺壁，不与胸膜粘连。

(4）、建立空泡特征模型

空泡征表现为肺结节内的类圆形含气腔，在肺窗上呈低密度影，使用变量空泡体积百分比PE_va来表示

${\mathrm{PE}}_{\mathrm{va}}=\frac{V_{\mathrm{va}}}{V}---\left(3\right)$

其中，V为肺结节的体积，V_va为肺结节内空泡体积，由公式(8)计算得到。

在二维图像中，取结节边界以内的部分为感兴趣区域，对感兴趣区域进行阈值分割得到空泡边界，然后根据空泡边界求得空泡体积。

(5)建立边界清晰度特征模型

边界清晰程度是描述肺结节边界清晰-模糊程度的参数。在临床应用中，特别是磨玻璃样结节的良恶性鉴别中，发挥重要的作用。在本专利中，拟利用肺结节边界点的径向梯度和GR来表示结节边界清晰程度。如果径向梯度和越大，表明结节的边界越清晰。

$\mathrm{GR}=\frac{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}(f(j,\mathrm{boundary}\left(j\right)+1)-f(j,\mathrm{boundary}(j))}{N}---\left(4\right)$

其中，f为结节的二维灰度图像，j为边界中第几个边界点，boundary(j)为结节边界第j个边界点对应的纵坐标，N为边界点的数量。

(6)建立磨玻璃样特征模型

磨玻璃样在CT图像上表现为界限较清楚、密度较淡影。本专利中，根据边界清晰程度与平均密度来度量结节的磨玻璃样。

$G_{\mathrm{go}}={\mathrm{\&omega;}}_1\frac{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}(f(j,\mathrm{boundary}\left(j\right)+1)-f(j,\mathrm{boundary}\left(j\right))}{N}+{\mathrm{\&omega;}}_2\frac{\underset{(x,y)\&Element;\mathrm{RE}}{\mathrm{\&Sigma;}}f(x,y)}{N_{\mathrm{RE}}}---\left(5\right)$

其中，f为结节的二维灰度图像，boundary(j)为结节边界第j个边界点对应的纵坐标，N为边界点的数量，RE为二维图像中结节分割边界以内的区域，N_RE为RE区域像素点个数，ω₁为边界清晰程度调节系数，ω₂为平均密度的调节系数，ω₁，ω₂使二者处于相近的数量级。

(7)建立钙化含量征特征模型

肿瘤内部由于坏死区域营养不良或者肿瘤异位内分泌导致钙盐沉着，形成钙化。在二维图像中表现为结节中有一个比较亮的区域。钙化的分布、形态及含量对分辨肿瘤的良恶性也很重要。恶性病灶的钙化形态多表现为偏心性、针尖样，少数为中心性，且含量不足10％；良性病灶的钙化形态多表现为层状、中心性、爆米花样及多发弥漫性，且含量超过10％。本专利采用钙化含量百分比PR_ca来量化结节的钙化程度。

${\mathrm{PR}}_{\mathrm{ca}}=\frac{V_{\mathrm{ca}}}{V}\&times;100\%---\left(6\right)$

其中，V为肺结节的体积，V_ca为肺结节内钙化部分体积，由公式(8)计算得到。

在二维图像中，取结节边界以内的部分为感兴趣区域，对感兴趣区域进行动态规划分割得到钙化边界，然后根据钙化边界求得钙化体积。

2.3建立图像形状及纹理的特征模型

(1)、获取结节表面积

在本专利中，由于螺旋扫描技术产生的视点近似均匀的分布在结节表面。因此，在结节分割后的图像中，距结节中心越远表面积越大。距离结节中心为r个像素的视点所在区域表面积是距结节中心1个像素视点所在区域表面积的r²倍。

$S=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}s\&times;{\left(\mathrm{boundary}\left(j\right)\right)}^2---\left(7\right)$

其中，N为结节边界点的数量，s为距结节中心1个像素视点所在表面积

boundary(j)为结节边界第j个边界点对应的纵坐标，由二维图像的生成原理可知，二维图像中结节边界点的纵坐标也表示边界点距离结节中心的距离。

(2)获取结节体积

肺结节的体积常用于估计肺结节的倍增时间，它是医生判断肺结节良、恶性的重要指标。在本专利中，由于螺旋扫描技术产生的视点近似均匀的分布在结节表面。因此，在结节分割后的图像中，距结节中心越远体积越大。距离结节中心为r个像素的视点所在体积是距结节中心1个像素视点所在体积的r³倍。

$V=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}v\&times;{\left(\mathrm{boundary}\left(j\right)\right)}^3---\left(8\right)$

其中，N为结节边界点的数量，v为距结节中心1个像素视点所在体积

boundary(j)为二维结节边界第j个边界点对应的纵坐标，由二维图像的生成原理可知，二维图像中结节边界点的纵坐标也表示边界点距离结节中心的距离。

(3、)获取结节球形度

球形度是描述结节近似球形的程度。它在[0，1]范围内，当球形度越接近1，则说明该结节越近似于球形。

$D=\frac{\sqrt[6]{36\mathrm{\&pi;}}\sqrt[3]{V}}{\sqrt{S}}---\left(9\right)$

其中，V为肺结节的体积，由公式(8)计算可得，S为肺结节的表面积，由公式(7)计算可得。

(4)、获取结节最大直径

最大直径d_max是描述结节大小的一个重要特征，在临床上有着重要的应用价值。在本专利中，由于螺旋扫描技术产生的视点近视对称的分布在结节表面。因此，在结节分割后的图像中，过结节中心，且过表面上下、左右对称视点的射线长度为最大直径d_max。

(5)、获取结节平均密度

$M_n=\frac{\underset{(x,y)\&Element;\mathrm{RE}}{\mathrm{\&Sigma;}}f(x,y)}{N_{\mathrm{RE}}}---\left(10\right)$

其中，f(x,y)为二维图像位于(x,y)点的灰度值，RE为二维图像中结节分割边界以内的区域，N_RE为RE区域像素点个数，M_n表示结节平均密度。

(6、)获取结节密度方差

$V_n=\sqrt{\frac{1}{N_{\mathrm{RE}}}\underset{(x,y)\&Element;\mathrm{RE}}{\mathrm{\&Sigma;}}{(f(x,y)-M_n)}^2}---\left(11\right)$

其中，f(x,y)为二维图像位于(x,y)点的灰度值，RE为二维图像中结节分割边界以内的区域，N_RE为RE区域像素点个数，M_n表示结节平均密度，V_n表示结节密度方差。

(7)、获取结节傅里叶描述子FD

傅里叶描述子是物体边界信号的频域分析的结果，它也是物体形状边界曲线的傅里叶变换系数，实验表明，傅里叶描述子相较于各种形状识别方法而言都具有最佳的形状识别性能。在本专利中，结节边界是在二维图像上通过动态规划获得。

$\mathrm{FD}\left(u\right)=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\mathrm{boundary}\left(j\right)e^{-j^{\&prime;}2\mathrm{\&pi;u}(j-1)/N}---\left(12\right)$

其中，FD(u)为结节傅里叶描述子，u=0,1,2,…,N-1，N为结节边界点的数量，boundary(j)为结节边界第j个边界点对应的纵坐标，j'表示复系数，由二维图像的生成原理可知，二维图像中结节边界点的纵坐标也表示边界点距离结节中心的距离。

(8)获取平均标准化半径

提取平均标准化半径之前需要对半径进行归一化，

$d_{\mathrm{mr}}=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_{\mathrm{no}}\left(j\right)---\left(13\right)$

其中，d_mr为平均标准化半径，d_no(j)为结节边界点距离结节中心的归一化距离，N为结节边界点的数量。

(9)获取标准化半径方差

$\mathrm{\&sigma;}=\sqrt{\frac{1}{N-1}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{(d_{\mathrm{no}}\left(j\right)-d_{\mathrm{mr}})}^2}---\left(14\right)$

其中，d_mr为平均标准化半径，d_no(j)为结节边界点距离结节中心的归一化距离，N为结节边界点的数量。

(10)获取标准化半径的熵E(entropy of the normalized radial)

$E=-\underset{k=1}{\overset{100}{\mathrm{\&Sigma;}}}{p}_k\left(\log \left(p_k\right)\right)---\left(15\right)$

其中，式中的p_k表示概率值，这个概率值是通过标准直方图来进行计算得出的。标准直方图的计算方法是将最大半径与最小半径之间的差值平均划分为100等分，然后将所有的半径分别装入自己所属的区段，最后计算出每个区段所占的比例，即p_k，它表示在第k个区段内的半径数所占的比例。

(11)、获取面积比率A_r

$A_r=\frac{1}{d_{\mathrm{mr}}N}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}(d_{\mathrm{no}}\left(j\right)-d_{\mathrm{mr}})---\left(16\right)$

其中，d_no(j)-d_mr=0，

d_mr为平均标准化半径，d_no(j)为结节边界点距离结节中心的归一化距离，N为结节边界点的数量。

(12)、获取粗糙度R(boundary roughness)

$\mathrm{Rough}=\frac{1}{N}\underset{j=1}{\overset{N-1}{\mathrm{\&Sigma;}}}|d_{\mathrm{no}}\left(j\right)-d_{\mathrm{no}}(j+1)|---\left(17\right)$

其中，d_no(j)为结节边界点距离结节中心的归一化距离，N为结节边界点的数量。

(13)、获取灰度共生矩阵统计量

共生矩阵用两个位置的像素的联合概率密度来定义，它不仅反映亮度的分布特性，也反映具有同样亮度或接近亮度的像素之间的位置分布特性，是有关图像亮度变化的二阶统计特征。它是定义一组纹理特征的基础。

一幅图像的灰度共生矩阵能反映出图像灰度关于方向、相邻间隔、变化幅度的综合信息，它是分析图像的局部模式和它们排列规则的基础。本专利采用灰度共生矩阵统计量对结节的二维图像进行分析，二维图像采用螺旋扫描技术生成。在二维图像中，水平相邻的两像素点距离结节中心越近，在三维结节中两像素点间的实际距离越近。

设f(x,y)为一幅二维数字图象，其大小为M×N，灰度级别为Ng，则满足一定空间关系的灰度共生矩阵为：

P(m,n)=#{(x1,y1),(x2,y2)∈M×N|f(x1,y1)=m,f(x2,y2)=n}    (18)

其中#(x)表示集合X中的元素个数，显然P为Ng×Ng的矩阵，(x1,y1)、(x2,y2)为图像中的两点，对应的灰度级分别为m，n。若(x1,y1)与(x2,y2)间距离为dist，两者与坐标横轴的夹角为θ，则可以得到各种间距及角度的灰度共生矩阵P(m,n,dist,θ)。

在提取灰度共生矩阵的特征参数之前，要作正规化处理。令p(m,n,dist,θ)=P(m,n,dist,θ)/H，其中H是正规化常数，是灰度共生矩阵中全部元素之和。

Haralick等人定义的14个用于纹理分析的灰度共生矩阵特征参数如下：

●获取角二阶矩(Angular second moment)

$\mathrm{ASM}=\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}{p}^2(m,n,\mathrm{dist},\mathrm{\&theta;})---\left(19\right)$

其中，式中q为选取的灰度级Ng的大小。

角二阶矩也称能量，是灰度共生矩阵元素值的平方和，反映了图像灰度分布均匀程度和纹理粗细度。如果共生矩阵的所有值均相等，则ASM值小；相反，如果其中一些值大而其余值小，则ASM值大。当共生矩阵中元素集中分布时，此时ASM值大。ASM值大表明一种较均一和规则变化的纹理模式。

●获取对比度(Contrast)

$\mathrm{Constrast}=\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}\left[{(m-n)}^2{p}^2(m,n,\mathrm{dist},\mathrm{\&theta;})\right]---\left(20\right)$

其中，式中q为选取的灰度级Ng的大小。

对比度反映了图像的清晰度和纹理沟纹深浅的程度。纹理沟纹越深，其对比度越大，视觉效果越清晰；反之，对比度小，则沟纹浅，效果模糊。灰度差即对比度大的象素对越多，这个值越大。灰度共生矩阵中远离对角线的元素值越大，Constrast越大。

●获取相关性(Correlation)

$\mathrm{Correlation}=\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}(m\&times;n\&times;p(m,n,\mathrm{dist},\mathrm{\&theta;})-u_1\&times;u_2)/d_1{\&times;d}_2---\left(21\right)$

其中，式中q为选取的灰度级Ng的大小

$\begin{array}{cc}{u}_1=\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}m\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}p(m,n,\mathrm{dist},\mathrm{\&theta;})& u_2=\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}j\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}p(m,n,\mathrm{dist},\mathrm{\&theta;})\end{array}$

$\begin{array}{cc}{d_1}^2=\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}{(m-u_1)}^2\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}p(m,n,\mathrm{dist},\mathrm{\&theta;})& {d_2}^2=\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}{(n-u_2)}^2\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}p(m,n,\mathrm{dist},\mathrm{\&theta;})\end{array}$

它度量空间灰度共生矩阵元素在行或列方向上的相似程度，因此，相关值大小反映了图像中局部灰度相关性。当矩阵元素值均匀相等时，相关值就大;相反，如果矩阵像元值相差很大则相关值小。如果图像中有水平方向纹理，则水平方向矩阵的Correlation大于其余矩阵的Correlation值。

●获取熵(Entropy)

$\mathrm{Entropy}=-\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}p(m,n,\mathrm{dist},\mathrm{\&theta;})\&times;\mathrm{lgp}(m,n,\mathrm{dist},\mathrm{\&theta;})---\left(22\right)$

其中，式中q为选取的灰度级Ng的大小。

熵是图像具有的信息量的度量，纹理信息也属于图像的信息，是一个随机性的度量，当共生矩阵中所有元素有最大的随机性、空间共生矩阵中所有值几乎相等时，共生矩阵中元素分散分布时，熵较大。它表示了图像中纹理的非均匀程度或复杂程度。

●获取方差(Variance)

$\mathrm{Variance}=\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}{(m-\mathrm{mean})}^{2}p(m,n,\mathrm{dist},\mathrm{\&theta;})---\left(23\right)$

其中，式中q为选取的灰度级Ng的大小，mean为p(m,n,dist,θ)的均值。

●获取方差和(Sum of variance)

$\mathrm{SOV}=\underset{l=2}{\overset{2q}{\mathrm{\&Sigma;}}}{(l-\mathrm{SOA})}^2{p}_x\left(l\right)---\left(24\right)$

其中，式中q为选取的灰度级Ng的大小，SOA由公式(25)计算得到，

$p_x\left(l\right)=\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}p{(m,n,\mathrm{dist},\mathrm{\&theta;})}_{|m+n|=l}l=\mathrm{2,3},\&CenterDot;\&CenterDot;\&CenterDot;,2q.$

●获取均值和(Sum of average)

$\mathrm{SOA}=\underset{l=2}{\overset{2q}{\mathrm{\&Sigma;}}}{l\&times;p}_x\left(l\right)---\left(25\right)$

其中，式中q为选取的灰度级Ng的大小，

$p_x\left(l\right)=\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}p{(m,n,\mathrm{dist},\mathrm{\&theta;})}_{|m+n|=l}l=\mathrm{2,3},\&CenterDot;\&CenterDot;\&CenterDot;,2q.$

均值和是图像区域内像素点平均灰度值的度量，反映了图像的明暗深浅，适用于灰度图像。

●获取逆差矩(Inverse difference moment)

$\mathrm{IDM}=\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}p(m,n,\mathrm{dist},\mathrm{\&theta;})/(1+{(m-2)}^2)---\left(26\right)$

其中，式中q为选取的灰度级Ng的大小。

逆差矩又称为局部平稳，它是图像纹理局部变化的度量，反映了纹理的规则程度。纹理越规则，IDM就越大；反之亦然。

●获取差的方差(Variance of difference)

$\mathrm{VOD}=\underset{l=0}{\overset{q-1}{\mathrm{\&Sigma;}}}{[l-\underset{l=0}{\overset{q-1}{\mathrm{\&Sigma;}}}l\&times;p_y\left(l\right)]}^2\&times;p_y\left(l\right)---\left(27\right)$

其中，式中q为选取的灰度级Ng的大小，

$p_y\left(l\right)=\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}p{(m,n,\mathrm{dist},\mathrm{\&theta;})}_{|m+n|=l}l=\mathrm{0,1},\&CenterDot;\&CenterDot;\&CenterDot;,q-1.$

差的方差是邻近像素对灰度值差异的方差，对比越强烈，VOD越大；反之亦然。

●获取和熵(Sum of entropy)

$\mathrm{SOE}=-\underset{l=2}{\overset{2q}{\mathrm{\&Sigma;}}}{p}_x\left(l\right)\&times;\log \left[p_x\left(l\right)\right]---\left(28\right)$

其中，式中q为选取的灰度级Ng的大小，

$p_x\left(l\right)=\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}p{(m,n,\mathrm{dist},\mathrm{\&theta;})}_{|m+n|=l}l=\mathrm{2,3},\&CenterDot;\&CenterDot;\&CenterDot;,2q.$

●获取差熵(Difference of entropy)

$\mathrm{DOE}=-\underset{l=0}{\overset{q-1}{\mathrm{\&Sigma;}}}{p}_y\left(l\right)\&times;\log \left[p_y\left(l\right)\right]---\left(29\right)$

其中，式中q为选取的灰度级Ng的大小，

$p_y\left(l\right)=\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}p{(m,n,\mathrm{dist},\mathrm{\&theta;})}_{|m+n|=l}l=\mathrm{0,1},\&CenterDot;\&CenterDot;\&CenterDot;,q-1.$

●获取聚类阴影(Clustering shadow)

$\mathrm{CS}=-\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}{[(m-u_1)+(n-u_2)]}^{3}p(m,n,\mathrm{dist},\mathrm{\&theta;})---\left(30\right)$

其中，式中q为选取的灰度级Ng的大小，

$\begin{array}{cc}{u}_1=\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}m\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}p(m,n,\mathrm{dist},\mathrm{\&theta;})& u_2=\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}n\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}p(m,n,\mathrm{dist},\mathrm{\&theta;})\end{array}.$

●获取显著聚类(Significant clustering)

$\mathrm{SC}=-\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}{[(m-u_1)+(n-u_2)]}^{4}p(m,n,\mathrm{dist},\mathrm{\&theta;})---\left(31\right)$

其中，式中q为选取的灰度级Ng的大小，

$\begin{array}{cc}{u}_1=\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}m\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}p(m,n,\mathrm{dist},\mathrm{\&theta;})& u_2=\underset{m=0}{\overset{q-1}{\mathrm{\&Sigma;}}}n\underset{n=0}{\overset{q-1}{\mathrm{\&Sigma;}}}p(m,n,\mathrm{dist},\mathrm{\&theta;})\end{array}.$

●获取最大概率(Maximum probability)

$\mathrm{MP}=\underset{m,n}{\mathrm{MAX}}\left(p(m,n,\mathrm{dist},\mathrm{\&theta;})\right)---\left(32\right)$

2.4结合患者的相关信息

(1)、患者年龄；

(2)、患者性别；

(3)、患者是否患过其他癌症；

(4)、患者是否痰中咯血。

步骤3、构建模糊分类器

虽然现在用于肺结节良恶性鉴别分类器有很多，但是这些分类器给出的结果都是确切的，即要么判断为恶性，要么判断为良性。在临床实践中，我们有时并不能给出确切的结论，只能给出一种结论的可能性有多大。所以在本专利中，我们将使用一种能够给出模糊结果的分类器，即模糊C-均值聚类分类器，并用马氏距离代替欧氏距离，克服欧氏距离不考虑量纲的问题。

模糊C-均值聚类(Fuzzy C-Means，FCM)算法以最小类内平方误差和为聚类准则，与C-均值聚类算法不同之处在于不将样本分成明确的子集，而是计算每个样本属于各模糊子集的隶属度。采用迭代法优化目标函数，获得对数据集的模糊分类。在数学上可以表示为目标函数求极值问题，即求下式中目标函数J的最小值：

$J(U,Z)=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{C}{\mathrm{\&Sigma;}}}{\left({\mathrm{\&mu;}}_{\mathrm{ij}}\right)}^m{d}^2(x_j,z_i)---\left(33\right)$

因为要将候选结节分成两类，所以类别数目C值为2。U={μ_ij}是C×N维模糊分类矩阵。m是一个常数，用于控制分类的隶属度。d²(x_j,z_i)表示样本x_j到第i类中心z_i的欧式距离，表示为d²(x_j,z_i)=||x_j-z_i||²。Z={z₁,z₂,...,z_c}是C维聚类中心矩阵，表示为

1≤i≤C。μ_ij表示样本x_j属于聚类中心z_i的隶属度，表示为

在分类过程中，不断迭代来更新隶属度函数和聚类中心。当满足迭代终止条件(1)迭代次数小于100，或(2)目标代价函数值J小于1.2×10⁵，或(3)连续两次迭代产生的聚类中心的差值小于0.1时，停止迭代取得最终的聚类中心z_i；

传统的模糊聚类算法使用欧氏距离作为相似度度量。欧氏距离在计算时，不考虑样本的量纲，导致变化大的数据和变化小的数据一起计算。变化大的数据对计算出的距离贡献大，容易产生偏差，而马氏距离正好可以克服这个问题，其计算不受特征单位的影响。

设s维空间的样本集X={x₁,x₂,…,x_n}，样本x到样本集X的马氏距离如公式(34)所示。

$D_M\left(x\right)=\sqrt{{(x-\mathrm{\&mu;})}^T{\mathrm{\&Sigma;}}^{-1}(x-\mathrm{\&mu;})}---\left(34\right)$

式中μ--样本集X的均值向量;∑--为协方差矩阵，

在计算马氏距离的过程中，若协方差矩阵为奇异矩阵，则协方差矩阵的逆阵不存在，将马氏距离无法直接求得。根据矩阵理论，任一秩为r的实对称半正定矩阵∑都可按公式(33)分解。

∑=A^TG^-1A    (35)

式中G--为r×r的对角阵，它由∑的非0特征值构成，G是非奇异的；

A--为r×m矩阵，由G中特征值所对应的特征向量构成，它张成了∑的非退化子空间，且AA^T为r×r的单位阵。

通过对公式(35)的分解，就可根据G的逆来求∑的伪逆矩阵，如公式(36)所示。

∑^-1=A^TG^-1A    (36)

由于马氏距离可以克服欧氏距离的缺点，可以用马氏距离来代替欧氏距离应用于模糊C均值聚类算法中。

将公式(34)代入到传统模糊C均值聚类算法的全局代价函数(33)中，得到基于马氏距离的模糊C均值聚类算法的全局代价函数，如公式(37)所示，约束条件如公式(38)、(39)所示：

$J_{\mathrm{FCM}}(U,Z)=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{C}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^m{d}^2(x_j-z_i)=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{C}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^m{(x_j-z_i)}^T{\mathrm{\&Sigma;}}^{-1}(x_j-z_i)---\left(37\right)$

$\begin{array}{cc}0\&le;{\mathrm{\&mu;}}_{\mathrm{ij}}\&le;1& \&ForAll;i\&Element;\end{array}\{\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,c\},j\&Element;\{\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n\}---\left(38\right)$

$\begin{array}{cc}\underset{i=1}{\overset{C}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}=1& \&ForAll;j\&Element;\{\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n\}\end{array}---\left(39\right)$

依据GK算法，聚类中心、协方差矩阵和隶属度值的计算公式分别如公式(40)、(41)、(42)所示。

$z_i=\frac{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^m{x}_j}{\underset{i=1}{\overset{C}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^m},1\&le;i\&le;C---\left(40\right)$

${\mathrm{\&Sigma;}}_i=\frac{\underset{j=1}{\overset{l}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^m(x_j-p_i){(x_j-p_i)}^T}{\underset{j=1}{\overset{l}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^m},1\&le;i\&le;C---\left(41\right)$

${\mathrm{\&mu;}}_{\mathrm{ij}}=\frac{D_{\mathrm{ij}}^{-2/(m-1)}}{\underset{k=1}{\overset{C}{\mathrm{\&Sigma;}}}{D}_{\mathrm{kj}}^{-2/(m-1)}},1\&le;i\&le;C,1\&le;j\&le;l---\left(42\right)$

式中D_ij--为马氏距离 $D_{\mathrm{ij}}=\sqrt{{(x_j-z_i)}^T{\mathrm{\&Sigma;}}_I^{-1}(x_j-z_i)}.$

在模糊C均值聚类算法中，模糊参数m被视为平滑因子，控制着样本隶属类别的模糊程度，其选择影响到模糊C均值聚类的性能，因此一个合适的模糊参数m值对模糊聚类的性能起着至关重要的作用。通过总结前人在模糊参数m方面做的研究，在本聚类算法中，m值取2。

模糊C均值聚类算法的算法流程：(如图8)

(1)设置目标函数精度ε，模糊指数m(m通常为2)，最大迭代次数T_m；

(2)初始化模糊聚类中心；

(3)由步骤(3)-(5)更新模糊划分矩阵U={μ_ij}和聚类中心Z={z₁,z₂,...,z_c}；

(4)若|J(t)-J(t-1)|<ε或t>T_m，则结束聚类；t←t+1并转到(3)步；

(5)由U={μ_ij}得到各像素点分类结果。

步骤4、评价分类性能

目前诊断实验的准确性评估指标有很多，包括灵敏度、特异度、假阴性率、假阳性率等评价指标。其中灵敏度(真阳性率，Sensitivity)和特异性(真阴性率，Specificity)分别反映了真实情况为有病时诊断系统发现疾病的能力，与真实情况为无病时诊断系统排除疾病的能力。为避免发病率对评估结果的影响，通常综合两者进行评估。但其缺点是随着诊断界值的改变，灵敏度和特异性会发生相应的同步变化，即当灵敏度提高时，特异性被降低；而当灵敏度降低时，特异性反而会被提高。本发明采用的样本数据来自肺部图像数据库联盟(LIDC)，同时采用敏感性、准确性、特异性和假阳性四个指标来分别评价分类性能。

表1：诊断实验评价形式

(5)根据敏感性(Sensitivity)(又称为真阳性率)，判断对结节样本恶性的检测性能。

$\mathrm{Sensitivity}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}---\left(43\right)$

其中，TP为真阳性结节的数量，FP为假阳性结节的数量。

(6)根据准确性(Accuracy)，判断所有的检测结果中，检测恶性的样本占全部样本的比率。

$\mathrm{Accuracy}=\frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{FP}+\mathrm{TN}+\mathrm{FN}}---\left(44\right)$

其中，TP为真阳性结节的数量，FP为假阳性结节的数量，TN为假阴性结节的数量，FN为真阴性结节的数量。

(7)根据特异性(Specificity)，判断对结节样本良性的检测性能。

$\mathrm{Specificity}=\frac{\mathrm{TN}}{\mathrm{TN}+\mathrm{FP}}---\left(45\right)$

其中，FP为假阳性结节的数量，TN为假阴性结节的数量。

(8)根据假阳性率，判断检测结果中被错误检测为良性结节样本所占的比率。

$\mathrm{FPF}=\frac{\mathrm{FP}}{\mathrm{TN}+\mathrm{FP}}---\left(46\right)$

其中，TP为真阳性结节的数量，FP为假阳性结节的数量，TN为假阴性结节的数量。

Claims (9)

1.一种基于多维信息的肺结节良恶性鉴别方法，其特征在于：包括如下步骤： 

步骤1、三维结节的二维表示； 

步骤2、建立特征模型； 

步骤3、构建模糊分类器； 

步骤4、评价分类性能。 

2.根据权利要求1所述的基于多维信息的肺结节良恶性鉴别方法，其特征在于：所述的步骤1、三维结节的二维表示具体为：采用螺旋扫描技术生成在候选结节表面上有序均匀分布的视点，在生成有序、均匀分布视点之后，以候选结节的中心为起始点向视点方向均匀的指定长度的射线，按顺时针方向将这些射线上的像素垂直地排列得到三维候选结节的二维图像。 

3.根据权利要求1所述的基于多维信息的肺结节良恶性鉴别方法，其特征在于：所述的步骤2、建立特征模型建立过程包括如下步骤： 

步骤2.2建立结节影像学特征模型； 

步骤2.3建立图像形状及纹理的特征模型； 

步骤2.4结合患者的相关信息。 

4.根据权利要求3所述的基于多维信息的肺结节良恶性鉴别方法，其特征在于：所述的步骤2.2建立结节影像学特征模型；包括如下步骤： 

(1)建立分叶特征模型 

采用弦长与弦距比、特征点个数、最小凹角、可疑分叶点数、最大凹弧线比、分叶等级、最小凸角、最大凸弧线比来描述分叶征； 

首先特征点检测，方法如下：a，根据二维图像分割边界，求各边界点与其相邻后边的边界点的斜率，获取各边界点斜率；b，判断特征点，具体为相邻两边界点斜率符号为正负则前点为凸；为负正则为凹；否则向下检测；将特征点依次放入集合A；c，优化特征点集合，检测相邻凹点之间的距离小于一定阈值且中间凸点也小于一定阈值，满足两条件则将三个特征点中的后两个从集合A中删除； 

a表示A中的特征点，a_i表示A中的第i个特征点，特征点集A中三个相邻的特征点分别为a_i-1，a_i，a_i+1；a_i-1、a_i+1具有相同的凹凸性，与a_i相反；a_i-1与a_i+1之间的距离为弦长d_ci；a_i到直线a_i-1a_i+1的距离为弦距d_vi；其中，d表示距离，下角标c表示弦长，下角标v表示弦距，弦长与弦距比为

其中，R表示弦长与弦距比，R_i表示a_i对应的弦长与弦距比， 以a_i为顶点所形成的夹角度数为α_i； 

●获取平均凸弧线比 

R_m=mean{R_i|a_i为凸点}    (1) 

其中，R为弦长与弦距比，下角标m表示平均，R_m表示平均凸弧线比，a表示A中的特征点，a_i为特征点集A中第i个特征点； 

●统计特征点个数 

N_f=card A 

其中，A为特征点集合，N表示个数，下角标f标识为特征点，N_f表示特征点个数； 

●获取最小凹角与最小凸角 

相邻三个特征点a_i-1，a_i，a_i+1对应的坐标分别为(x_i-1,y_i-1)，(x_i,y_i)，(x_i+1,y_i+1)；通过余弦定理可求得a_i为顶点所形成的夹角α_i； 

最小凹角大小：α_c=min{α_i|a_i为凹点} 

最小凸角大小：α_t=min{α_i|a_i为凸点} 

其中，a_i为特征点集A中第i个特征点，x_i，y_i分别表示特征点a_i在二维图像上的横坐标和纵坐标，α表示角度，α_i表示以a_i为顶点所形成的夹角度数，下角标c表示凹角，t表示凸角，α_c表示最小凹角，α_t表示最小凸角； 

●统计可疑分叶点数 

N_s=card{R_i>0.5} 

其中，R_i表示a_i对应的弦长与弦距比，N表示个数，下角标s标识为可疑的分叶，N_s表示特征点个数； 

●获取最大凹弧线比与最大凸弧线比 

最大凹弧线比：R_c=max{R_i|a_i为凹点} 

最大凸弧线比：R_t=max{R_i|a_i为凸点} 

其中，a_i为特征点集A中第i个特征点，R_i表示a_i对应的弦长与弦距比，R为弦长与弦距比，下角标c表示凹角，t表示凸角，R_c表示最大凹弧线比，R_t表示最大凸弧线比； 

(2)建立毛刺特征模型 

具体步骤如下： 

①对二维图像进行伽马校正； 

②利用卷积函数求图像竖直边缘Ix，水平边缘Iy，然后根据竖直边缘Ix和水平边缘Iy来求得边缘强度sqrt(Ix.^2＋Iy.^2)和边缘斜率Iy./Ix； 

③将图像每m′×n′个像素分到一个单元中，对于大小为M′×N′的二维灰度图像，可以分成至少两个单元，这里每个单元大小为16*16，图像大小为32*412； 

④对于每个单元求其梯度方向直方图，通常取Q个方向，也就是每360/Q分到一个方向，方向大小按像素边缘强度加权，最后归一化直方图； 

⑤每2*2个单元合成一个块，所以这里就有(2-1)*(25-1)=24个块； 

⑥每个块中都有2*2*9个特征，这里一共有24个块，所以24个36维的特征构成总特征； 

(3)建立胸膜凹陷特征模型 

胸膜凹陷征表现为肺内结节与邻近胸膜之间的一条或几条线状影，在结节的二维图像中，如果结节存在胸膜凹陷，则相应列上会出现上下两组边界，当出现胸膜凹陷征时，拟使用动态规划算法搜索2D图像的上边界，利用上边界所含区域的大小PC来衡量胸膜凹陷的程度； 

其中，d为结节中心到肺边界的距离，d_max为结节的最大直径，如果则该结节可能为与胸膜粘连的结节；如果d在[d_max/2,2d_max]范围内，则该结节可能存在胸膜凹陷征，用搜索边界上部区域面积PC来度量胸膜凹陷程度；如果d>2d_max，则该结节远离肺壁，不与胸膜粘连； 

(4)建立空泡特征模型 

空泡征表现为肺结节内的类圆形含气腔，在肺窗上呈低密度影，使用变量空泡体积百分 比PE_va来表示 

其中，V为肺结节的体积，V_va为肺结节内空泡体积，由公式(8)计算得到； 

在二维图像中，取结节边界以内的部分为感兴趣区域，对感兴趣区域进行阈值分割得到空泡边界，然后根据空泡边界求得空泡体积； 

(5)建立边界清晰度特征模型 

利用肺结节边界点的径向梯度和GR来表示结节边界清晰程度，如果径向梯度和越大，表明结节的边界越清晰， 

其中，f为结节的二维灰度图像，j为边界中第几个边界点，boundary(j)为结节边界第j个边界点对应的纵坐标，N为边界点的数量； 

(6)建立磨玻璃样特征模型 

根据边界清晰程度与平均密度来度量结节的磨玻璃样， 

其中，f为结节的二维灰度图像，boundary(j)为结节边界第j个边界点对应的纵坐标，N为边界点的数量，RE为二维图像中结节分割边界以内的区域，N_RE为RE区域像素点个数，ω₁为边界清晰程度调节系数，ω₂为平均密度的调节系数，ω₁,ω₂使二者处于相近的数量级； 

(7)建立钙化含量征特征模型 

采用钙化含量百分比PR_ca来量化结节的钙化程度， 

其中，V为肺结节的体积，V_ca为肺结节内钙化部分体积，由公式(8)计算得到； 

在二维图像中，取结节边界以内的部分为感兴趣区域，对感兴趣区域进行动态规划分割得到钙化边界，然后根据钙化边界求得钙化体积。 

5.根据权利要求3所述的基于多维信息的肺结节良恶性鉴别方法，其特征在于：所 述的步骤2.3建立图像形状及纹理的特征模型；包括如下步骤： 

(1)获取结节表面积 

由于螺旋扫描技术产生的视点近似均匀的分布在结节表面，在结节分割后的图像中，距结节中心越远表面积越大，距离结节中心为r个像素的视点所在区域表面积是距结节中心1个像素视点所在区域表面积的r²倍， 

其中，N为结节边界点的数量，s为距结节中心1个像素视点所在表面积

boundary(j)为结节边界第j个边界点对应的纵坐标，由二维图像的生成原理可知，二维图像中结节边界点的纵坐标也表示边界点距离结节中心的距离； 

(2)获取结节体积 

肺结节的体积常用于估计肺结节的倍增时间，它是医生判断肺结节良、恶性的重要指标，在本专利中，由于螺旋扫描技术产生的视点近似均匀的分布在结节表面，因此，在结节分割后的图像中，距结节中心越远体积越大，距离结节中心为r个像素的视点所在体积是距结节中心1个像素视点所在体积的r³倍， 

其中，N为结节边界点的数量，v为距结节中心1个像素视点所在体积

boundary(j)为二维结节边界第j个边界点对应的纵坐标，由二维图像的生成原理可知，二维图像中结节边界点的纵坐标也表示边界点距离结节中心的距离； 

(3)获取结节球形度 

球形度是描述结节近似球形的程度，在[0，1]范围内，当球形度越接近1，则说明该结节越近似于球形； 

其中，V为肺结节的体积，由公式(8)计算可得，S为肺结节的表面积，由公式(7)计算可 得； 

(4)获取结节最大直径 

最大直径d_max是描述结节大小的一个重要特征，由于螺旋扫描技术产生的视点近视对称的分布在结节表面，在结节分割后的图像中，过结节中心，且过表面上下、左右对称视点的射线长度为最大直径d_max； 

(5)获取结节平均密度 

其中，f(x,y)为二维图像位于(x,y)点的灰度值，RE为二维图像中结节分割边界以内的区域，N_RE为RE区域像素点个数，M_n表示结节平均密度； 

(6)获取结节密度方差 

其中，f(x,y)为二维图像位于(x,y)点的灰度值，RE为二维图像中结节分割边界以内的区域，N_RE为RE区域像素点个数，M_n表示结节平均密度，V_n表示结节密度方差； 

(7)获取结节傅里叶描述子 

结节边界是在二维图像上通过动态规划获得； 

其中，FD(u)为结节傅里叶描述子，u=0,1,2,…,N-1，N为结节边界点的数量，boundary(j)为结节边界第j个边界点对应的纵坐标，j'表示复系数，由二维图像的生成原理可知，二维图像中结节边界点的纵坐标也表示边界点距离结节中心的距离； 

(8)获取平均标准化半径 

提取平均标准化半径之前需要对半径进行归一化，

其中，d_mr为平均标准化半径，d_no(j)为结节边界点距离结节中心的归一化距离，N为结节边界点的数量； 

(9)获取标准化半径方差 

其中，d_mr为平均标准化半径，d_no(j)为结节边界点距离结节中心的归一化距离，N为结节边界点的数量； 

(10)获取标准化半径的熵E(entropy of the normalized radial) 

其中，式中的p_k表示概率值，这个概率值是通过标准直方图来得出的，标准直方图的获取方法是将最大半径与最小半径之间的差值平均划分为100等分，然后将所有的半径分别装入自己所属的区段，最后得到每个区段所占的比例，即p_k，它表示在第k个区段内的半径数所占的比例； 

(11)获取面积比率A_r

其中，d_no(j)-d_mr=0，

d_mr为平均标准化半径，d_no(j)为结节边界点距离结节中心的归一化距离，N为结节边界点的数量； 

(12)获取粗糙度R(boundary roughness) 

其中，d_no(j)为结节边界点距离结节中心的归一化距离，N为结节边界点的数量； 

(13)获取灰度共生矩阵统计量 

采用灰度共生矩阵统计量对结节的二维图像进行分析，二维图像采用螺旋扫描技术生成，在二维图像中，水平相邻的两像素点距离结节中心越近，在三维结节中两像素点间的实际距离越近， 

设f(x,y)为一幅二维数字图象，其大小为M×N，灰度级别为Ng，则满足一定空间关系的灰度共生矩阵为： 

P(m,n)=#{(x1,y1),(x2,y2)∈M×N|f(x1,y1)=m,f(x2,y2)=n}    (18) 

其中#(x)表示集合X中的元素个数，显然P为Ng×Ng的矩阵，(x1,y1)、(x2,y2)为图像中的两点，对应的灰度级分别为m，n，若(x1,y1)与(x2,y2)间距离为dist，两者与坐标横轴的夹角为θ，则可以得到各种间距及角度的灰度共生矩阵P(m,n,dist,θ)； 

在提取灰度共生矩阵的特征参数之前，要作正规化处理，令p(m,n,dist,θ)=P(m,n,dist,θ)/H，其中H是正规化常数，是灰度共生矩阵中全部元素之和； 

定义用于纹理分析的灰度共生矩阵特征参数如下： 

●获取角二阶矩 

其中，式中q为选取的灰度级Ng的大小； 

角二阶矩也称能量，是灰度共生矩阵元素值的平方和，反映了图像灰度分布均匀程度和纹理粗细度，如果共生矩阵的所有值均相等，则ASM值小；相反，如果其中一些值大而其余值小，则ASM值大，当共生矩阵中元素集中分布时，此时ASM值大，ASM值大表明一种较均一和规则变化的纹理模式； 

●获取对比度 

其中，式中q为选取的灰度级Ng的大小， 

对比度反映了图像的清晰度和纹理沟纹深浅的程度，纹理沟纹越深，其对比度越大，视觉效果越清晰；反之，对比度小，则沟纹浅，效果模糊，灰度差即对比度大的象素对越多，这个值越大，灰度共生矩阵中远离对角线的元素值越大，Constrast越大； 

●获取相关性 

其中，式中q为选取的灰度级Ng的大小； 

用于度量空间灰度共生矩阵元素在行或列方向上的相似程度，相关值大小反映了图像中局部灰度相关性，当矩阵元素值均匀相等时，相关值就大;相反，如果矩阵像元值相差很大则相关值小，如果图像中有水平方向纹理，则水平方向矩阵的Correlation大于其余矩阵的Correlation值； 

●获取熵 

其中，式中q为选取的灰度级Ng的大小； 

熵是图像具有的信息量的度量，纹理信息也属于图像的信息，是一个随机性的度量，当共生矩阵中所有元素有最大的随机性、空间共生矩阵中所有值几乎相等时，共生矩阵中元素分散分布时，熵较大，它表示了图像中纹理的非均匀程度或复杂程度； 

●获取方差 

其中，式中q为选取的灰度级Ng的大小，mean为p(m,n,dist,θ)的均值； 

●获取方差和 

其中，式中q为选取的灰度级Ng的大小，SOA由公式(25)计算得到， 

●获取均值和 

其中，式中q为选取的灰度级Ng的大小， 

均值和是图像区域内像素点平均灰度值的度量，反映了图像的明暗深浅，适用于灰度图像； 

●获取逆差矩 

其中，式中q为选取的灰度级Ng的大小； 

逆差矩又称为局部平稳，它是图像纹理局部变化的度量，反映了纹理的规则程度，纹理越规则，IDM就越大；反之亦然； 

●获取差的方差 

其中，式中q为选取的灰度级Ng的大小， 

差的方差是邻近像素对灰度值差异的方差，对比越强烈，VOD越大；反之亦然； 

●获取和熵 

其中，式中q为选取的灰度级Ng的大小， 

●获取差熵 

其中，式中q为选取的灰度级Ng的大小， 

●获取聚类阴影 

其中，式中q为选取的灰度级Ng的大小， 

●获取显著聚类 

其中，式中q为选取的灰度级Ng的大小， 

●获取最大概率 

6.根据权利要求3所述的基于多维信息的肺结节良恶性鉴别方法，其特征在于：所述的步骤2.4结合患者的相关信息，包括如下信息： 

(1)患者年龄； 

(2)患者性别； 

(3)患者是否患过其他癌症； 

(4)患者是否痰中咯血。 

7.根据权利要求1所述的基于多维信息的肺结节良恶性鉴别方法，其特征在于：所述的步骤3、构建模糊分类器；具体如下： 

通过模糊C-均值聚类分类器，并用马氏距离代替欧氏距离构建； 

求下式中目标函数J的最小值： 

因为要将候选结节分成两类，所以类别数目C值为2。U={μ_ij}是C×N维模糊分类矩阵。N表示样本数量。m是一个常数，用于控制分类的隶属度。d²(x_j,z_i)表示样本x_j到第i类中心z_i的欧式距离，表示为d²(x_j,z_i)=||x_j-z_i||²。Z={z₁,z₂,...,z_c}是C维聚类中心矩阵，表示为 

μ_ij表示样本x_j属于聚类中心z_i的隶属度，表示为

在分类过程中，不断迭代来更新隶属度函数和聚类中心。当满足迭代终止条件(1)迭代次数小于100，或(2)目标代价函数值J小于1.2×10⁵，或(3)连续两次迭代产生的聚类中心的差值小于0.1时，停止迭代取得最终的聚类中心z_i； 

设s维空间的样本集X={x₁,x₂,…,x_n}，样本x到样本集X的马氏距离如公式(34)所示； 

式中μ--样本集X的均值向量; 

∑--为协方差矩阵，

在计算马氏距离的过程中，若协方差矩阵为奇异矩阵，则协方差矩阵的逆阵不存在，将马氏距离无法直接求得，根据矩阵理论，任一秩为r的实对称半正定矩阵∑都可按公式(33)分解； 

∑=A^TG^-1A    (35) 

式中G--为r×r的对角阵，它由∑的非0特征值构成，G是非奇异的； 

A--为r×m矩阵，由G中特征值所对应的特征向量构成，它张成了∑的非退化子空间，且AA^T为r×r的单位阵； 

通过对公式(35)的分解，就可根据G的逆来求∑的伪逆矩阵，如公式(36)所示； 

∑^-1=A^TG^-1A   (36) 

用马氏距离来代替欧氏距离应用于模糊C均值聚类算法中； 

将公式(34)代入到传统模糊C均值聚类算法的全局代价函数(33)中，得到基于马氏距离的模糊C均值聚类算法的全局代价函数，如公式(37)所示，约束条件如公式(38)、(39)所示： 

依据GK算法，聚类中心、协方差矩阵和隶属度值的计算公式分别如公式(40)、(41)、(42)所示； 

式中D_ij--为马氏距离

8.根据权利要求7所述的基于多维信息的肺结节良恶性鉴别方法，其特征在于：所述的模糊C均值聚类算法的算法如下： 

(1)设置目标函数精度ε，模糊指数m，最大迭代次数T_m； 

(2)初始化模糊聚类中心； 

(3)由步骤(3)-(5)更新模糊划分矩阵U={μ_ij}和聚类中心Z={z₁,z₂,...,z_c}；； 

(4)若|J(t)-J(t-1)|<ε或t>T_m，则结束聚类；t←t+1并转到(3)步； 

(5)由U={μ_ij}得到各像素点分类结果。 

9.根据权利要求1所述的基于多维信息的肺结节良恶性鉴别方法，其特征在于：所述的步骤4、评价分类性能具体如下： 

采用的样本数据来自肺部图像数据库联盟，同时采用敏感性、准确性、特异性和假阳性四个指标来分别评价分类性能，具体为： 

(1)根据敏感性又称为真阳性率，判断对结节样本恶性的检测性能； 

(2)根据准确性，判断所有的检测结果中，检测恶性的样本占全部样本的比率； 

(3)根据特异性，判断对结节样本良性的检测性能； 

(4)根据假阳性率，判断检测结果中被错误检测为良性结节样本所占的比率； 

上式中：TP为真阳性结节的数量，FP为假阳性结节的数量，TN为假阴性结节的数量，FN为真阴性结节的数量。 

Images