基于双视角信息融合的乳腺CR图像的计算机确定方法 【技术领域】
本发明属于医学影像的计算机辅助分析技术的应用领域，具体涉及到一种基于双视角信息融合的乳腺CR图像的计算机确定方法。
背景技术
乳腺癌是目前世界上女性发病率最高的恶性肿瘤之一。据统计，全球每年约有120万妇女被确诊患有乳腺癌，约有50万妇女死于乳腺癌。乳腺癌以前主要是在北美、欧洲和大洋洲地区比较高发，据世界卫生组织(World HealthOrganization，WHO)估计，欧美发达国家平均不到十位女性中会有一位患上乳腺癌。但近些年来，乳腺癌在我国的发病率节节攀升，每年的增长率达到3％，发病年龄也呈逐渐年轻化的趋势，较西方国家大约要提前10至15年。在北京、天津、上海等经济较发达的城市，乳腺癌已经超过肺癌，成为女性恶性肿瘤中的头号杀手，因此乳腺癌的诊断与治疗越来越受到政府与国民的关注。
尽管乳腺癌的死亡率很高，但如果能尽早发现，治愈率将高达80～90％，到晚期就只剩下40％左右的治愈率。而实际中，80％的患者确诊时已经是癌症中晚期了，为避免细胞转移，医生大都采取手术切除，这样，不仅增加了患者的生理痛苦、经济负担，还将使患者忍受缺失乳房的终生遗憾。因此，许多专家呼吁，女性要想在症状显现之前就能检测到疾病的存在，在做好自我检查的基础上，还需要定期到医院进行乳腺癌检测。
在乳腺癌检测方面虽然存在多种方法，但目前而言，国内外都将乳腺钼靶X线摄影技术当作常规的检测手段，它能较清晰地表现乳腺癌的征象——微钙化与肿块，对乳腺癌的早期检测效果较好。从05年开始，由卫生部组织实施的“百万妇女乳腺普查工程”在中国全面铺开，大批适龄妇女得到了免费的乳腺癌检查，乳腺钼靶摄影也是其重要的检测手段。为了让乳腺钼靶图像的处理、保存与查询更加方便快捷，并融入各级医疗机构所使用的HIS/PACS(HospitalInformation System/Picture Archiving and Communication System)系统之中，需要将钼靶图像数字化，国内通常采用计算机X线成像(ComputedRadiography，CR)技术来实现。目前，国内许多医院的放射科医师都是面对乳腺钼靶CR图像进行乳腺癌筛查与诊断工作的。
由于乳腺钼靶图像的阅读分析依赖于放射科医师的专业技术和经验，对医师素质的要求很高，而在普查或体检过程中往往会产生海量的钼靶摄片，其倍增的阅读工作量将让各级医疗机构的放射科不堪重负；同时，病灶有时会隐藏在各种组织结构内部，并非总是显而易见，如果让医师们长时间读片，可能会产生瞬间疏忽和视觉疲劳，根据美国相关医疗机构的统计，医师们由于瞬间疏忽和视觉疲劳而遗漏的病灶可占总病灶数的10～30％。
为了提高乳腺钼靶CR图像的阅读效率，通常可以采用计算机辅助诊断(Computer-aided Diagnosis，CAD)系统来帮助放射科医师。CAD系统并不是要创造出一个技艺高超的机器医师，而只是希望为医师提供辅助意见、减轻医师的工作量以及提高诊断的精确性。它会采用一系列的计算机视觉与机器学习方法，将乳腺钼靶图像中的疑似病灶(包括微钙化和肿块等)自动检测出来，作为一种初筛结果或者第二观感提供给医师，为医师的最终诊断提供参考。构建医学影像特别是乳腺CR图像的CAD系统在现代的医疗保健模式中已经显得越来越重要，成为国内外研究与开发的热点。
在传统的乳腺CAD系统中，自动检测通常基于当前图像的局部信息，故而当病灶在当前图像中不够明显时，检测效果就无法达到很好。而如果能将同一乳房两个视角(MLO：内外侧斜位；CC：头尾轴位)的图像信息结合起来，对同一病灶在不同图像中的特征进行融合，就能在病灶不够明显的情况下依旧取得不错的检测效果，这就是本发明的研究重点。
【发明内容】
本发明的目的在于提供一种基于双视角信息融合的乳腺CR图像的计算机确定方法。该方法能将同一乳房两个视角的乳腺CR图像的信息结合起来，通过计算机视觉、机器学习以及特征融合技术，将那些隐藏在乳腺组织结构内部的疑似病灶区域检测出来，提供给放射科医师作参考。
本发明方法的具体步骤是：
步骤(1)读入一幅待分析的乳腺CR图像；
步骤(2)使用图像分割技术(Snake算法)提取乳房区域；
步骤(3)在乳房区域内寻找疑似肿块可能存在的位置；
步骤(4)将乳房区域内的疑似肿块进行手动分割；
步骤(5)提取已分割地疑似肿块的特征(包括紧凑度，偏心率，矩特征，Gabor特征、分维数)，用于鉴定疑似肿块的真伪；
步骤(6)在乳房区域内寻找疑似微钙化可能存在的位置；
步骤(7)使用分水岭算法对乳房区域内的疑似微钙化进行分割；
步骤(8)提取已分割的疑似微钙化的特征(包括紧凑度，偏心率，矩特征，Gabor特征、分维数)，用于鉴定疑似微钙化的真伪；
步骤(9)确定乳腺钼靶X线摄影的两个视角中同一病灶的图像元素间的匹配关系；
步骤(10)对两个视角中同一病灶的特征信息进行融合(Multi-Agent方法)，确定病灶良恶性；
步骤(11)生成计算机辅助诊断系统检测结果报告。
本发明所提供的基于双视角信息融合的乳腺CR图像计算机辅助诊断方法由一组功能模块构成，它们包括：图像管理模块、乳房区域提取模块、肿块检测模块、微钙化簇检测模块、双视角病灶匹配模块、病灶决策模块和计算机辅助诊断报告模块。其中：
图像管理模块负责通过磁盘、网络或其他应用程序(包括HIS/PACS)的接口读取或存储乳腺钼靶CR图像数据，查询、登记或更新与乳腺CR图像所对应的病历数据。
乳房区域提取模块负责将读入的乳腺CR图像中的乳房区域提取出来，将其他的标识信息等都剔除。
肿块检测模块分成疑似肿块搜索模块，疑似肿块分割模块和疑似肿块鉴定模块三个子模块：
疑似肿块搜索模块负责在整个乳房区域内寻找疑似肿块可能存在的位置。
疑似肿块分割模块负责将已知大致部位的疑似肿块从乳房组织背景中提取出来。
疑似肿块鉴定模块负责鉴定已分割的疑似肿块，剔除假阳性。
微钙化簇检测模块分成疑似微钙化搜索模块、疑似微钙化分割模块、疑似微钙化鉴定模块和微钙化聚簇模块四个子模块：
疑似微钙化搜索模块负责在整个乳房区域内寻找疑似微钙化可能存在的位置。
疑似微钙化分割模块负责将已知大致部位的疑似微钙化从乳房组织背景中提取出来。
疑似微钙化鉴定模块负责鉴定已分割的疑似微钙化，剔除假阳性。
微钙化聚簇模块负责将原本孤立的微钙化根据彼此的相互距离聚合起来，形成微钙化簇。
双视角病灶匹配模块分成双视角肿块匹配模块和双视角微钙化簇匹配模块两个子模块：
双视角肿块匹配模块负责确定两个视角中同一肿块的图像元素之间的匹配关系。
双视角微钙化簇匹配模块负责确定两个视角中同一微钙化簇的图像元素之间的匹配关系。
病灶决策模块分成基于双视角信息融合的肿块决策模块和基于双视角信息融合的微钙化簇决策模块两个子模块：
基于双视角信息融合的肿块决策模块负责将同一肿块在两个视角中的特征信息进行融合，确定肿块的良恶性。
基于双视角信息融合的微钙化簇决策模块负责将同一微钙化簇在两个视角中的特征信息进行融合，确定微钙化簇的良恶性。
计算机辅助诊断报告模块负责将计算机辅助诊断系统的检测结果制作成报告，呈报给放射科医师参考。
采用本发明的系统框架与技术方案，能在没有放射科医师监督的情况下自动地对乳腺CR图像中的乳房区域进行检测，搜索并分割出疑似的病灶(包括肿块与微钙化簇)区域，并将双视角的病灶信息进行融合，最终判断出病灶的良恶性，作为初筛结果或者第二观感提供给医师作参考，从而有效地提高基于乳腺CR图像的乳腺癌筛查工作的效率与精度，减轻相应的工作负荷。目前我国各级医疗机构中乳腺CR机覆盖面很广，本发明所提供的基于双视角信息融合的乳腺CR图像计算机辅助诊断系统将CAD技术与乳腺CR图像的特点相结合，能快速准确地对疑似病灶进行定位与鉴定，对于技术还有待改善的实习医师来说，是摄片阅读诊断的良好辅助，对于经验丰富的放射科医师来说，也能提供很有裨益的参考意见。
【具体实施方式】
一种基于双视角信息融合的乳腺CR图像计算机确定方法的具体工作步骤与实施方式如下：
(1)读入一幅待分析的乳腺CR图像。
(2)使用一种基于最优阈值与分水岭的分割技术，将乳房区域分离出来，并剔除剩余的标识区域以及背景中的非零区域。
(2.a)第一步是剔除背景中灰度大于零且又明显小于乳房与标识区域的像素，采用的方法是最优阈值法，这是一种迭代方法。首先计算乳腺X线图像的平均灰度作为背景阈值的初始值：
T0=1NΣ0≤k≤Mk×h(k),]]>N=Σ0≤k≤Mh(k)---(1)]]>
T0表示初始阈值，k表示灰度值，h(k)表示灰度为k的象素的个数(即灰度直方图中横坐标为k时的纵坐标取值)，N表示图像中象素的总数，M表示图像的最大灰度级减去1。
然后，分别计算所有灰度高于该阈值的象素的平均灰度和所有灰度不高于该阈值的象素的平均灰度，并求其均值。
Tu=1Σk>T0h(k)Σk>T0k×h(k),]]>Td=1Σk≤T0h(k)Σk≤T0k×h(k),]]>T1＝(Tu+Td)/2 (2)
如果T0和T1非常接近或者相等，则认为T1即可作为迭代阈值的最优结果；如果T0和T1有较大差距，则将T0的取值更新成T1，继续计算公式(2)，继续判断T0和T1是否接近，如此迭代，直至得到最优结果。计算得到最优阈值之后，就可以对全图实施一次向下阈值化(即将灰度值小于等于阈值的象素灰度置0，其他象素不变)，剔除背景中的非零区域。
(2.b)第二步是使用经典的种子填充方法，将图像中剩余的连通域(包括乳房、边框和标识等区域)分别提取出来，由于乳房区域的面积远远大于其他区域，因此根据面积进行筛选，就可以剔除其他区域，只留下乳房区域。
(2.c)第三步是使用经典的分水岭方法，将乳房区域的粗分割边缘进行精细化。先利用一阶差分算子计算出乳房区域在垂直与水平方向上的梯度，生成相应的梯度模值图。在该图中，乳房区域内部和背景区域是两个很大的聚水盆地，它们的海拔高度是0，处于两者之间的就是边缘带中的分水岭(其海拔高度就是梯度模值)。然后，往两个聚水盆地中注水，使其水面抬升。这时，两个聚水盆地的表面积随着水面的抬升而增大，向分水岭漫堤。当分别来自两个聚水盆地的水将要突破某处的分水岭堤坝而汇聚在一起时，就在此处建立宽度为一个象素的永久性堤坝，即分水线。随着两个聚水盆地的水不断上涨和趋向汇聚，最后，整条边缘带分水岭都变成了永久性堤坝，这条线就是区域的分界，也就是乳房区域的精细化边界。
(3)使用一种双模式定位方法，在乳房区域内寻找肿块可能存在的位置。
隐藏在乳房区域内部的肿块主要可以分成两种类型，一种是处于脂肪组织区域或稀疏的腺体组织区域的肿块，它相对于周围的背景(脂肪或稀疏腺体)来说比较明显，对比度较高；另一种是藏于致密的腺体组织内部的肿块，它与周围的背景(致密腺体)差别不够明显，对比度较低。对于这两种特性不同的肿块，不宜采用统一的标准进行检测，因此使用双模式定位方法。
(3.a)使用一种迭代的脂肪组织剥离技术来提取乳房区域内部的高灰度区域。在迭代的每一步中，步骤(2.a)中介绍的最优阈值法和被用来计算剥离阈值，通过向下阈值化，灰度低于剥离阈值的脂肪组织像素被剔除，留下高灰度的腺体组织和肿块区域的像素；经过几步剥离，这些高灰度区域之间的脂肪组织像素就被大批剥离，使得残留区域之间相互分离；如果一次剥离过程从某个残留区域中所剔除的像素与该区域面积的比值低于一个经验阈值(说明脂肪组织像素已经被剥离殆尽)，该区域就被视作已成功分离，此后不再接受剥离。在这些残留区域中，面积最大的是腺体组织区域(因腺体彼此相连，分布较广)，而其他则是处在脂肪组织区域的局部高灰度团块。上面介绍的第一类肿块如果存在，通常就位于这些局部高灰度团块之中，通过这种模式的检测方法即可完成定位；而第二类肿块则往往深藏在腺体组织区域内部，需要使用第二种模式的定位方法。
(3.b)使用一种基于离散小波分解的高频零模值区域(黑洞)检测方法来搜索腺体组织区域内的第二类肿块。小波分解是一种经典的频域分析技术，它具有很好的光滑性和局部性，能将信号根据分辨率逐层分解成高频段与低频段，从而进行多分辨率分析(Multi-resolution Analysis，MRA)。在数字图像领域，信号是离散二维的，因此使用的是二维的离散小波变换(Discrete WaveletTransform，DWT)。利用二维DWT技术，腺体组织区域的信息将被逐层分解(每层分解成低频子带LL和三个高频子带LH、HL、HH，通常将三个高频子带模值对位累加，构成一个统一的高频合并带)，该区域内部的灰度起伏会在高频合并带的对应位置呈现出一定模值的信号。
而在第二类肿块的邻域环境中，虽然肿块与腺体组织之间不存在很高的对比度，但腺体的纹理丰富，灰度起伏较大，而肿块内部致密，灰度一致性较高。因此对应到高频合并带中，肿块所在的区域小波系数模值极低，而周围的腺体区域的小波系数相对较高，通过在小波分解的高频合并带中搜索零模值区域(黑洞)可以对疑似的第二类肿块进行定位。由于在高频合并带中存在许多孤立的零模值点，可以先使用中值滤波技术对高频合并带进行处理，抑制孤立的零模值点，然后采用经典的区域生长技术在高频合并带中提取零模值区域，并根据一个经验阈值，将面积太小的零模值区域剔除，只留下大块的零模值区域，作为疑似肿块进行重建。由于小波分解第一层解析出来的是最细微的高频信号，而肿块的尺寸一般不会太小，因此零模值检测通常在小波分解的第三、四层进行，超过四层之后，对应区块尺寸过大，也不宜用来检测肿块。
(4)使用一种基于高斯金字塔的多分辨率梯度矢量流(Gradient VectorFlow，GVF)Snake技术，将已知大致位置的疑似肿块从乳房背景中提取出来。
高斯金字塔是对图像进行高斯卷积和下采样得到的一组分辨率由高到低(低分辨率到高分辨率)的图像。在低分辨率下，图像尺寸小且粗糙，噪声较少；在高分辨率下，图像尺寸大而精细，噪声、伪影较多。GVF Snake是经典的Snake分割模型的改进，将GVF作为Snake模型中的外力：
E=∫∫μ(ux2+uy2+vx2+vy2)+|▿f|2|v-▿f|2dxdy---(3)]]>
式中为拉普拉斯算子，f为边界图，v为梯度矢量流场。
多分辨率GVF Snake算法首先对图像进行高斯金字塔分解。在低分辨率图像上设置初始轮廓作为GVF Snake的初始条件，并进行迭代，由于此时尺寸较小且噪声较少，可以快速得到一个相对接近真实目标的轮廓，将该结果轮廓作为较高分辨率图像上GVF Snake的初始条件，可进一步得到一个更准确的轮廓，依次类推，最终可在最高分辨率图像(原图像)上获取较为理想的目标轮廓，完成疑似肿块的分割。
(5)从已分割出来的疑似肿块中抽取一组灰度、形状与纹理特征，然后使用一种神经网络分类器来鉴定疑似肿块的真伪。
本发明所采用的神经网络分类器是多层感知器网络(Multi-LayerPerceptrons Networks，MLPN)，这是一种在模式识别领域很常用的技术，采用多层前馈的网络结构，使用反向传播(Back Propagation，BP)技术进行训练，因此也常被称作BP神经网络。所用的MLPN采用了三层网络结构，隐层节点设置成12个，输入节点数就是特征的个数，输出节点为1个，其理想输出为：1为真，0为假。而所采用的特征组则共有8种特征，涉及到灰度、形状和纹理等，它们包括：区域的平均灰度、区域的灰度标准差、边缘带的对比度(肿块边缘带中，属于肿块区域的像素的平均灰度与外周像素的平均灰度的差值)、边界的梯度方向熵、区域面积、区域的紧凑度(面积与周长的比值)、区域的狭长度(最长的肿块区域穿心线段与其垂直方向的肿块区域穿心线段之间的比值)以及区域的Hausdorff分形维数(通过经典的盒计数法求解)。
这里提到的肿块边缘带是指在肿块边界线两侧有一定宽度的过渡带，通过将二值化的肿块区域分别实施一定次数的形态学膨胀和腐蚀，然后用膨胀的结果减去腐蚀的结果即可得到；肿块边界梯度方向熵则是对于肿块边界线上的每个点，先利用一阶差分算子计算出垂直梯度与水平梯度，然后融合得到梯度矢量的方向，并将其规整到一系列标准化方向(以肿块区域的几何中心为起点，从零角度开始，每隔3度向外作一条射线，共有120条这样的射线，每条射线的方向就是标准化方向)上，最后根据标准化方向i上的边界点分布比例Pi计算熵值，即：
Pi=NiM,]]>M=Σi=0119Ni,]]>E=-Σi=0119PilogPi---(4)]]>
当从肿块区域几何中心发出的每条射线与肿块边界线相交时，交点与起点之间的线段就是中心射线段，它可以用来度量每个标准化方向上从肿块中心到边界的距离，而将方向相反(间隔180度，即相隔60条射线)的两条中心射线段连接起来，就形成了一条穿过肿块区域几何中心且分别到达两个边界的线段，也就是计算肿块狭长度时所提到的肿块区域穿心线段，这样的线段共有60条。
(6)使用一种基于离散小波分解的高频提取技术，在乳房区域内寻找微钙化可能存在的位置。
微钙化通常表现为乳房区域内的局部高对比度小块，在小波分解的第二、三层的高频合并带中表现为具有较高模值的点或小块。因此，通过在相应分解层的高频合并带中搜索模值大于一个经验阈值的小波系数，就能获得乳房区域内的局部高对比度位置，当作疑似微钙化。
(7)使用与步骤(4)相同的技术，将已知大致部位的疑似微钙化从乳房组织背景中提取出来。
(8)从已分割出来的疑似微钙化中抽取一组灰度与形状特征，然后使用与步骤(5)相同的分类器来鉴定疑似微钙化的真伪。
这里所采用的特征组共有8个特征，涉及到灰度和形状等，它们包括：区域的平均灰度、区域的灰度标准差、边缘带的对比度、区域的面积、区域的紧凑度、边界的凹陷点比率(微钙化边界线上的凹陷点个数与微钙化边界周长之间的比值)、区域的空洞比率(微钙化区域内部的空洞个数与微钙化区域面积之间的比值)和区域的开运算腐蚀率(将二值化的微钙化区域进行一定深度的开运算后的区域面积与原始面积之间的比值)。
(9)使用经典的基于二维空间欧氏距离的合并聚类(即最小生成树，Minimum Spanning Tree，MST)方法，将原本孤立的微钙化聚合起来，形成微钙化簇。
(10)使用一种基于双视角交线坐标估测的弹性匹配技术，确定两个视角中同一病灶(包括肿块与微钙化簇)的图像元素之间的匹配关系。
乳腺钼靶X线摄影的两个视角(MLO和CC)之间接近正交，只有在这两个平面交线(Z轴)方向的坐标才存在关联性，在对来自两个视角的同一病灶的图像元素进行匹配时，只能依靠Z轴的坐标。因此，首先必须先确定乳腺CR图像中Z轴的位置。在CC视角中，由于拍摄角度是垂直方向的，乳房区域近似于一个水平轴对称的半椭圆形，只要在乳房区域作一条水平中分线即可得到Z轴。而在MLO视角中，拍摄角度是斜向的，且由于重力原因乳房本身会产生一定的下垂，因此乳房区域虽然也接近椭圆形，但其对称轴方向未知。因此，先将MLO视角中的乳房区域边缘点提取出来，然后使用经典的椭圆Hough变换技术来获取该乳房区域的最佳拟合椭圆，这个椭圆的长短轴中指向乳房内侧的轴就是MLO视角中的Z轴。无论在哪个视角中，Z轴的原点都在Z轴与乳房区域边界线的交点位置，这个位置通常就在乳头处或离乳头很近。
对于第一个视角A中的任何一个病灶(包括肿块与微钙化簇)，首先测量其Z轴坐标zi，然后到另一个视角B中，将Z轴坐标处于[zi-Δz，zi+Δz]区间(这里的Δz是一个经验常数)的横截带作为搜索范围，在其中寻找对应的病灶：如果存在一个同类型的病灶(即肿块对应肿块，微钙化簇对应微钙化簇)，就实施匹配；如果存在多个同类型的病灶，则选择Z轴坐标最接近zi的病灶进行匹配；如果不存在同类型的病灶，则选择步骤(5)或步骤(8)中识别为假的同类型病灶中最接近真(即其MLPN输出数值最大)的那个进行匹配，如果假的同类型病灶也不存在，则放弃匹配。由于在病灶匹配时，并非严格遵守Z轴坐标相同的原则，而是放宽到一定的邻域，且当邻域内病灶无法匹配时还允许识别为假的病灶参与匹配，因此该方法被称作弹性匹配。
(11)通过一组灰度、形状与纹理特征和多分类器融合技术，将同一病灶(包括肿块与微钙化簇)在两个视角中的特征信息进行融合，确定病灶的良恶性。
这里使用的特征组用来判别病灶的良恶性，因此与步骤(5)和步骤(8)中的特征组有所区别。用来表征肿块的特征组共有10种特征，涉及到灰度、形状和纹理等，它们包括：区域的平均灰度、区域的灰度标准差、边缘带的对比度、边界的梯度方向熵、边界的平均离心距离、边界的离心距离标准差(步骤(5)中提到的每个标准化角度上的中心射线段长度所组成的序列的标准差)、区域的面积、区域的紧凑度、区域的狭长度和区域的Hausdorff分形维数。用来表征微钙化簇的特征组共有6种特征，主要涉及到微钙化簇的形态，它们包括：簇中微钙化个数、簇面积、微钙化的平均面积、微钙化的平均紧凑度、所有微钙化的离心距离标准差的平均值和所有微钙化的边界梯度标准差的平均值。
对于在步骤(10)中已匹配成功的每一组病灶，分别在两个视角中将其对应的特征测量出来，放入每个视角所对应的MLPN分类器中，然后将两个视角的MLPN分类器的判别结果进行融合，获得最终的结论。本发明所采用的是一种基于多代理的多分类器融合技术：它将两个视角所对应的MLPN分类器的输出结果(一个处于[0，1]区间的浮点数置信度)作为两个代理的初始意见，然后经过两个代理之间的信息交互，利用两个MLPN分类器对样本的置信度及分类结果之间的相关性，自适应地调节两个代理之间的意见，不断迭代，直到输出最优化的信息融合结果。
对于那些在步骤(10)中未匹配成功的病灶，则只在其对应的那个视角中抽取特征，并放入本视角对应的MLPN分类器中进行判别。
(12)将计算机辅助诊断系统的检测结果制成报告，提供给放射科医师参考。