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

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

Info

Publication number
CN106327480B
CN106327480B CN201610846458.8A CN201610846458A CN106327480B CN 106327480 B CN106327480 B CN 106327480B CN 201610846458 A CN201610846458 A CN 201610846458A CN 106327480 B CN106327480 B CN 106327480B
Authority
CN
China
Prior art keywords
image
thyroid
gray
gradient
value
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201610846458.8A
Other languages
English (en)
Other versions
CN106327480A (zh
Inventor
彭文献
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhejiang University ZJU
Original Assignee
Zhejiang Medical College
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhejiang Medical College filed Critical Zhejiang Medical College
Priority to CN201610846458.8A priority Critical patent/CN106327480B/zh
Publication of CN106327480A publication Critical patent/CN106327480A/zh
Application granted granted Critical
Publication of CN106327480B publication Critical patent/CN106327480B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/0002Inspection of images, e.g. flaw detection
    • G06T7/0012Biomedical image inspection
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10081Computed x-ray tomography [CT]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30004Biomedical image processing
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Quality & Reliability (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Health & Medical Sciences (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Processing (AREA)

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图像上移动,获得
Figure GDA0004156745200000051
种像素点对组合,Ng为图像灰度级别值,i∈[0,Ng-1],j∈[0,Ng-1],改变d和θ,θ为像素点对连线与水平线的夹角,统计甲状腺CT图像中的像素点对的灰度值,构成灰度共生矩阵P(i,j,θ,d),
Figure GDA0004156745200000052
其中#{x}是集合x中的所有元素的个数;然后根据灰度共生矩阵提取13维纹理特征,分别为角二阶矩、相关度、熵、对比度、逆差矩、平均和、和熵、方差和、方差、差平均、惯性、差方差、差熵。
作为优选,步骤二中所述的13维灰度共生纹理特征的公式为:
角二阶矩:
Figure GDA0004156745200000053
角二阶矩是灰度共生矩阵中所有像素的平方和,表现了图像灰度的分布均匀程度和纹理的粗细程度。
相关度:
Figure GDA0004156745200000054
其中,
Figure GDA0004156745200000061
Figure GDA0004156745200000062
Figure GDA0004156745200000063
相关度反映了灰度共生矩阵中的元素在水平方向和竖直方向的变化情况,当灰度共生矩阵中像素的元素分布均匀且数值接近,则相关度值较高。对于有方向性的纹理图像,如果特定θ方向的相关度较强,则可以反映了图像纹理的主要走向。
熵:
Figure GDA0004156745200000064
熵是灰度共生矩阵中信息量的表征,反映了灰度共生矩阵中元素随机性的程度,可以反映图像中纹理的复杂程度。当图像中所有像素分布均匀,没有任何纹理时,灰度共生矩阵的元素几乎全为0,则熵值较小;当图像中纹理非常复杂,像素分布杂乱时,灰度共生矩阵中的元素会存在较大的差异,则熵值较大。
对比度:
Figure GDA0004156745200000065
对比度是灰度共生矩阵相对主对角线的惯性矩。若灰度共生矩阵中的元素距离主对角线越远且元素值越大,则对比度较大,反之亦然。而距离主对角线越远,则图像中纹理基元的灰度差值越大;此类元素越多,则反映图像中此类纹理基元越明显。对比度可以反映图像中纹理基元强度,即沟纹的深浅程度。纹理基元表现越强,沟纹越深,则对比度越大;反之亦然。
逆差矩:
Figure GDA0004156745200000071
逆差矩也可称为局部平稳,可以反映图像纹理的同质性。当图像的纹理在不同区域变化较小时,图像的像素点对灰度更为接近,则逆差矩较大。
平均和:
Figure GDA0004156745200000072
其中
Figure GDA0004156745200000073
且|i+j|=k,k=2,3,…,2Ng
平均和反映了图像中像素点对的平均灰度值,反映了图像纹理的明暗程度。
和熵:
Figure GDA0004156745200000074
其中
Figure GDA0004156745200000075
且|i+j|=k,k=2,3,…,2Ng
和熵反映了图像中像素点对的分布变化,像素点对的变化越随机,分布范围越大,则和熵值越大。
方差和:
Figure GDA0004156745200000076
其中
Figure GDA0004156745200000081
且|i+j|=k,k=2,3,…,2Ng
方差和是图像中像素点对的变化程度,反映了图像纹理的周期性。图像中纹理的周期越大,则方差和越大。
方差:
Figure GDA0004156745200000082
其中m是灰度共生矩阵P(i,j,θ,d)所有元素的均值;
灰度共生矩阵的方差反映了图像像素点对的变化,图像中纹理分布越复杂,则方差越大。
差平均:
Figure GDA0004156745200000083
其中
Figure GDA0004156745200000084
且|i-j|=k,k=0,1,…,Ng-1;
差平均反映了图像的像素点对的平均灰度差异,如果图像纹理基元具有较强的灰度对比,并且分布空间较广,则差平均值较大。
惯性:
Figure GDA0004156745200000085
惯性描绘了图像纹理基元的对比度。如果图像中的像素点对具有较大的灰度差异,像素点对的数量较多,则惯性值较大。
差方差:
Figure GDA0004156745200000091
其中
Figure GDA0004156745200000092
且|i-j|=k,k=0,1,…,Ng-1;
差方差是图像中像素点对的灰度差的方差,描绘了图像纹理的对比度。当图像的纹理越显著,对比度越强强烈,则差方差越大。
差熵:
Figure GDA0004156745200000093
其中
Figure GDA0004156745200000094
且|i-j|=k,k=0,1,…,Ng-1。
作为优选,步骤二中获取所述的15维灰度梯度共生纹理特征过程为:对图像f求取梯度图像g,其中图像f大小为Nx×Ny,图像f的最高灰度级为N,灰度范围为[0,N-1],则有
Figure GDA0004156745200000095
其中gx和gy分别为边缘检测算子,对图像f在水平方向和竖直方向上计算得到梯度图像g,对灰度图像f和梯度图像g分别进行归一化处理,对应得到图像F和图像G,并且使得归一化后的图像G灰度分布范围为[0,Ng-1],图像F的灰度分布范围为[0,Nf-1],计算公式如下:
F=[f×(Nf-1)/fmax]+1,G=[g×(Ng-1)/gmax]+1,
其中fmax是灰度图像f中的灰度最大值,gmax是梯度图像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)/(Ng×Nf)
其中,i=1,2,…,Nf,j=1,2,…,Ng,根据灰度梯度共生矩阵提取15维纹理特征,分别为小梯度优势、大梯度优势、灰度分布的不均匀性、梯度分布的不均匀性、能量、灰度平均值、梯度平均值、灰度标准差、梯度标准差、相关性、灰度熵、梯度熵、混合熵、差分矩、逆差矩。对于灰度图像,不仅可以使用灰度信息统计分析图像的纹理特征,而且可以使用灰度的邻域变化信息分析图像纹理。梯度就是灰度的邻域变化信息,构成了图像的边缘轮廓。灰度梯度共生矩阵可以综合图像的灰度和梯度信息,更为全面的反映图像中每个像素点的灰度和梯度的分布特征,反映像素点与其邻域各像素点的空间位置关系。
作为优选,步骤二中所述的15维灰度梯度共生纹理特征的公式为:
小梯度优势:
Figure GDA0004156745200000101
其中
Figure GDA0004156745200000102
小梯度优势反映了图像中小梯度的分布强度。当图像的灰度变化较平缓,图像大部分区域的梯度值较小,则小梯度优势值大。
大梯度优势:
Figure GDA0004156745200000103
其中
Figure GDA0004156745200000104
大梯度优势反映了图像中的大梯度的分布强度。如果图像中的灰度变化非常剧烈,大梯度在图像中分布较广泛,则大梯度优势值大。
灰度分布的不均匀性:
Figure GDA0004156745200000111
其中
Figure GDA0004156745200000112
如果图像中不同区域的灰度差别比较大,则灰度分布的不均匀性值较大。
梯度分布的不均匀性:
Figure GDA0004156745200000113
其中
Figure GDA0004156745200000114
如果图像中不同区域的梯度变化较大,则梯度分布的不均匀性值较大。
能量:
Figure GDA0004156745200000115
能量反映了图像纹理的强度和密集程度。如果图像中灰度和梯度变化比较强烈,且较为密集,则能量值较大。
灰度平均值:
Figure GDA0004156745200000116
灰度平均值反映了图像中的灰度强度,即图像的整体亮度。
梯度平均值:
Figure GDA0004156745200000117
梯度平均值反映了图像中的梯度的平均强度。
灰度标准差:
Figure GDA0004156745200000121
灰度标准差反映了图像中灰度的变化强度
梯度标准差:
Figure GDA0004156745200000122
梯度标准差反映了图像中梯度的变化强度。如果图像中的梯度大小较为接近,则梯度标准差值较小,反之亦然。
相关性:
Figure GDA0004156745200000123
相关性反映了图像中灰度和梯度的相互关系。如果图像中纹理基元的分布比较规则,那么灰度和梯度的变化存在一致性,则相关性值较大。
灰度熵:
Figure GDA0004156745200000124
灰度熵反映了图像中灰度分布的规律性,即图像中的平均信息量。
梯度熵:
Figure GDA0004156745200000125
梯度熵反映了图像中梯度分布的规律性,即梯度图像的平均信息量。
混合熵:
Figure GDA0004156745200000131
混合熵反映了灰度-梯度共生矩阵的平均信息量。如果图像中灰度和梯度的变化较复杂,则混合熵值较大。
差分矩:
Figure GDA0004156745200000132
差分矩描绘了图像纹理基元的对比度。如果图像中像素点对的灰度较高并且灰度变化很小,那么图像的差分矩较大。
逆差矩:
Figure GDA0004156745200000133
逆差矩,即局部平稳,可以反映图像纹理的同质性。当图像纹理基元的灰度与梯度非常接近,并且覆盖图像的大部分区域,则逆差矩较大。
本发明中,步骤一中对甲状腺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图像上移动,获得
Figure GDA0004156745200000171
种像素点对组合,Ng为图像灰度级别值,i∈[0,Ng-1],j∈[0,Ng-1],改变d和θ,θ为像素点对连线与水平线的夹角,统计甲状腺CT图像中的像素点对的灰度值,构成灰度共生矩阵P(i,j,θ,d),
Figure GDA0004156745200000181
其中#{x}是集合x中的所有元素的个数;预先设定d值,分别计算0°、45°、90°和135°四个方向的灰度共生矩阵,如果灰度范围为[0,Ng-1],则每个方向的灰度共生矩阵的尺寸是Ng×Ng。本实施例中,预先设定d=1,分别计算0°、45°、90°和135°的灰度共生矩阵,然后基于每个方向的矩阵计算纹理特征,将四个方向的纹理特征求取平均值,得到旋转不变的纹理特征。然后根据灰度共生矩阵提取13维纹理特征,分别为角二阶矩、相关度、熵、对比度、逆差矩、平均和、和熵、方差和、方差、差平均、惯性、差方差、差熵。
本实施例中,步骤二中所述的13维灰度共生纹理特征的公式为:
角二阶矩:
Figure GDA0004156745200000182
相关度:
Figure GDA0004156745200000183
其中,
Figure GDA0004156745200000184
Figure GDA0004156745200000185
Figure GDA0004156745200000186
熵:
Figure GDA0004156745200000191
对比度:
Figure GDA0004156745200000192
逆差矩:
Figure GDA0004156745200000193
平均和:
Figure GDA0004156745200000194
其中
Figure GDA0004156745200000195
且|i+j|=k,k=2,3,…,2Ng
和熵:
Figure GDA0004156745200000196
其中
Figure GDA0004156745200000197
且|i+j|=k,k=2,3,…,2Ng
方差和:
Figure GDA0004156745200000198
其中
Figure GDA0004156745200000201
且|i+j|=k,k=2,3,…,2Ng
方差:
Figure GDA0004156745200000202
其中m是灰度共生矩阵P(i,j,θ,d)所有元素的均值;
差平均:
Figure GDA0004156745200000203
其中
Figure GDA0004156745200000204
且|i-j|=k,k=0,1,…,Ng-1;
惯性:
Figure GDA0004156745200000205
差方差:
Figure GDA0004156745200000206
其中
Figure GDA0004156745200000207
且|i-j|=k,k=0,1,…,Ng-1;
差熵:
Figure GDA0004156745200000208
其中
Figure GDA0004156745200000211
且|i-j|=k,k=0,1,…,Ng-1。
本实施例中,步骤二中获取所述的15维灰度梯度共生纹理特征过程为:对图像f求取梯度图像g,其中图像f大小为Nx×Ny,图像f的最高灰度级为N,灰度范围为[0,N-1],则有
Figure GDA0004156745200000212
其中gx和gy分别为边缘检测算子,本实施例中采用尺寸为3*3的Sobel算子求取梯度图像,Sobel算子是水平方向和竖直方向的两个边界检测算子,其分别为:
Figure GDA0004156745200000213
对图像f在水平方向和竖直方向上计算得到梯度图像g,对灰度图像f和梯度图像g分别进行归一化处理,对应得到图像F和图像G,并且使得归一化后的图像G灰度分布范围为[0,Ng-1],图像F的灰度分布范围为[0,Nf-1],计算公式如下:
F=[f×(Nf-1)/fmax]+1,G=[g×(Ng-1)/gmax]+1,
其中fmax是灰度图像f中的灰度最大值,gmax是梯度图像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)/(Ng×Nf)
其中,i=1,2,…,Nf,j=1,2,…,Ng,根据灰度梯度共生矩阵提取15维纹理特征,分别为小梯度优势、大梯度优势、灰度分布的不均匀性、梯度分布的不均匀性、能量、灰度平均值、梯度平均值、灰度标准差、梯度标准差、相关性、灰度熵、梯度熵、混合熵、差分矩、逆差矩。
本实施例中,步骤二中所述的15维灰度梯度共生纹理特征的公式为:
小梯度优势:
Figure GDA0004156745200000221
其中
Figure GDA0004156745200000222
大梯度优势:
Figure GDA0004156745200000223
其中
Figure GDA0004156745200000224
灰度分布的不均匀性:
Figure GDA0004156745200000225
其中
Figure GDA0004156745200000226
梯度分布的不均匀性:
Figure GDA0004156745200000227
其中
Figure GDA0004156745200000228
能量:
Figure GDA0004156745200000231
灰度平均值:
Figure GDA0004156745200000232
梯度平均值:
Figure GDA0004156745200000233
灰度标准差:
Figure GDA0004156745200000234
梯度标准差:
Figure GDA0004156745200000235
相关性:
Figure GDA0004156745200000236
灰度熵:
Figure GDA0004156745200000237
梯度熵:
Figure GDA0004156745200000241
混合熵:
Figure GDA0004156745200000242
差分矩:
Figure GDA0004156745200000243
逆差矩:
Figure GDA0004156745200000244
通过以上的纹理特征提取,得到28维特征参数,参数总结如下表所示:
Figure GDA0004156745200000245
总之,以上所述仅为本发明的较佳实施例,凡依本发明申请专利范围所作的均等变化与修饰,皆应属本发明专利的涵盖范围。

Claims (9)

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值异常的判断结果;
采用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,优化计算运行结束。
2.根据权利要求1所述的一种甲状腺CT图像异常密度的检测方法,其特征在于:最优高密度阈值和最优低密度阈值分别为140-160HU和40-60HU。
3.根据权利要求1所述的一种甲状腺CT图像异常密度的检测方法,其特征在于:最优高密度阈值和最优低密度阈值分别为145-155HU和45-55HU。
4.根据权利要求1所述的一种甲状腺CT图像异常密度的检测方法,其特征在于:最优高密度阈值和最优低密度阈值分别为145-150HU和45-50HU。
5.根据权利要求1所述的一种甲状腺CT图像异常密度的检测方法,其特征在于:最优高密度阈值为140HU或145HU或150HU或155HU或160HU,最优低密度阈值为40HU或45HU或50HU或55HU或60HU。
6.根据权利要求1或3或5所述的一种甲状腺CT图像异常密度的检测方法,其特征在于:步骤二中获取所述的13维灰度共生纹理特征过程为:
在甲状腺CT图像中任取一像素点A(x,y),然后获得与像素点A距离为d的另一个像素点B(x+a,y+b),像素点A和B形成一个点对,记录这一点对的灰度值(i,j),固定a和b的值,使像素点A(x,y)在甲状腺CT图像上移动,获得
Figure QLYQS_1
种像素点对组合,Ng为图像灰度级别值,i∈[0,Ng-1],j∈[0,Ng-1],改变d和θ,θ为像素点对连线与水平线的夹角,统计甲状腺CT图像中的像素点对的灰度值,构成灰度共生矩阵P(i,j,θ,d),
Figure QLYQS_2
其中#{X}是集合X中的所有元素的个数;然后根据灰度共生矩阵提取13维纹理特征,分别为角二阶矩、相关度、熵、对比度、逆差矩、平均和、和熵、方差和、方差、差平均、惯性、差方差、差熵。
7.根据权利要求6所述的一种甲状腺CT图像异常密度的检测方法,其特征在于:步骤二中所述的13维灰度共生纹理特征的公式为:
角二阶矩:
Figure QLYQS_3
相关度:
Figure QLYQS_4
其中,
Figure QLYQS_5
Figure QLYQS_6
Figure QLYQS_7
熵:
Figure QLYQS_8
对比度:
Figure QLYQS_9
逆差矩:
Figure QLYQS_10
平均和:
Figure QLYQS_11
其中
Figure QLYQS_12
且|i+j|=k,k=2,3,…,2Ng
和熵:
Figure QLYQS_13
其中
Figure QLYQS_14
且|i+j|=k,k=2,3,…,2Ng
方差和:
Figure QLYQS_15
其中
Figure QLYQS_16
且|i+j|=k,k=2,3,…,2Ng
方差:
Figure QLYQS_17
其中m是灰度共生矩阵P(i,j,θ,d)所有元素的均值;
差平均:
Figure QLYQS_18
其中
Figure QLYQS_19
且|i-j|=k,k=0,1,…,Ng-1;
惯性:
Figure QLYQS_20
差方差:
Figure QLYQS_21
其中
Figure QLYQS_22
且|i-j|=k,k=0,1,…,Ng-1;
差熵:
Figure QLYQS_23
其中
Figure QLYQS_24
且|i-j|=k,k=0,1,…,Ng-1。
8.根据权利要求1或3或5或7所述的一种甲状腺CT图像异常密度的检测方法,其特征在于:步骤二中获取所述的15维灰度梯度共生纹理特征过程为:对图像f求取梯度图像g,其中图像f大小为Nx×Ny,图像f的最高灰度级为N,灰度范围为[0,N-1],则有
Figure QLYQS_25
其中gx和gy分别为边缘检测算子,对图像f在水平方向和竖直方向上计算得到梯度图像g,对灰度图像f和梯度图像g分别进行归一化处理,对应得到图像F和图像G,并且使得归一化后的图像G灰度分布范围为[0,Ng-1],图像F的灰度分布范围为[0,Nf-1],计算公式如下:
F=[f×(Nf-1)/fmax]+1,G=[g×(Ng-1)/gmax]+1,
其中fmax是灰度图像f中的灰度最大值,gmax是梯度图像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)/(Ng×Nf)
其中,i=1,2,…,Nf,j=1,2,…,Ng,根据灰度梯度共生矩阵提取15维纹理特征,分别为小梯度优势、大梯度优势、灰度分布的不均匀性、梯度分布的不均匀性、能量、灰度平均值、梯度平均值、灰度标准差、梯度标准差、相关性、灰度熵、梯度熵、混合熵、差分矩、逆差矩。
9.根据权利要求8所述的一种甲状腺CT图像异常密度的检测方法,其特征在于:步骤二中所述的15维灰度梯度共生纹理特征的公式为:
小梯度优势:
Figure QLYQS_26
其中
Figure QLYQS_27
大梯度优势:
Figure QLYQS_28
其中
Figure QLYQS_29
灰度分布的不均匀性:
Figure QLYQS_30
其中
Figure QLYQS_31
梯度分布的不均匀性:
Figure QLYQS_32
其中
Figure QLYQS_33
能量:
Figure QLYQS_34
灰度平均值:
Figure QLYQS_35
梯度平均值:
Figure QLYQS_36
灰度标准差:
Figure QLYQS_37
梯度标准差:
Figure QLYQS_38
相关性:
Figure QLYQS_39
灰度熵:
Figure QLYQS_40
梯度熵:
Figure QLYQS_41
混合熵:
Figure QLYQS_42
差分矩:
Figure QLYQS_43
逆差矩:
Figure QLYQS_44
CN201610846458.8A 2016-09-23 2016-09-23 一种甲状腺ct图像异常密度的检测方法 Expired - Fee Related CN106327480B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610846458.8A CN106327480B (zh) 2016-09-23 2016-09-23 一种甲状腺ct图像异常密度的检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610846458.8A CN106327480B (zh) 2016-09-23 2016-09-23 一种甲状腺ct图像异常密度的检测方法

