CN106327480A - 一种甲状腺ct图像异常密度的检测方法 - Google Patents

Abstract

本发明公开了一种甲状腺CT图像异常密度的检测方法，包括获取CT图像、读取CT图像像素CT值、获得CT值异常的判断结果等步骤。本发明采用判断算法对甲状腺CT图像自动实施寻找异常密度所在的位置区域，自动判断得出甲状腺CT图像是否存在异常密度值的结果，实现了快速而及时地为医生提供异常密度信息并自动帮助医生确定异常密度值所在的位置区域；又通过阈值优化算法对甲状腺CT图像的CT值阈值进行最优值计算，进一步提高了识别异常密度区域的准确度，既有效地降低了医生的漏诊率，又大大减轻了医生的工作量；该检测方法由于不存在人为主观因素影响，避免了人为主观因素对病理检查与其他检查的干扰，大大提高了检测系统的运行效率和医生的工作效率。

Description

一种甲状腺CT图像异常密度的检测方法

技术领域

本发明涉及一种甲状腺CT图像异常密度的检测方法。

背景技术

目前，随着CT的应用普及和甲状腺CT检查的数量增加，日常放射科医生阅片图像数量大幅度增加。

在由CT扫描仪器扫描获得的甲状腺CT图像上，正常甲状腺密度均匀，CT值一般在90-120HU范围，但是在某些情况下甲状腺内部的密度会发生变化，如出现高密度或低密度。高密度物质一般是甲状腺内部异常组织发生了钙化，使得局部密度变大，CT值升高，往往大于120HU，达到200HU，甚至更高；低密度物质一般是甲状腺内部出现了液体样成分，如腺体组织坏死、产生囊肿等病变，CT值降低，CT值为0-40HU左右。

一直以来，对于依据甲状腺CT图像判断甲状腺内部是否产生了过高或过低密度物质以提示甲状腺发生了异常的检测方法，普遍采用医生直接通过肉眼读片的方法来判断甲状腺CT图像上是否存在甲状腺内部产生了过高或过低密度的物质，从而判断甲状腺是否发生了异常，这种依靠医生人为主观判断的检测方法，由于医生在知识背景、专业水平及能力、从业经验、视觉疲劳程度以及情绪影响等各种因素的综合作用下，容易得出错误的结论，从而引起漏诊或误诊，使得患者不能得到及时而准确的治疗。

发明内容

本发明解决的技术问题是，克服现有技术的缺陷，提供一种采用判断算法并通过阈值优化算法实现了更准确地提示医生而能有效降低医生工作量和减少漏诊率的甲状腺CT图像异常密度的检测方法。

为了解决上述技术问题，本发明通过下述技术方案得以解决：

一种甲状腺CT图像异常密度的检测方法，包括以下步骤：

步骤一：通过输入模块对由CT扫描仪器扫描获取的甲状腺CT原始图像进行轮廓分割，自动提取甲状腺CT横断面图像，再将甲状腺CT横断面图像经过滤波降噪处理，获得用于检测异常密度值的甲状腺CT图像，然后将甲状腺CT图像导入到检测系统中；

步骤二：通过纹理特征提取模块对步骤一中导入的甲状腺CT图像进行图像纹理分析，获取甲状腺CT图像的图像特征，依据所述的图像特征读取甲状腺CT图像各个像素的CT值；其中，所述的图像特征包含13维灰度共生纹理特征和15维灰度梯度共生纹理特征；

步骤三：先通过阈值优化算法计算并输出与甲状腺CT图像像素的CT值相对应的最优低密度阈值和最优高密度阈值，再将甲状腺CT图像中每2×2像素的CT值分别与最优低密度阈值和最优高密度阈值进行比较，然后通过判断算法对所述的每2×2像素所形成的点阵区域进行判断并获得判断结果，最后将判断结果输出并提示医生；其中，判断算法包括：若所述的每2×2像素中的每个像素的CT值均小于最优低密度阈值，则得出所述的每2×2像素所形成的点阵区域为低密度异常区域的判断结果，从而得出甲状腺CT图像在低密度异常区域存在CT值异常的判断结果；若所述的每2×2像素中的每个像素的CT值均大于最优高密度阈值，则得出所述的每2×2像素所形成的点阵区域为高密度异常区域的判断结果，从而得出甲状腺CT图像在高密度异常区域存在CT值异常的判断结果。本发明中，通过输出低密度异常区域和/或高密度异常区域被点亮的甲状腺CT图像，从而提示医生甲状腺CT图像中的被点亮的区域为甲状腺CT图像存在CT值异常的区域。

为了检测当前的甲状腺CT图像中相邻的每2×2像素同时大于最优高密度阈值或小于最优低密度阈值，以每2×2像素相邻的四个像素作为比较对象的目的是为了防止个别像素由于噪声影响准确度。

作为优选，采用MATLAB作为计算分析工具对步骤三中的阈值优化算法进行运算，具体步骤如下：

步骤A：初始化低密度阈值low为40HU，最高准确率的初始值ACC_best为0；

步骤B：初始化高密度阈值high为140HU；

步骤C：根据与当前的甲状腺CT图像像素的CT值相对应的当前的低密度阈值low和当前的高密度阈值high的值，计算出当前的准确率ACC、当前的敏感度SEN、当前的特异性SPC、当前的阳性预测率PPV、当前的阴性预测率NPV；

步骤D：比较当前的准确率ACC与最高准确率ACC_best的值：若当前的准确率ACC小于最高准确率ACC_best，则直接进入下一步骤E；若当前的准确率ACC大于最高准确率ACC_best，先记录数据，该数据包括low_best、high_best、ACC_best、SEN_best、SPC_best、PPV_best和NPV_best，再分别把当前的低密度阈值low、高密度阈值high、当前的准确率ACC、当前的敏感度SEN、当前的特异性SPC、当前的阳性预测率PPV、当前的阴性预测率NPV赋值给low_best、high_best、ACC_best、SEN_best、SPC_best、PPV_best和NPV_best，并进入下一步骤E；

步骤E：比较当前的高密度阈值high与160HU：若高密度阈值high小于160HU，则高密度阈值high自动加1，然后返回到步骤C；若高密度阈值high大于160HU，则进入下一步骤F；

步骤F：比较低密度阈值low与60HU：若低密度阈值low小于60HU，则低密度阈值low自动加1，然后返回到步骤C；若低密度阈值low大于60HU，则进入下一步骤G；

步骤G：输出最优高密度阈值high_best和最优低密度阈值low_best，优化计算运行结束。

