说明书三维超声乳腺全容积图像感兴趣区域的自动提取方法
技术领域
本发明属于图像处理技术领域，具体为三维超声乳腺全容积图像感兴趣区域的自动提取方法。
背景技术
超声成像因具有无创、实时、可重复性强、费用低廉等优点，在临床上有重要应用。与传统的手持式的二维超声成像相比，ABVS具有全新的成像模式，可以标准化自动扫描乳腺，对图像进行数字化处理，避免使用者的个体差异；ABVS可以进行全乳扫描，较常规超声，增加了重建的冠状切面，从而可以提供比二维图像更多的信息，有很好的可重复性。
由于肿瘤的体积相对于整个ABVS图像来说较小，直接进行肿瘤的分割准确率低。因此，通常需要使用者在几百幅横断面图像中手动标注肿瘤的中心位置或者感兴趣的区域，以进行进一步的分析。这种人工标定的方法非常耗时，并且依赖于使用者的经验。
针对这一问题，本发明提出了一种全自动提取ABVS感兴趣区域的方法。该方法不需要使用者提前标记肿瘤，可以自动寻找感兴趣的区域，并且最终得到包含感兴趣区域的最小椭球。将本发明的方法应用到ABVS图像的自动分析系统中，可以提高整个系统的准确性。
发明内容
本发明的目的是提出一种自动提取三维超声乳腺全容积图像中的感兴趣区域的方法。
本发明提出一种三维超声乳腺全容积图像感兴趣区域的自动提取方法，其具体步骤为：
1．首先将DICOM格式的ABVS图像根据三维方向上的像素点间的距离进行图像重建，使之与实际乳腺图像大小相对应；重建以后得到820幅横断面的图像片、750幅矢状面的图像片；根据扫描深度的不同，得到98~294幅冠状面的图像片；在每个切面上，相邻两幅图像之间的距离为0.2 mm，每幅切面图像上的相邻像素点间的距离也为0.2 mm；
2．根据乳腺在冠状面图像上一般为椭圆形的特点，计算重建后ABVS图像的前十幅冠状面图像片的最小值映射图像，对其阈值处理后，用霍夫椭圆变换获取乳腺在冠状面上的模板，并将其应用于所有冠状面图像片；利用ABVS图像冠状面和横断面三维坐标间的关系，将冠状面上的乳腺位置投影至横断面，确定横断面内乳腺位置，初步缩小横断面感兴趣区域的搜索范围，去除乳腺外噪声、伪影等干扰的影响；
3. 对横断面的820幅图像，把每十幅图像通过最小值映射合并为一幅新的图像，合并后得到82幅图像；对其中每幅图像采用基于最大方向相位方法，结合肿瘤位置特性和纹理特点，提取感兴趣候选区域，得到相应区域的二值图像；
4. 由步骤3，得到一系列连续横断面感兴趣候选区域的二值图像；根据肿瘤的连续性，若某一片的连通区域与前后相邻片无重叠，则称为无关片，先加以去除；对剩下的片，根据其连续性及连通区域中心位置进行分组，每组表示为一个可疑肿瘤；
5. 对经步骤4处理的各组横断面图像片进行分类，提取其二值图像和对应灰度图像的形状、纹理等特征，输入至逻辑回归分类器，得到每组可能为肿瘤的概率，概率最大者为真实的肿瘤区域所对应的横断面图像片；
6．在步骤5确定的肿瘤区域连续片的二值图像中，找到最大的连通区域，确定包含该连通区域的最小椭圆；根据矢状面和横断面图像间的坐标对应关系，找到椭圆中心横坐标对应的矢状面图像片，采用步骤3中提取感兴趣候选区域的方法得到矢状面图像片的感兴趣候选区域，进而确定矢状面上包含感兴趣候选区域的最小椭圆；利用矢状面和横断面的两个椭圆，确定ABVS图像包含感兴趣区域的最小椭球的三个主轴长，从而得到肿瘤的感兴趣区域。
下面就本发明方法的各个步骤涉及的相关技术细节作进一步的具体描述。
关于步骤1. ABVS图像是由乳腺自动全容积成像仪器沿乳腺横断面方向扫描得到的连续图像片重建得到的，图像格式一般为DICOM格式，有三个正交的平面：横断面、矢状面和冠状面；横断面和矢状面上的超声图像与传统的手持式超声仪器获取的二维超声图像比较相近，冠状面上的超声图像可以观察到乳腺的大致轮廓；DICOM文件中包含图像的全部信息，包括图像像素点的灰度值以及图像的像素点间距。因此，可以通过DICOM文件信息的读取，重建三维的ABVS图像。ABVS图像的扫描深度一般在20~60 mm，根据扫描深度的不同，重建后的ABVS图像的大小为(98-294)×750×820。原始的ABVS图像的三个切面示意图如图1所示，重建后横断面、矢状面和冠状面的图像如图2所示。
根据从DICOM文件中读取到的ABVS图像像素点之间的距离与实际的乳腺图像的关系，重建ABVS图像，其相邻像素点间距为0.2 mm。
关于步骤2. 该步骤是为了获取乳腺的模板，排除乳腺周围因接触不良等原因产生的伪影、噪声等干扰。其实施方法为：首先取表示乳腺皮肤表层的冠状面前10幅图像片，计算前10幅冠状面图像片的最小值映射图像，然后将这十幅图像按最小值映射的原理进行合并。对合并后的图像采用OSTU的阈值处理[1]方法得到二值图像；采用形态学的膨胀腐蚀方法得到乳腺的大致轮廓，然后用霍夫椭圆变换的方法检测找到乳腺对应的椭圆的位置，生成一幅二值模板图像，椭圆内部的值为1，椭圆外部的值为0；对于ABVS图像来说，其冠状面图像的纵坐标对应横断面图像的横坐标，因此，可以利用冠状面图像生成的模板图，确定横断面图像上乳腺区域，去除乳腺外的噪声、伪影干扰的影响。
为了降低霍夫椭圆变换的运算量，根据先验知识，对椭圆的长短轴、方向和长短轴之比进行限定。这里限定椭圆的主轴长度在400~800之间，长短轴之比大于0.75，椭圆长轴与x轴正方向的夹角在(0,π)之间。椭圆的方程为：
(1)
其中a表示椭圆长轴，b表示椭圆的短轴，θ表示椭圆长轴与x轴正方向的夹角，(x0,y0)为椭圆的中心。
检测到冠状面上包含乳腺的椭圆以后，利用ABVS图像冠状面和横断面三维坐标间的关系，将冠状面上的乳腺位置投影至横断面，确定横断面内乳腺位置，初步缩小横断面感兴趣候选区域的搜索范围，去除乳腺外噪声、伪影等干扰的影响。该部分模板提取的结果如图3所示。利用冠状面的模板作用于整幅ABVS图像，得到的三个切面的图像如图4所示。
关于步骤3. 该步骤的目的是提取横断面图像片的感兴趣区域。首先，将横断面的连续图像每十幅进行最小值映射合并，减少运算量。然后，对合并后的每幅横断面图像进行感兴趣候选区域的提取。
二维乳腺超声肿瘤图像感兴趣候选区域的提取过程为：
1) 图像预处理。
a) 超声图像的斑点噪声比较严重，所以先进行各向异性斑点噪声消除(SRAD)。对于图5(a)所示的图像，滤波后的结果如图5(b)所示。
b) 减小灰度值的范围[2]，利用线性归一化公式(2)对图像中像素点的灰度值进行调整，结果如图5(c)所示。
(2)
其中Ln为灰度级数，lbound和ubound分别取为Q(0.5)和Q(0.95)，Q是直方图累积分布的分位点函数。
c) 增强低回声区域[3]，采用公式(3)表示的自适应Z型方程进行灰度变换：
(3)
其中za，zb，zc决定了Z型函数的形状。za和zc决定了曲线非线性变换的范围，za一般设置为20，zc为图像的均值，而zb则确定了曲线的倾斜程度，根据灰度分布的斜度来得到：
(4)
(5)
设经过预处理以后得到的图像为，如图5(d)所示。
2) 求取的最大方向相位图PMO。根据文献[4]中所提到的方法，在图像的频域沿着六个方向(0°, 30°, 60°, 90°, 120°, 150°)对图像进行滤波，最终提取具有最大能量方向的相位信息。这六个方向覆盖了整个频谱(0~360°)。
对图像I，为了计算PMO矩阵，首先需要二维的log-gabor滤波器对图像滤波，这里二维log-gabor滤波器的传递函数是对数尺度下的高斯函数，将其分解为两个部分：径向滤波器和角度滤波器。对于每一个方向θ0，构建的滤波器为：
(6)
其中大括号中前一部分表示径向滤波器，后一部分为角度滤波器。ω为径向坐标；θ为角度坐标；ω0为滤波器中心频率，设置四个不同尺度的中心频率，其值分别为1/3、1/(3×1.7)、1/(3×1.72)、1/(3×1.7-3)；σθ决定以θ为中心的角度带宽，将其设置为30°；决定径向带宽，通常将其设置为0.55，以保证滤波效果的同时避免混叠。
图像在经过不同中心频率、不同角度的总共24个log-gabor滤波器滤波后，经快速傅里叶变换(FFT)的逆变换到24个图像矩阵(x为方向，s表示中心频率尺度)，再计算每个角度上的相位矩阵LPA：
(7)
其中表示在角度θ时的相位矩阵；n表示不同中心频率滤波器的数量；为尺度s下的增强相位矩阵，其计算式为：
(8)
(9)
其中imag(x,s)为的虚部，real(x,s)为的实部。
经过六个角度滤波后，得到六个LPA特征矩阵。然后，将这六个矩阵合并为一个特征矩阵。这里，提取某一点上具有最大能量角度的相位特征值作为该点最终的相位特征值，从而获得一个新的相位特征矩阵。因为局部的能量刻画了结构信息，所以具有最大局部能量的角度最接近于边界方向。于是，相位特征PMO在每个像素点上的定义为：
(10)
其中表示在角度θ上的能量矩阵，其计算式为：
(11)
这里imag(x,s)和real(x,s)分别表示图像经过log-gabor滤波器后的虚部和实部。ρ是按式(5)求出的具有最大局部能量的角度，表示在ρ角度上四个尺度中心频率的相位矩阵的叠加。
经过上面的过程，可以得到图像I的最大方向相位图PMO。对如图6(a)所示的预处理后的图像，其最大方向相位图如图6(b)所示。然后，为了突出肿瘤区域，去除背景噪声，将PMO与256-相乘，对相乘后的图像采用5×5的中值滤波器进行滤波, 结果如图6(c)所示。但是，经过最大方向相位提取以后，PMO的大部分像素灰度值范围为[0,0.5]，因此整幅图像看起来比较模糊灰暗。采用下面两个公式进行灰度校正：
(12)
(13)
综上，得到灰度图像的最大方向相位图。为了使其中灰度值大的区域表示肿瘤区域，对最大方向相位图取反，最后所得图像如图6(d)所示。
3) 由于肿瘤类型多种多样，其超声图像的纹理和形状差异也比较大，有些肿瘤内部灰度不均匀，但是肿瘤的边缘有连续的相位；有些肿瘤内部灰度较为均匀，但是边缘的相位不明显或者与背景中的相位混叠。因此，在感兴趣候选区域的选取中，需要结合灰度图像和最大方向相位图的信息。用OSTU阈值法对预处理后的图像I、归一化的PMO和PMO+I三幅图像分别进行阈值分割处理，得到含有很多个连通区域的二值图像。去除图像中面积小于300个像素点和与边缘的连通率大于50%的噪声区域。计算剩下所有的连通区域的能量函数E：
(14)
其中Compact表示连通区域的紧致度，其值是连通区域的面积与其最小凸边形的面积之比；Area为连通区域的面积；Centdis为连通区域中心与图像中心的距离；Eccentricity为连通区域最小外接椭圆的离心率；为区域与边缘重合率，其值为区域与边缘重合的像素点的数量占图像边缘像素点的总数量之比。w1~w4分别表示区域的紧致度、面积、中心距和离心率对区域能量函数值的影响大小，一般取其为1，也可以根据具体的超声肿瘤图像数据库进行一定的调整，调整后范围在1±0.5。w为边缘重合率的权重，其取值范围一般为(0,0.5)。
经过上述过程，分别得到预处理后的图像I、归一化的PMO和PMO+I三幅图像的感兴趣候选区域。计算三幅图像感兴趣候选区域的能量函数E，选择E最大的区域为图像最终的感兴趣候选区域。图像I、归一化的PMO和PMO+I三幅图像及其阈值分割和区域选择的结果如图7所示。
关于步骤4. 经过上一步骤，得到一系列连续横断面的感兴趣候选区域的二值图像。肿瘤在横断面的图像片是连续的，其对应的二值图像中连通区域有重叠，如图8所示。在图8中，相邻两幅图像片的感兴趣候选区域是彼此重叠的，且感兴趣候选区域的面积先由小变大，再由大变小。根据肿瘤的连续性可以除去二值图像连通区域与前后邻接片无重叠的无关图像片。判断某一图像片是否为无关图像片的依据为：
(15)
其中BWk表示第k幅横断面图像的二值图像，m和n分别为图像轴向和横向的像素点数。当R(k)<50时，认为第k幅图像的感兴趣候选区域与相邻片无重叠，在图像集合中去除第k幅图像。作为例子，图9给出了 ABVS图像的连续三幅横断面图像及其所对应的感兴趣候选区域的二值图像，其中第二幅图像的感兴趣候选区域与前后两幅图像皆无重叠，因此将其从该ABVS的图像集中去除。
对剩下的图像片，先根据肿瘤区域在横断面图像片中的连续性，将连续的图像片分为一组；然后在每组连续的图像片中，计算相邻片的感兴趣候选区域中心之间的距离，将距离在50个像素点以内的图像片归为一组。这样，最终划分得到的每组表示一个可疑肿瘤的连续横断面切片。在这些组中除了真正的肿瘤区域外，还包含了乳头后方阴影和一些因探头接触不良等原因导致的低回声区域。
关于步骤5. 该步骤的目的是提取每组连续二维图像对应的二值图像和灰度图像的形状和纹理特征，并利用这些特征进行分类。
通过计算可以得到每幅图像片的形状和纹理特征，但是因为每组中包含数量不等的图像片，每组图像片获取的每个特征都有与图像片数量相对应的多个值，因此需要对其进行归纳，得到一个可以用来描述每组图像整体特征的特征值。
1) 形状特征：
组内横断面图像片中的感兴趣候选区域的纵横比：一般来说肿瘤的纵横比小于1，而一些非肿瘤区比如乳腺后方阴影区域的纵横比要大于1，因此可以将感兴趣区域的纵横比作为区分肿瘤和非肿瘤的一个特征。该参数通常用肿瘤的轴向最长径Lshort和横向最长径Llong的比值来表示：
(16)
计算组内每幅图像感兴趣区域的纵横比Mabr，取其中值作为这一组的纵横比。
连续图像片间感兴趣区域的重叠率：由图8可以看到，肿瘤的连续横断面图像片的感兴趣区域的面积由小变大，再由大变小，一般来说，位于中间位置的图像片具有面积最大的感兴趣区域。
定义相邻两幅图像的感兴趣区域重叠率为：
(17)
其中BWk表示组内第k幅图像的二值图像。式(17)将重叠率定义为第k幅图像和第k+1幅图像感兴趣候选区域的重叠面积与第k+1幅图像感兴趣候选区域的面积之比。这样，对于一个含有l幅图像的可疑肿瘤图像组，我们可以得到l-1个重叠率。为了对这l-1个特征进行归纳，分别求其均值、标准差、乘积和梯度均值，作为分类的特征。
2) 纹理特征：
纹理特征包括灰度共生矩阵[5]特征和灰度特征。
灰度共生矩阵[5]特征包括熵、对比度、相关度、能量和同质性。为了计算的简便，先求每幅二值图像中感兴趣候选区域的最小外接矩形，然后将最小外接矩形分别沿纵向和横向向外扩展20个像素，得到包含感兴趣候选区域的矩形。计算位于该矩形内的灰度图像的纹理特征，然后计算组内所有图像片的灰度共生矩阵特征的均值，作为用于分类的特征。
灰度特征包括：a) 感兴趣候选区域的灰度均值和标准差；b) 感兴趣候选区域与包含感兴趣候选区域的矩形内背景区域的灰度均值之比和标准差之比，一般来说，感兴趣区域内部的灰度均值低于背景区域灰度均值，灰度标准差也低于周围背景区域灰度标准差；c) 感兴趣区域前后方的声影特征MeAP，用以区分乳头后方阴影等低回声区域，感兴趣区域前后方的声影特征MeAP的计算式为：
(18)
其中有：
(19)
(20)
(21)
这里m和n分别为图像轴向和横向的像素点数，Ih和BWh分别表示82幅ABVS横断面图像片中的第h幅图像和它的感兴趣候选区域的二值图像，MBW表示面积最大的感兴趣候选区域，MeA为MBW中感兴趣候选区域对应的可疑肿瘤前方灰度均值，MeP为MBW中感兴趣候选区域对应的可疑肿瘤后方灰度均值，Me为可疑肿瘤所有感兴趣候选区域的灰度均值。
综上，最终提取的特征有五个形状特征，十个纹理特征。将这些特征输入至逻辑回归分类器进行分类，得到每组可疑肿瘤区域可能为肿瘤的概率，选择概率最大者为肿瘤。为了验证方法的有效性，采用十次交叉验证的方法。
关于步骤6. 该步骤的目的是通过获取的肿瘤区域的连续片的二值图像和灰度图像，确定包含感兴趣区域的椭球的三个主轴长，从而得到椭球形状的感兴趣区域。
首先，找到位于肿瘤区域中间位置的面积最大的感兴趣候选区域，用形态学的膨胀法将感兴趣候选区域向外扩展10个像素点，计算扩展后的区域的最小外接椭圆，椭圆的长轴长为a，短轴长为b。由图1可知，横断面上的横轴坐标对应矢状面图像片的序号，设感兴趣区域外接椭圆的中心点坐标为(x0，y0)，则取第x0片矢状面图像片，在该矢状面图像的y轴坐标在(y0-b/2，y0+b/2)的区域内利用上述基于最大方向相位的感兴趣候选区域提取方法获取矢状面感兴趣候选区域。计算该感兴趣候选区域向外扩展10个像素后区域的外接椭圆及其长轴长c。这样，就得到了包含感兴趣区域的椭球的三个主轴长a、b和c。
椭球的方程为：
(22)
其中(x0,y0,z0)为椭球中心点的坐标，a、b、c分别为椭球在三个正交方向上的轴长。ABVS图像及其椭球形感兴趣区域示意图如图10所示，实际ABVS图像最终得到的感兴趣区域如图11所示。
附图说明
图1. ABVS图像横断面、矢状面和冠状面示意图，其中Horizontal plane表示横断面示意图，Sagittal plane表示矢状面示意图，Coronal plane表示冠状面示意图。
图2. 重建后ABVS图像在横断面、矢状面和冠状面的视图。方形线框内包含乳腺肿瘤。
图3. 冠状面模板获取，椭圆形线框内为提取的模板。
图4. 应用椭圆模板限制感兴趣区域搜索范围：(a)原始图像的三个切面，线框内为采用椭圆形模板后各切面的图；(b)应用椭圆模板后的三个切面。
图5. 横断面图像预处理。其中，(a)原始图像；(b)SRAD滤波后图像；(c)灰度范围减小的图像；(d)低回声区域增强图像。
图6. 横断面图像最大方向相位图。其中，(a)原始图像；(b) PMO图像；(c) PMO×(256-I)中值滤波后的图像；(d)灰度调整后图像。
图7. 感兴趣区域的提取。其中，(a)、(b)、(c)分别为I、PMO、I+PMO；(d)、(e)、(f)分别为三幅图像的OSTU二值图像；(g)、(h)、(i)分别为三幅图像得到的感兴趣候选区域；矩形线框内为最终感兴趣候选区域。
图8. 肿瘤区域连续图像及其对应的感兴趣候选区域。左侧为连续的原始图像片，右侧为其对应的二值图像，白色区域部分为感兴趣候选区域。
图9. 非肿瘤区域连续图像的感兴趣候选区域。左侧为连续的原始图像片，右侧为其对应的二值图像，白色区域为感兴趣候选区域。
图10. ABVS图像及其椭球形感兴趣区域示意图。其中，(a)横断面感兴趣区域图；(b)矢状面感兴趣区域图；(c)冠状面感兴趣区域图。
图11. 肿瘤感兴趣区域提取结果：(a)肿瘤最大感兴趣区域三个切面的图像；(b)提取出的感兴趣区域三个切面的图像。
具体实施方式
对本发明提出的三维超声乳腺全容积成像(ABVS)中感兴趣区域的自动提取方法进行测试。ABVS图像取自西门子公司的ACUSON S2000TM超声仪器。该系统装配有宽带线性探头(14L5BV)，可以获得15.4 cm×16.8 cm×(2~6) cm的乳腺容积图像。本实验中共采集了15例ABVS图像，每一例有98~294幅冠状面，820幅横断面图像，750幅矢状面图像。
首先，对原始的ABVS图像进行重建，重建后图像的三个切面(横断面、矢状面、冠状面)如图2所示，其中在冠状面图像上可以看到乳房的大致轮廓。由图2可知，乳房轮廓一般接近于椭圆形，因此在冠状面的图像上用霍夫变换找到表示乳房的椭圆，如图3所示。在接下来的处理中，仅需考虑处理椭圆内的区域，可以减小运算量，排除乳房外噪声、伪影等的干扰。图4所示为原始各切面原始图像与应用椭圆模板以后各切面的图像。
图5所示的是对横断面图像进行预处理的过程。图5(a)所示为原始图像。针对超声图像的对比度低、斑点噪声严重的特点，首先采用各向异性斑点噪声消除，其结果如图5(b)所示。然后，减小灰度范围，实际上是对位于0~255中间的灰度值进行一定的拉伸，其结果如图5(c)所示。最后，增强低回声区域灰度值，其结果如图5(d)所示。
图6所示为横断面图像的最大方向相位图。由于有些肿瘤的内部灰度值差异较大，如果只利用灰度图进行阈值分割，无法得到完整的肿瘤区域，因此，增加相位信息。图7所示为I、PMO和I+PMO三幅图像感兴趣区域的提取以及最后从三者中选取的感兴趣候选区域。这里，最后选取的感兴趣候选区域是通过对图像I+PMO进行阈值分割和区域选择得到的。
图8和图9分别是非肿瘤区域和肿瘤区域横断面连续图像及其对应的感兴趣候选区域。由此可以看出，肿瘤区域图像的感兴趣区域是连续的，彼此重叠的；而非肿瘤区域图像的感兴趣区域是彼此不重叠的。根据这个特点可以去除掉一部分的无关图像片。
图11为肿瘤感兴趣区域提取结果。对去除无关片的图像进行分组，然后计算每组图像片的纹理和形状特征，并且利用逻辑回归分类器进行分类，采用十次交叉验证的方法验证其有效性。
在实验中，共对15例ABVS图像进行处理，其中14例都可以得到含有肿瘤区域的组，准确率达到93.3%。利用含有肿瘤区域的组中的图像片的灰度图像和二值图像，结合矢状面的图像，计算包含肿瘤区域的最小椭球的三个正交方向的主轴长，得到肿瘤的椭球形感兴趣区域，其示意图如图10所示。图11(a)为ABVS图像的三个切面，最终得到的实际感兴趣区域的三个切面如图11(b)所示。采用椭球表示感兴趣的区域，得到的椭球比立方体区域更接近真实的肿瘤区域。
综上，本发明适合用于三维乳腺全容积图像中肿瘤感兴趣区域的提取，整个提取的过程是完全自动的、准确的，可以实现肿瘤区域的定位，不再依赖使用者的标注，更加客观和高效。
参考文献
[1] N. Otsu, “A threshold selection method from gray level histograms,” IEEE Trans. Syst., Man Cybern., vol. 9, pp. 62–66, 1979.
[2] B. Liu, H.D. Cheng, J.H. Huang, J.W. Tian, X.L. Tang, J.F. Liu, “Fully automatic and segmentation-robust classification of breast tumors based on local texture analysis of ultrasound images,” Pattern Recognit., vol.43, pp. 280–298, 2010.
[3] M. Xian, Y.T. Zhang, H.D. Cheng, “Fully automatic segmentation of breast ultrasound images based on breast characteristics in space and frequency domains,” Pattern Recognit., vol.48, pp. 485–497, 2015.
[4] J. Shan, H.D. Cheng, Y. Wang, “Completely automated segmentation approach for breast ultrasound images using multiple-domain feature,” Ultrasound in Med. & Biol., vol.38, pp. 262–275, 2012.
[5] R. M. Haralick, K. Shanmugam, I. H. Dinstein, “Textural features for image classification,” IEEE Trans. Syst., Man Cybern., vol.3, pp. 10–621, 1973。