Publications (2)

Publication Number Publication Date
CN106327480A CN106327480A (zh) 2017-01-11
CN106327480B true CN106327480B (zh) 2023-06-27

Family

ID=57819986

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610846458.8A Expired - Fee Related CN106327480B (zh) 2016-09-23 2016-09-23 一种甲状腺ct图像异常密度的检测方法

Country Status (1)

Country Link
CN (1) CN106327480B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107633519A (zh) * 2017-09-27 2018-01-26 合肥美亚光电技术股份有限公司 单调性图像特征辨别系统的参数异常诊断方法及系统
CN108831532B (zh) * 2018-06-15 2022-03-29 北京大学第一医院 一种核医学甲状腺显像图像处理方法及系统
CN111657985A (zh) * 2020-06-30 2020-09-15 杭州依图医疗技术有限公司 肺部影像信息的处理方法、显示方法及可读存储介质
CN112465824B (zh) * 2021-01-28 2021-08-03 之江实验室 基于pet/ct图像亚区影像组学特征的肺腺鳞癌诊断装置
CN117078671B (zh) * 2023-10-13 2023-12-12 陕西秒康医疗科技有限公司 一种甲状腺超声影像智能分析系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102419864A (zh) * 2011-09-05 2012-04-18 东软集团股份有限公司 一种提取脑部ct图像骨骼方法及装置
CN104000619A (zh) * 2014-06-16 2014-08-27 彭文献 一种甲状腺ct图像计算机辅助诊断系统及方法
JP2015226274A (ja) * 2014-05-29 2015-12-14 京セラドキュメントソリューションズ株式会社 画像読取装置、画像形成装置、画像処理方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102419864A (zh) * 2011-09-05 2012-04-18 东软集团股份有限公司 一种提取脑部ct图像骨骼方法及装置
JP2015226274A (ja) * 2014-05-29 2015-12-14 京セラドキュメントソリューションズ株式会社 画像読取装置、画像形成装置、画像処理方法
CN104000619A (zh) * 2014-06-16 2014-08-27 彭文献 一种甲状腺ct图像计算机辅助诊断系统及方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Xia S ."A novel computational ct image analysis method for classifying nodules from normal thyroid tissue".《Iet International Conference on Biomedical Image &amp Signal Processing》.2016,第1-6页. *
耿欢."基于CT 影像的肺组织分割方法综述".《计算机应用研究》.2016,(第7期),第1-7页. *