作为优选，最优高密度阈值和最优低密度阈值分别为140-160HU和40-60HU。本发明中，最优高密度阈值和最优低密度阈值是基于试验基础上获得的经验值来优化计算，每次增加1HU，比较利用前后两个不同参数计算得到的判断准确率ACC，经过20乘20＝400次循环得到最高的ACC值，同时记录相应的参数和判断结果。

作为优选，最优高密度阈值和最优低密度阈值分别为145-155HU和45-55HU。

作为优选，最优高密度阈值和最优低密度阈值分别为145-150HU和45-50HU。

作为优选，最优高密度阈值为140HU或145HU或150HU或155HU或160HU，最优低密度阈值为40HU或45HU或50HU或55HU或60HU。本发明中，阈值的选择非常关键，直接影响到该检测方法的准确性。本发明采用多种阈值组合的方法进行比较。当设定高密度阈值时，应考虑甲状腺内钙化的程度和范围，有些钙化程度不深，钙化处的CT值不超过140HU，也有的钙化完全，钙化处的CT值超过200HU，如果高密度阈值取得过高，就难以检测出一般的钙化值，容易造成高密度的疏漏，所以高密度阈值分别取140HU、145HU、150HU、155HU、160HU；当设定低密度阈值时，也和组织坏死液化的程度有关，分别取40HU、45HU、50HU、55HU、60HU。

作为优选，步骤二中获取所述的13维灰度共生纹理特征过程为：

在甲状腺CT图像中任取一像素点A(x，y)，然后获得与像素点A距离为d的另一个像素点B(x+a，y+b)，像素点A和B形成一个点对，记录这一点对的灰度值(i，j)，固定a和b的值，使像素点A(x，y)在甲状腺CT图像上移动，获得种像素点对组合，N_g为图像灰度级别值，i∈[0,N_g-1]，j∈[0,N_g-1]，改变d和θ，θ为像素点对连线与水平先的夹角，统计甲状腺CT图像中的像素点对的灰度值，构成灰度共生矩阵P(i,j,θ,d)，

$P(i,j,\mathrm{\θ},d)=\#\{(x,y),(x+dsin\mathrm{\θ},y+dcos\mathrm{\θ})\left|\begin{array}{c}F(x,y)=i\\ F(x+dsin\mathrm{\θ},y+dcos\mathrm{\θ})=j\end{array}\right.\}$

其中#{x}是集合x中的所有元素的个数；然后根据灰度共生矩阵提取13维纹理特征，分别为角二阶矩、相关度、熵、对比度、逆差矩、平均和、和熵、方差和、方差、差平均、惯性、差方差、差熵。

作为优选，步骤二中所述的13维灰度共生纹理特征的公式为：

角二阶矩：

$f_1=\underset{i=1}{\overset{N_g}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}{\[P(i,j,\mathrm{\θ},d)\]}^2;$

角二阶矩是灰度共生矩阵中所有像素的平方和，表现了图像灰度的分布均匀程度和纹理的粗细程度。

相关度：

$f_2=\underset{i=1}{\overset{N_g}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}\[i\×j\×p(i,j,\mathrm{\θ},d)-{\mathrm{\μ}}_1\×{\mathrm{\μ}}_2\]/({\mathrm{\σ}}_1\×{\mathrm{\σ}}_2)$

其中，

${\mathrm{\μ}}_1=\underset{i=1}{\overset{N_g}{\Σ}}i\[\underset{j=1}{\overset{N_g}{\Σ}}P(i,j,\mathrm{\θ},d)\],{\mathrm{\μ}}_2=\underset{i=1}{\overset{N_g}{\Σ}}j\[\underset{j=1}{\overset{N_g}{\Σ}}P(i,j,\mathrm{\θ},d)\],$

${{\mathrm{\σ}}_1}^2=\underset{i=1}{\overset{N_g}{\Σ}}{(i-{\mathrm{\μ}}_1)}^2\[\underset{j=1}{\overset{N_g}{\Σ}}P(i,j,\mathrm{\θ},d)\],$

${{\mathrm{\σ}}_2}^2=\underset{j=1}{\overset{N_g}{\Σ}}{(j-{\mathrm{\μ}}_2)}^2\[\underset{i=1}{\overset{N_g}{\Σ}}P(i,j,\mathrm{\θ},d)\];$

相关度反映了灰度共生矩阵中的元素在水平方向和竖直方向的变化情况，当灰度共生矩阵中像素的元素分布均匀且数值接近，则相关度值较高。对于有方向性的纹理图像，如果特定θ方向的相关度较强，则可以反映了图像纹理的主要走向。

熵：

$f_3=-\underset{i=1}{\overset{N_g}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}P(i,j,\mathrm{\θ},d)\×{\log }_{2}P(i,j,\mathrm{\θ},d);$

熵是灰度共生矩阵中信息量的表征，反映了灰度共生矩阵中元素随机性的程度，可以反映图像中纹理的复杂程度。当图像中所有像素分布均匀，没有任何纹理时，灰度共生矩阵的元素几乎全为0，则熵值较小；当图像中纹理非常复杂，像素分布杂乱时，灰度共生矩阵中的元素会存在较大的差异，则熵值较大。

对比度：

$f_4=\underset{i=1}{\overset{N_g}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}{(i-j)}^2\×{\[P(i,j,\mathrm{\θ},d)\]}^2;$

对比度是灰度共生矩阵相对主对角线的惯性矩。若灰度共生矩阵中的元素距离主对角线越远且元素值越大，则对比度较大，反之亦然。而距离主对角线越远，则图像中纹理基元的灰度差值越大；此类元素越多，则反映图像中此类纹理基元越明显。对比度可以反映图像中纹理基元强度，即沟纹的深浅程度。纹理基元表现越强，沟纹越深，则对比度越大；反之亦然。

逆差矩：

$f_5=\underset{i=1}{\overset{N_g}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}\frac{P(i,j,\mathrm{\θ},d)}{1+{(i-j)}^2};$

逆差矩也可称为局部平稳，可以反映图像纹理的同质性。当图像的纹理在不同区域变化较小时，图像的像素点对灰度更为接近，则逆差矩较大。

平均和：

$f_6=\underset{k=2}{\overset{2N_g}{\Σ}}{\mathrm{kP}}_{x+y}\left(k\right),$

其中且|i+j|＝k，k＝2,3,…,2N_g；

平均和反映了图像中像素点对的平均灰度值，反映了图像纹理的明暗程度。

和熵：

$f_7=-\underset{k=2}{\overset{2N_g}{\Σ}}{P}_{x+y}\left(k\right)\[{\log }_2{P}_{x+y}\left(k\right)\],$