Also Published As

Publication number Publication date
CN106327480A (zh) 2017-01-11

Similar Documents

Publication Publication Date Title
CN106327480B (zh) 一种甲状腺ct图像异常密度的检测方法
US20240096048A1 (en) System and Method for the Visualization and Characterization of Objects in Images
Chang et al. Computer‐aided diagnosis for classifying benign versus malignant thyroid nodules based on ultrasound images: a comparison with radiologist‐based assessments
Zhang et al. Dual-mode artificially-intelligent diagnosis of breast tumours in shear-wave elastography and B-mode ultrasound using deep polynomial networks
US11669964B2 (en) Systems and methods to facilitate review of liver tumor cases
Banerjee et al. Automated 3D segmentation of brain tumor using visual saliency
Khan et al. J Pathol Inform
TWI552013B (zh) 對患者組織超音波衰減影像中可疑關注區域分類之方法、裝置和電腦程式產品
Garnavi et al. Automatic segmentation of dermoscopy images using histogram thresholding on optimal color channels
Abbas et al. Skin tumor area extraction using an improved dynamic programming approach
Khan et al. HyMaP: A hybrid magnitude-phase approach to unsupervised segmentation of tumor areas in breast cancer histology images
WO2014113786A1 (en) Quantitative predictors of tumor severity
JP2007007440A (ja) 医療画像において腫瘤や実質組織変形をコンピュータを用いて検出する自動化した方法と装置
Belkacem-Boussaid et al. Automatic detection of follicular regions in H&E images using iterative shape index
US20140140629A1 (en) Methods for processing target pattern, method for generating classification system of target patterns and method for classifying detected target patterns
Kaur et al. Computer-aided diagnosis of renal lesions in CT images: a comprehensive survey and future prospects
CN111127404A (zh) 一种医疗影像轮廓快速提取方法
Beheshti et al. Classification of abnormalities in mammograms by new asymmetric fractal features
Janan et al. RICE: A method for quantitative mammographic image enhancement
CN113689424B (zh) 可自动识别图像特征的超声检查系统及识别方法
Arikidis et al. A two‐stage method for microcalcification cluster segmentation in mammography by deformable models
US20160100789A1 (en) Computer-aided diagnosis system and computer-aided diagnosis method
Huang et al. Neural network analysis applied to tumor segmentation on 3D breast ultrasound images
Chang et al. Quantitative analysis of breast echotexture patterns in automated breast ultrasound images
Gill et al. Automatic region growing segmentation for medical ultrasound images

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20230627