其中且|i+j|＝k，k＝2,3,…,2N_g；

和熵反映了图像中像素点对的分布变化，像素点对的变化越随机，分布范围越大，则和熵值越大。

方差和：

$f_8=\underset{k=2}{\overset{2N_g}{\Σ}}{(k-f_6)}^2{P}_{x+y}\left(k\right),$

其中且|i+j|＝k，k＝2,3,…,2N_g；

方差和是图像中像素点对的变化程度，反映了图像纹理的周期性。图像中纹理的周期越大，则方差和越大。

方差：

其中m是灰度共生矩阵P(i,j,θ,d)所有元素的均值；

灰度共生矩阵的方差反映了图像像素点对的变化，图像中纹理分布越复杂，则方差越大。

差平均：

$f_{10}=\underset{k=0}{\overset{N_g-1}{\Σ}}{\mathrm{kP}}_{x-y}\left(k\right),$

其中且|i-j|＝k，k＝0,1,…,N_g-1；

差平均反映了图像的像素点对的平均灰度差异，如果图像纹理基元具有较强的灰度对比，并且分布空间较广，则差平均值较大。

惯性：

$f_{11}=\underset{i=1}{\overset{N_g}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}{(i-j)}^{2}P(i,j,\mathrm{\θ},d);$

惯性描绘了图像纹理基元的对比度。如果图像中的像素点对具有较大的灰度差异，像素点对的数量较多，则惯性值较大。

差方差：

$f_{12}=\underset{k=0}{\overset{N_g-1}{\Σ}}{\[k-\underset{k=0}{\overset{N_g-1}{\Σ}}{\mathrm{kP}}_{x-y}\left(k\right)\]}^2\×P_{x-y}\left(k\right),$

其中且|i-j|＝k，k＝0,1,…,Ng-1；

差方差是图像中像素点对的灰度差的方差，描绘了图像纹理的对比度。当图像的纹理越显著，对比度越强强烈，则差方差越大。

熵差：

$f_{13}=-\underset{k=0}{\overset{N_g}{\Σ}}{P}_{x-y}\left(k\right)\[{\log }_2\left(P_{x-y}\right(k\left)\right)\],$

其中且|i-j|＝k，k＝0,1,…,N_g-1。

作为优选，步骤二中获取所述的15维灰度梯度共生纹理特征过程为：对图像f求取梯度图像g，其中图像f大小为N_x×N_y，图像f的最高灰度级为N，灰度范围为[0，N-1]，则有其中g_x和g_y分别为边缘检测算子，对图像f在水平方向和竖直方向上计算得到梯度图像g，对灰度图像f和梯度图像g分别进行归一化处理，对应得到图像F和图像G，并且使得归一化后的图像G灰度分布范围为[0,N_g-1]，图像F的灰度分布范围为[0,N_f-1]，计算公式如下：

F＝[f×(N_f-1)/f_max]+1，G＝[g×(N_g-1)/g_max]+1，

其中f_max是灰度图像f中的灰度最大值，g_max是梯度图像g的灰度最大值；

对图像f中任一像素点(x，y)，归一化处理后的灰度值F(x，y)为i，且归一化后的梯度值G(x，y)为j，得到灰度梯度共生矩阵中坐标为(i，j)的元素，记为h(i，j)，

h(i,j)＝#{(x,y)|F(x,y)＝i,G(x,y)＝j}

其中#{X}是集合X中的所有元素的个数；对h(i，j)进行归一化处理，得到

H(i,j)＝h(i,j)/(N_g×N_f)

其中，i＝1,2,…,N_f，j＝1,2,…,N_g，根据灰度梯度共生矩阵提取15维纹理特征，分别为小梯度优势、大梯度优势、灰度分布的不均匀性、梯度分布的不均匀性、能量、灰度平均值、梯度平均值、灰度标准差、梯度标准差、相关性、灰度熵、梯度熵、混合熵、差分矩、逆差矩。对于灰度图像，不仅可以使用灰度信息统计分析图像的纹理特征，而且可以使用灰度的邻域变化信息分析图像纹理。梯度就是灰度的邻域变化信息，构成了图像的边缘轮廓。灰度梯度共生矩阵可以综合图像的灰度和梯度信息，更为全面的反映图像中每个像素点的灰度和梯度的分布特征，反映像素点与其邻域各像素点的空间位置关系。

作为优选，步骤二中所述的15维灰度梯度共生纹理特征的公式为：

小梯度优势：

其中

小梯度优势反映了图像中小梯度的分布强度。当图像的灰度变化较平缓，图像大部分区域的梯度值较小，则小梯度优势值大。

大梯度优势：

其中

大梯度优势反映了图像中的大梯度的分布强度。如果图像中的灰度变化非常剧烈，大梯度在图像中分布较广泛，则大梯度优势值大。

灰度分布的不均匀性：

其中

如果图像中不同区域的灰度差别比较大，则灰度分布的不均匀性值较大。

梯度分布的不均匀性：

其中

如果图像中不同区域的梯度变化较大，则梯度分布的不均匀性值较大。

能量：

$t_5=\underset{i=1}{\overset{N_f}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}{\[H(i,j)\]}^2;$

能量反映了图像纹理的强度和密集程度。如果图像中灰度和梯度变化比较强烈，且较为密集，则能量值较大。

灰度平均值：

$t_6=\underset{i=1}{\overset{N_f}{\Σ}}i\[\underset{j=1}{\overset{N_g}{\Σ}}H(i,j)\];$

灰度平均值反映了图像中的灰度强度，即图像的整体亮度。

梯度平均值：

$t_7=\underset{j=1}{\overset{N_g}{\Σ}}j\[\underset{i=1}{\overset{N_f}{\Σ}}H(i,j)\];$

梯度平均值反映了图像中的梯度的平均强度。

灰度标准差：

$t_8={\{\underset{i=1}{\overset{N_f}{\Σ}}{(i-t_6)}^2\[\underset{j=1}{\overset{N_g}{\Σ}}H(i,j)\]\}}^{1/2};$

灰度标准差反映了图像中灰度的变化强度

梯度标准差：

$t_9={\{\underset{j=1}{\overset{N_g}{\Σ}}{(j-t_7)}^2\[\underset{i=1}{\overset{N_f}{\Σ}}H(i,j)\]\}}^{1/2};$

梯度标准差反映了图像中梯度的变化强度。如果图像中的梯度大小较为接近，则梯度标准差值较小，反之亦然。

相关性：

$t_{10}=\frac{1}{t_8{t}_9}\underset{i=1}{\overset{N_f}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}(i-t_6)(j-t_7)H(i,j);$

相关性反映了图像中灰度和梯度的相互关系。如果图像中纹理基元的分布比较规则，那么灰度和梯度的变化存在一致性，则相关性值较大。

灰度熵：

$t_{11}=-\underset{i=1}{\overset{N_f}{\Σ}}\{\[\underset{j=1}{\overset{N_g}{\Σ}}H(i,j)\]{\log }_2\[\underset{j=1}{\overset{N_g}{\Σ}}H(i,j)\]\};$

灰度熵反映了图像中灰度分布的规律性，即图像中的平均信息量。

梯度：

$t_{12}=-\underset{j=1}{\overset{N_g}{\Σ}}\{\[\underset{i=1}{\overset{N_f}{\Σ}}H(i,j)\]{\log }_2\[\underset{i=1}{\overset{N_f}{\Σ}}H(i,j)\]\};$

梯度熵反映了图像中梯度分布的规律性，即梯度图像的平均信息量。

混合熵：

$t_{13}=-\underset{i=1}{\overset{N_f}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}H(i,j)\[{\log }_{2}H(i,j)\];$

混合熵反映了灰度-梯度共生矩阵的平均信息量。如果图像中灰度和梯度的变化较复杂，则混合熵值较大。

差分矩：

$t_{14}=\underset{i=1}{\overset{N_f}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}{(i-j)}^{2}H(i,j);$

差分矩描绘了图像纹理基元的对比度。如果图像中像素点对的灰度较高并且灰度变化很小，那么图像的差分矩较大。

逆差矩：

$t_{15}=\underset{i=1}{\overset{N_f}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}\frac{1}{1+{(i-j)}^2}H(i,j);$

逆差矩，即局部平稳，可以反映图像纹理的同质性。当图像纹理基元的灰度与梯度非常接近，并且覆盖图像的大部分区域，则逆差矩较大。

本发明中，步骤一中对甲状腺CT图像进行轮廓分割以及自动提取的过程为：先由CT扫描仪器对甲状腺肿块患者进行扫描，获取得到甲状腺CT原始图像，再在甲状腺CT原始图像中选取肿块最大横径的横断面图像，在所述的最大横径的横断面图像中将含有肿块的区域边缘手动勾画，然后再将勾画的区域所在的图像提取出来，并设定图像灰度级为256级，将所有提取的不同的区域图像构成病变甲状腺图像集。

本发明由于采用了以上技术方案，具有显著的技术效果：采用判断算法对甲状腺CT图像自动实施寻找异常密度所在的位置区域，自动判断得出甲状腺CT图像是否存在异常密度值的结果，实现了快速而及时地为医生提供异常密度信息并自动帮助医生确定异常密度值所在的位置区域；又通过阈值优化算法对甲状腺CT图像的CT值阈值进行最优值计算，进一步提高了识别异常密度区域的准确度，既有效地降低了医生的漏诊率，又大大减轻了医生的工作量；该检测方法由于不存在人为主观因素影响，避免了人为主观因素对病理检查与其他检查的干扰，大大提高了检测系统的运行效率和医生的工作效率。

附图说明

图1为本发明甲状腺CT图像异常密度的检测方法实施例的系统功能流程图。

图2为本发明阈值优化算法实施例的程序流程图。

图3为本发明当最优高密度阈值为150HU时当前的低密度阈值与准确率之间的关系图。

图4为本发明当最优低密度阈值为45HU时当前的高密度阈值与准确率之间的关系图。

图5为本发明获取13维灰度共生纹理特征过程中灰度矩阵像素点对实施例的结构示意图。

具体实施方式

下面结合附图与实施例对本发明作进一步详细描述。

一种甲状腺CT图像异常密度的检测方法，如图1-5所示，包括以下步骤：

步骤一：通过输入模块对由CT扫描仪器扫描获取的甲状腺CT原始图像进行轮廓分割，自动提取甲状腺CT横断面图像，再将甲状腺CT横断面图像经过滤波降噪处理，获得用于检测异常密度值的甲状腺CT图像，然后将甲状腺CT图像导入到检测系统中；

通过CT扫描仪对患者扫描采集甲状腺CT原始图像，CT扫描仪采用西门子Sensation 16层螺旋CT，采集甲状腺肿块患者的CT平扫横断面图像，图像格式为DICOM。CT设备扫描参数为球管电压120kV，管电流220mAs，层厚为2-3mm，层间距2-3mm，螺距为1-1.5，图像重建类型为B40，软组织显示窗，横断面图像的分辨率为512×512像素，每位患者10-15个横断面图像。然后在CT平扫横断面图像中，选取肿块所在最大横径的横断面图像，进行轮廓分隔，自动提取甲状腺CT横断面图像。

先在甲状腺CT原始图像中选取肿块最大横径的横断面图像，再使用microMRI软件在所述的最大横径的横断面图像中将含有肿块的区域边缘手动勾画，然后再将勾画的区域所在的图像提取出来，并设定图像灰度级为256级，将所有提取的不同的区域图像构成病变甲状腺图像集。在获得病变甲状腺图像集的同时，通过检测系统中存储的已得到证实的病理结果将病变甲状腺CT图像相应的标记为恶性或良性，例如良性肿块包括甲状腺腺瘤、甲状腺肿和桥本氏病灯良性病变，恶性肿块包括乳头状腺癌和滤泡状腺癌。

步骤二：通过纹理特征提取模块对步骤一中导入的甲状腺CT图像进行图像纹理分析，获取甲状腺CT图像的图像特征，依据所述的图像特征读取甲状腺CT图像各个像素的CT值；其中，所述的图像特征包含13维灰度共生纹理特征和15维灰度梯度共生纹理特征；

步骤三：先通过阈值优化算法计算并输出与甲状腺CT图像像素的CT值相对应的最优低密度阈值和最优高密度阈值，再将甲状腺CT图像中每2×2像素的CT值分别与最优低密度阈值和最优高密度阈值进行比较，然后通过判断算法对所述的每2×2像素所形成的点阵区域进行判断并获得判断结果，最后将判断结果输出并提示医生；其中，判断算法包括：若所述的每2×2像素中的每个像素的CT值均小于最优低密度阈值，则得出所述的每2×2像素所形成的点阵区域为低密度异常区域的判断结果，从而得出甲状腺CT图像在低密度异常区域存在CT值异常的判断结果；若所述的每2×2像素中的每个像素的CT值均大于最优高密度阈值，则得出所述的每2×2像素所形成的点阵区域为高密度异常区域的判断结果，从而得出甲状腺CT图像在高密度异常区域存在CT值异常的判断结果。

本实施例中，采用MATLAB作为计算分析工具对步骤三中的阈值优化算法进行运算，具体步骤如下：

步骤A：初始化低密度阈值low为40HU，最高准确率的初始值ACC_best为0；

步骤B：初始化高密度阈值high为140HU；

步骤C：根据与当前的甲状腺CT图像像素的CT值相对应的当前的低密度阈值low和当前的高密度阈值high的值，计算出当前的准确率ACC、当前的敏感度SEN、当前的特异性SPC、当前的阳性预测率PPV、当前的阴性预测率NPV；

步骤D：比较当前的准确率ACC与最高准确率ACC_best的值：若当前的准确率ACC小于最高准确率ACC_best，则直接进入下一步骤E；若当前的准确率ACC大于最高准确率ACC_best，先记录数据，该数据包括low_best、high_best、ACC_best、SEN_best、SPC_best、PPV_best和NPV_best，再分别把当前的低密度阈值low、高密度阈值high、当前的准确率ACC、当前的敏感度SEN、当前的特异性SPC、当前的阳性预测率PPV、当前的阴性预测率NPV赋值给low_best、high_best、ACC_best、SEN_best、SPC_best、PPV_best和NPV_best，并进入下一步骤E；

步骤E：比较当前的高密度阈值high与160HU：若高密度阈值high小于160HU，则高密度阈值high自动加1，然后返回到步骤C；若高密度阈值high大于160HU，则进入下一步骤F；

步骤F：比较低密度阈值low与60HU：若低密度阈值low小于60HU，则低密度阈值low自动加1，然后返回到步骤C；若低密度阈值low大于60HU，则进入下一步骤G；

步骤G：输出最优高密度阈值high_best和最优低密度阈值low_best，优化计算运行结束。

本实施例中，最优高密度阈值为140HU或145HU或150HU或155HU或160HU，最优低密度阈值为40HU或45HU或50HU或55HU或60HU。

本实施例中，依据阈值优化算法对高密度阈值和低密度阈值进行优化计算，若优化计算之后的最优高密度阈值和最优低密度阈值分别为151HU和57HU，则ACC、SEN、SPC、PPV、NPV对应的计算结果分别为0.8511、0.8060、0.8984、0.8926、0.8156，优于经验取值的最高准确率。

本实施例中，步骤二中获取所述的13维灰度共生纹理特征过程为：

如图5所示，在甲状腺CT图像中任取一像素点A(x，y)，然后获得与像素点A距离为d的另一个像素点B(x+a，y+b)，像素点A和B形成一个点对，记录这一点对的灰度值(i，j)，固定a和b的值，使像素点A(x，y)在甲状腺CT图像上移动，获得种像素点对组合，N_g为图像灰度级别值，i∈[0,N_g-1]，j∈[0,N_g-1]，改变d和θ，θ为像素点对连线与水平先的夹角，统计甲状腺CT图像中的像素点对的灰度值，构成灰度共生矩阵P(i,j,θ,d)，

$P(i,j,\mathrm{\θ},d)=\#\{(x,y),(x+dsin\mathrm{\θ},y+dcos\mathrm{\θ})\left|\begin{array}{c}F(x,y)=i\\ F(x+dsin\mathrm{\θ},y+dcos\mathrm{\θ})=j\end{array}\right.\}$

其中#{x}是集合x中的所有元素的个数；预先设定d值，分别计算0°、45°、90°和135°四个方向的灰度共生矩阵，如果灰度范围为[0,N_g-1]，则每个方向的灰度共生矩阵的尺寸是N_g×N_g。本实施例中，预先设定d＝1，分别计算0°、45°、90°和135°的灰度共生矩阵，然后基于每个方向的矩阵计算纹理特征，将四个方向的纹理特征求取平均值，得到旋转不变的纹理特征。然后根据灰度共生矩阵提取13维纹理特征，分别为角二阶矩、相关度、熵、对比度、逆差矩、平均和、和熵、方差和、方差、差平均、惯性、差方差、差熵。

本实施例中，步骤二中所述的13维灰度共生纹理特征的公式为：

角二阶矩：

$f_1=\underset{i=1}{\overset{N_g}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}{\[P(i,j,\mathrm{\θ},d)\]}^2;$

相关度：

$f_2=\underset{i=1}{\overset{N_g}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}\[i\×j\×p(i,j,\mathrm{\θ},d)-{\mathrm{\μ}}_1\×{\mathrm{\μ}}_2\]/({\mathrm{\σ}}_1\×{\mathrm{\σ}}_2)$

其中，

${\mathrm{\μ}}_1=\underset{i=1}{\overset{N_g}{\Σ}}i\[\underset{j=1}{\overset{N_g}{\Σ}}P(i,j,\mathrm{\θ},d)\],{\mathrm{\μ}}_2=\underset{i=1}{\overset{N_g}{\Σ}}j\[\underset{j=1}{\overset{N_g}{\Σ}}P(i,j,\mathrm{\θ},d)\],$

${{\mathrm{\σ}}_1}^2=\underset{i=1}{\overset{N_g}{\Σ}}{(i-{\mathrm{\μ}}_1)}^2\[\underset{j=1}{\overset{N_g}{\Σ}}P(i,j,\mathrm{\θ},d)\],$

${{\mathrm{\σ}}_2}^2=\underset{j=1}{\overset{N_g}{\Σ}}{(j-{\mathrm{\μ}}_2)}^2\[\underset{i=1}{\overset{N_g}{\Σ}}P(i,j,\mathrm{\θ},d)\];$

熵：

$f_3=-\underset{i=1}{\overset{N_g}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}P(i,j,\mathrm{\θ},d)\×{\log }_{2}P(i,j,\mathrm{\θ},d);$

对比度：

$f_4=\underset{i=1}{\overset{N_g}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}{(i-j)}^2\×{\[P(i,j,\mathrm{\θ},d)\]}^2;$

逆差矩：

$f_5=\underset{i=1}{\overset{N_g}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}\frac{P(i,j,\mathrm{\θ},d)}{1+{(i-j)}^2};$

平均和：

$f_6=\underset{k=2}{\overset{2N_g}{\Σ}}{\mathrm{kP}}_{x+y}\left(k\right),$

其中且|i+j|＝k，k＝2,3,…,2N_g；和熵：

$f_7=-\underset{k=2}{\overset{2N_g}{\Σ}}{P}_{x+y}\left(k\right)\[{\log }_2{P}_{x+y}\left(k\right)\],$

其中且|i+j|＝k，k＝2,3,…,2N_g；方差和：

$f_8=\underset{k=2}{\overset{2N_g}{\Σ}}{(k-f_6)}^2{P}_{x+y}\left(k\right),$

其中且|i+j|＝k，k＝2,3,…,2N_g；

方差：

其中m是灰度共生矩阵P(i,j,θ,d)所有元素的均值；

差平均：

$f_{10}=\underset{k=0}{\overset{N_g-1}{\Σ}}{\mathrm{kP}}_{x-y}\left(k\right),$

其中且|i-j|＝k，k＝0,1,…,N_g-1；

惯性：

$f_{11}=\underset{i=1}{\overset{N_g}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}{(i-j)}^{2}P(i,j,\mathrm{\θ},d);$

差方差：

$f_{12}=\underset{k=0}{\overset{N_g-1}{\Σ}}{\[k-\underset{k=0}{\overset{N_g-1}{\Σ}}{\mathrm{kP}}_{x-y}\left(k\right)\]}^2\×P_{x-y}\left(k\right),$

其中且|i-j|＝k，k＝0,1,…,N_g-1；

熵差：

$f_{13}=-\underset{k=0}{\overset{N_g}{\Σ}}{P}_{x-y}\left(k\right)\[{\log }_2\left(P_{x-y}\right(k\left)\right)\],$

其中且|i-j|＝k，k＝0,1,…,N_g-1。

本实施例中，步骤二中获取所述的15维灰度梯度共生纹理特征过程为：对图像f求取梯度图像g，其中图像f大小为N_x×N_y，图像f的最高灰度级为N，灰度范围为[0，N-1]，则有其中g_x和g_y分别为边缘检测算子，本实施例中采用尺寸为3*3的Sobel算子求取梯度图像，Sobel算子是水平方向和竖直方向的两个边界检测算子，其分别为：

对图像f在水平方向和竖直方向上计算得到梯度图像g，对灰度图像f和梯度图像g分别进行归一化处理，对应得到图像F和图像G，并且使得归一化后的图像G灰度分布范围为[0,N_g-1]，图像F的灰度分布范围为[0,N_f-1]，计算公式如下：

F＝[f×(N_f-1)/f_max]+1，G＝[g×(N_g-1)/g_max]+1，

其中f_max是灰度图像f中的灰度最大值，g_max是梯度图像g的灰度最大值；

对图像f中任一像素点(x，y)，归一化处理后的灰度值F(x，y)为i，且归一化后的梯度值G(x，y)为j，得到灰度梯度共生矩阵中坐标为(i，j)的元素，记为h(i，j)，

h(i,j)＝#{(x,y)|F(x,y)＝i,G(x,y)＝j}

其中#{X}是集合X中的所有元素的个数；对h(i，j)进行归一化处理，得到

H(i,j)＝h(i,j)/(N_g×N_f)

其中，i＝1,2,…,N_f，j＝1,2,…,N_g，根据灰度梯度共生矩阵提取15维纹理特征，分别为小梯度优势、大梯度优势、灰度分布的不均匀性、梯度分布的不均匀性、能量、灰度平均值、梯度平均值、灰度标准差、梯度标准差、相关性、灰度熵、梯度熵、混合熵、差分矩、逆差矩。

本实施例中，步骤二中所述的15维灰度梯度共生纹理特征的公式为：

小梯度优势：

其中

大梯度优势：

其中

灰度分布的不均匀性：

其中

梯度分布的不均匀性：

其中

能量：

$t_5=\underset{i=1}{\overset{N_f}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}{\[H(i,j)\]}^2;$

灰度平均值：

$t_6=\underset{i=1}{\overset{N_f}{\Σ}}i\[\underset{j=1}{\overset{N_g}{\Σ}}H(i,j)\];$

梯度平均值：

$t_7=\underset{j=1}{\overset{N_g}{\Σ}}j\[\underset{i=1}{\overset{N_f}{\Σ}}H(i,j)\];$

灰度标准差：

$t_8={\{\underset{i=1}{\overset{N_f}{\Σ}}{(i-t_6)}^2\[\underset{j=1}{\overset{N_g}{\Σ}}H(i,j)\]\}}^{1/2};$

梯度标准差：

$t_9={\{\underset{j=1}{\overset{N_g}{\Σ}}{(j-t_7)}^2\[\underset{i=1}{\overset{N_f}{\Σ}}H(i,j)\]\}}^{1/2};$

相关性：

$t_{10}=\frac{1}{t_8{t}_9}\underset{i=1}{\overset{N_f}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}(i-t_6)(j-t_7)H(i,j);$

灰度熵：

$t_{11}=-\underset{i=1}{\overset{N_f}{\Σ}}\{\[\underset{j=1}{\overset{N_g}{\Σ}}H(i,j)\]{\log }_2\[\underset{j=1}{\overset{N_g}{\Σ}}H(i,j)\]\};$

梯度：

$t_{12}=-\underset{j=1}{\overset{N_g}{\Σ}}\{\[\underset{i=1}{\overset{N_f}{\Σ}}H(i,j)\]{\log }_2\[\underset{i=1}{\overset{N_f}{\Σ}}H(i,j)\]\};$

混合熵：

$t_{13}=-\underset{i=1}{\overset{N_f}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}H(i,j)\[{\log }_{2}H(i,j)\];$

差分矩：

$t_{14}=\underset{i=1}{\overset{N_f}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}{(i-j)}^{2}H(i,j);$

逆差矩：

$t_{15}=\underset{i=1}{\overset{N_f}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}\frac{1}{1+{(i-j)}^2}H(i,j).$

通过以上的纹理特征提取，得到28维特征参数，参数总结如下表所示：

总之，以上所述仅为本发明的较佳实施例，凡依本发明申请专利范围所作的均等变化与修饰，皆应属本发明专利的涵盖范围。

Claims (10)

1.一种甲状腺CT图像异常密度的检测方法，其特征在于：包括以下步骤：

步骤一：通过输入模块对由CT扫描仪器扫描获取的甲状腺CT原始图像进行轮廓分割，自动提取甲状腺CT横断面图像，再将甲状腺CT横断面图像经过滤波降噪处理，获得用于检测异常密度值的甲状腺CT图像，然后将甲状腺CT图像导入到检测系统中；

步骤二：通过纹理特征提取模块对步骤一中导入的甲状腺CT图像进行图像纹理分析，获取甲状腺CT图像的图像特征，依据所述的图像特征读取甲状腺CT图像各个像素的CT值；其中，所述的图像特征包含13维灰度共生纹理特征和15维灰度梯度共生纹理特征；

步骤三：先通过阈值优化算法计算并输出与甲状腺CT图像像素的CT值相对应的最优低密度阈值和最优高密度阈值，再将甲状腺CT图像中每2×2像素的CT值分别与最优低密度阈值和最优高密度阈值进行比较，然后通过判断算法对所述的每2×2像素所形成的点阵区域进行判断并获得判断结果，最后将判断结果输出并提示医生；其中，判断算法包括：若所述的每2×2像素中的每个像素的CT值均小于最优低密度阈值，则得出所述的每2×2像素所形成的点阵区域为低密度异常区域的判断结果，从而得出甲状腺CT图像在低密度异常区域存在CT值异常的判断结果；若所述的每2×2像素中的每个像素的CT值均大于最优高密度阈值，则得出所述的每2×2像素所形成的点阵区域为高密度异常区域的判断结果，从而得出甲状腺CT图像在高密度异常区域存在CT值异常的判断结果。

2.根据权利要求1所述的一种甲状腺CT图像异常密度的检测方法，其特征在于：采用MATLAB作为计算分析工具对步骤三中的阈值优化算法进行运算，具体步骤如下：

步骤A：初始化低密度阈值low为40HU，最高准确率的初始值ACC_best为0；

步骤B：初始化高密度阈值high为140HU；

步骤C：根据与当前的甲状腺CT图像像素的CT值相对应的当前的低密度阈值low和当前的高密度阈值high的值，计算出当前的准确率ACC、当前的敏感度SEN、当前的特异性SPC、当前的阳性预测率PPV、当前的阴性预测率NPV；

步骤D：比较当前的准确率ACC与最高准确率ACC_best的值：若当前的准确率ACC小于最高准确率ACC_best，则直接进入下一步骤E；若当前的准确率ACC大于最高准确率ACC_best，先记录数据，该数据包括low_best、high_best、ACC_best、SEN_best、SPC_best、PPV_best和NPV_best，再分别把当前的低密度阈值low、高密度阈值high、当前的准确率ACC、当前的敏感度SEN、当前的特异性SPC、当前的阳性预测率PPV、当前的阴性预测率NPV赋值给low_best、high_best、ACC_best、SEN_best、SPC_best、PPV_best和NPV_best，并进入下一步骤E；

步骤E：比较当前的高密度阈值high与160HU：若高密度阈值high小于160HU，则高密度阈值high自动加1，然后返回到步骤C；若高密度阈值high大于160HU，则进入下一步骤F；

步骤F：比较低密度阈值low与60HU：若低密度阈值low小于60HU，则低密度阈值low自动加1，然后返回到步骤C；若低密度阈值low大于60HU，则进入下一步骤G；

步骤G：输出最优高密度阈值high_best和最优低密度阈值low_best，优化计算运行结束。

3.根据权利要求1或2所述的一种甲状腺CT图像异常密度的检测方法，其特征在于：最优高密度阈值和最优低密度阈值分别为140-160HU和40-60HU。

4.根据权利要求3所述的一种甲状腺CT图像异常密度的检测方法，其特征在于：最优高密度阈值和最优低密度阈值分别为145-155HU和45-55HU。

5.根据权利要求1或2或4所述的一种甲状腺CT图像异常密度的检测方法，其特征在于：最优高密度阈值和最优低密度阈值分别为145-150HU和45-50HU。

6.根据权利要求5所述的一种甲状腺CT图像异常密度的检测方法，其特征在于：最优高密度阈值为140HU或145HU或150HU或155HU或160HU，最优低密度阈值为40HU或45HU或50HU或55HU或60HU。

7.根据权利要求1或2或4或6所述的一种甲状腺CT图像异常密度的检测方法，其特征在于：步骤二中获取所述的13维灰度共生纹理特征过程为：

在甲状腺CT图像中任取一像素点A(x，y)，然后获得与像素点A距离为d的另一个像素点B(x+a，y+b)，像素点A和B形成一个点对，记录这一点对的灰度值(i，j)，固定a和b的值，使像素点A(x，y)在甲状腺CT图像上移动，获得种像素点对组合，N_g为图像灰度级别值，i∈[0，N_g-1]，j∈[0，N_g-1]，改变d和θ，θ为像素点对连线与水平先的夹角，统计甲状腺CT图像中的像素点对的灰度值，构成灰度共生矩阵P(i，j，θ，d)，

$P(i,j,\mathrm{\θ},d)=\#\{(x,y),(x+d\sin \mathrm{\θ},y+d\cos \mathrm{\θ})|\begin{array}{c}F(x,y)=i\\ F(x+dsin\mathrm{\θ},y+dcos\mathrm{\θ})=j\end{array}\}$

其中#{x｝是集合x中的所有元素的个数；然后根据灰度共生矩阵提取13维纹理特征，分别为角二阶矩、相关度、熵、对比度、逆差矩、平均和、和熵、方差和、方差、差平均、惯性、差方差、差熵。

8.根据权利要求7所述的一种甲状腺CT图像异常密度的检测方法，其特征在于：步骤二中所述的13维灰度共生纹理特征的公式为：

角二阶矩：

$f_1=\underset{i=1}{\overset{N_g}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}{\[P(i,j,\mathrm{\θ},d)\]}^2;$

相关度：

$f_2=\underset{i=1}{\overset{N_g}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}\[i\×j\×P(i,j,\mathrm{\θ},d)-{\mathrm{\μ}}_1\×{\mathrm{\μ}}_2\]/({\mathrm{\σ}}_1\×{\mathrm{\σ}}_2)$

其中，

${\mathrm{\μ}}_1=\underset{i=1}{\overset{N_g}{\Σ}}i\[\underset{j=1}{\overset{N_g}{\Σ}}P(i,j,\mathrm{\θ},d)\],{\mathrm{\μ}}_2=\underset{i=1}{\overset{N_g}{\Σ}}j\[\underset{j=1}{\overset{N_g}{\Σ}}P(i,j,\mathrm{\θ},d)\],$

${{\mathrm{\σ}}_1}^2=\underset{i=1}{\overset{N_g}{\Σ}}{(i-{\mathrm{\μ}}_1)}^2\[\underset{j=1}{\overset{N_g}{\Σ}}P(i,j,\mathrm{\θ},d)\],$

${{\mathrm{\σ}}_2}^2=\underset{j=1}{\overset{N_g}{\Σ}}{(j-{\mathrm{\μ}}_2)}^2\[\underset{i=1}{\overset{N_g}{\Σ}}P(i,j,\mathrm{\θ},d)\];$

熵：

$f_3=-\underset{i=1}{\overset{N_g}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}P(i,j,\mathrm{\θ},d)\×{\log }_{2}P(i,j,\mathrm{\θ},d);$

对比度：

$f_4=\underset{i=1}{\overset{N_g}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}{(i-j)}^2\×{\[P(i,j,\mathrm{\θ},d)\]}^2;$

逆差矩：

$f_5=\underset{i=1}{\overset{N_g}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}\frac{P(i,j,\mathrm{\θ},d)}{1+{(i-j)}^2};$

平均和：

$f_6=\underset{k=2}{\overset{2N_g}{\Σ}}{\mathrm{kP}}_{x+y}\left(k\right),$

其中且|i+j|＝k，k＝2，3，…，2N_g；

和熵：

$f_7=-\underset{k=2}{\overset{2N_g}{\Σ}}{P}_{x+y}\left(k\right)\[{\log }_2{P}_{x+y}\left(k\right)\],$

其中且|i+j|＝k，k＝2，3，…，2N_g；

方差和：

$f_8=\underset{k=2}{\overset{2N_g}{\Σ}}{(k-f_6)}^2{P}_{x+y}\left(k\right),$

其中且|i+j|＝k，k＝2，3，…，2N_g；

方差：

其中m是灰度共生矩阵P(i，j，θ，d)所有元素的均值；

差平均：

$f_{10}=\underset{k=0}{\overset{N_g-1}{\Σ}}{\mathrm{kP}}_{x-y}\left(k\right),$

其中且|i-j|＝k，k＝0，1，…，N_g-1；

惯性：

$f_{11}=\underset{i=1}{\overset{N_g}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}{(i-j)}^{2}P(i,j,\mathrm{\θ},d);$

差方差：

$f_{12}=\underset{k=0}{\overset{N_g-1}{\Σ}}{\[k-\underset{k=0}{\overset{N_g-1}{\Σ}}{\mathrm{kP}}_{x-y}\left(k\right)\]}^2\×P_{x-y}\left(k\right),$

其中且|i-j|＝k，k＝0，1，…，N_g-1；

熵差：

$f_{13}=-\underset{k=0}{\overset{N_g}{\Σ}}{P}_{x-y}\left(k\right)\[{\log }_2\left(P_{x-y}\right(k\left)\right)\],$

其中且|i-j|＝k，k＝0，1，…，N_g-1。

9.根据权利要求1或2或4或6或8所述的一种甲状腺CT图像异常密度的检测方法，其特征在于：步骤二中获取所述的15维灰度梯度共生纹理特征过程为：对图像f求取梯度图像g，其中图像f大小为N_x×N_y，图像f的最高灰度级为N，灰度范围为[0，N-1]，则有其中g_x和g_y分别为边缘检测算子，对图像f在水平方向和竖直方向上计算得到梯度图像g，对灰度图像f和梯度图像g分别进行归一化处理，对应得到图像F和图像G，并且使得归一化后的图像G灰度分布范围为[0，N_g-1]，图像F的灰度分布范围为[0，N_f-1]，计算公式如下：

F＝[f×(N_f-1)/f_max]+1，G＝[g×(N_g-1)/g_max]+1，

其中f_max是灰度图像f中的灰度最大值，g_max是梯度图像g的灰度最大值；

对图像f中任一像素点(x，y)，归一化处理后的灰度值F(x，y)为i，且归一化后的梯度值G(x，y)为j，得到灰度梯度共生矩阵中坐标为(i，j)的元素，记为h(i，j)，

h(i，j)＝#{(x，y)|F(x，y)＝i，G(x，y)＝j}

其中#{X}是集合X中的所有元素的个数；对h(i，j)进行归一化处理，得到

H(i，j)＝h(i，j)/(N_g×N_f)

其中，i＝1，2，…，N_f，j＝1，2，…，N_g，根据灰度梯度共生矩阵提取15维纹理特征，分别为小梯度优势、大梯度优势、灰度分布的不均匀性、梯度分布的不均匀性、能量、灰度平均值、梯度平均值、灰度标准差、梯度标准差、相关性、灰度熵、梯度熵、混合熵、差分矩、逆差矩。

10.根据权利要求9所述的一种甲状腺CT图像异常密度的检测方法，其特征在于：步骤二中所述的15维灰度梯度共生纹理特征的公式为：

小梯度优势：

其中

大梯度优势：

其中

灰度分布的不均匀性：

其中

梯度分布的不均匀性：

其中

能量：

$t_5=\underset{i=1}{\overset{N_f}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}{\[H(i,j)\]}^2;$

灰度平均值：

$t_6=\underset{i=1}{\overset{N_f}{\Σ}}i\[\underset{j=1}{\overset{N_g}{\Σ}}H(i,j)\];$

梯度平均值：

$t_7=\underset{j=I}{\overset{N_g}{\Σ}}j\[\underset{i=1}{\overset{N_f}{\Σ}}H(i,j)\];$

灰度标准差：

$t_8={\{\underset{i=1}{\overset{N_f}{\Σ}}{(i-t_6)}^2\[\underset{j=1}{\overset{N_g}{\Σ}}H(i,j)\]\}}^{1/2};$

梯度标准差：

$t_9={\{\underset{j=1}{\overset{N_g}{\Σ}}{(j-t_7)}^2\[\underset{i=1}{\overset{N_f}{\Σ}}H(i,j)\]\}}^{1/2};$

相关性：

$t_{10}=\frac{1}{t_8{t}_9}\underset{i=i}{\overset{N_f}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}(i-t_6)(j-t_7)H(i,j);$

灰度熵：

$t_{11}=-\underset{i=1}{\overset{N_f}{\Σ}}\{\[\underset{j=1}{\overset{N_g}{\Σ}}H(i,j)\]{\log }_2\[\underset{j=1}{\overset{N_g}{\Σ}}H(i,j)\]\};$

梯度：

$t_{12}=-\underset{j=1}{\overset{N_g}{\Σ}}\{\[\underset{i=1}{\overset{N_f}{\Σ}}H(i,j)\]{\log }_2\[\underset{i=1}{\overset{N_f}{\Σ}}H(i,j)\]\};$

混合熵：

$t_{13}=-\underset{i=1}{\overset{N_f}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}H(i,j)\[{\log }_{2}H(i,j)\];$

差分矩：

$t_{14}=\underset{i=1}{\overset{N_f}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}{(i-j)}^{2}H(i,j);$

逆差矩：

$t_{15}=\underset{i=1}{\overset{N_f}{\Σ}}\underset{j=1}{\overset{N_g}{\Σ}}\frac{1}{1+{(i-j)}^2}H(i,j).$