一种超声图像中穿刺活检针的检测方法

摘要

本发明公开了一种超声图像中穿刺活检针的检测方法，包括以下步骤：根据背景差分法从穿刺活检针插入后的超声图像中找出感兴趣区域；接着在感兴趣区域内选取一个窗口S，旋转该窗口S得到灰度投影矩阵GPM，并根据矩阵极值所在的行和列得出直线方程的参数；然后根据该直线方程的参数新选取一个窗口S0，旋转得到新的灰度投影矩阵NPGM，并根据矩阵极值所在的行和列得出再次直线方程参数，从而得出超声图像中穿刺活检针的位置。本发明通过对该检测方法整体流程以及各个步骤所基于的原理、相应处理所采用的具体计算过程等进行改进，能够提高超声图像(尤其是前列腺超声图像)中穿刺活检针的检测正确率，同时减少处理时间以满足实际临床的实时需求。

说明书

技术领域

本发明属于图像分析与处理中的目标检测领域，更具体地，涉及一种超声图像中穿刺活检针的检测方法，该检测方法是基于动态窗口和粗-精搜索的灰度投影积分算法，尤其适用于前列腺超声图像中活检针检测，是种快速精确的检测方法。

背景技术

目前，介入式诊断与治疗在现代医学影像中发挥着越来越重要的作用。而在介入式手术中如何准确快速的定位出手术器械的位置一直是一个热点和难点问题。在现代医学影像技术中，超声成像技术以其独有的非损伤性和经济性在临床手术中发挥了越来越重要的作用，已经成为医院诊断和治疗的常规手段，尤其在前列腺穿刺活检手术中，通常会作为医生观察病灶的首选。

在实际手术中，由于成像部位周围组织的影响会干扰到医生对活检针具体位置的判断，有可能会做出错误的诊断结果。精确确定穿刺活检针的检测方法可帮助医生方便地找出活检针的位置，辅助医生进行诊断和治疗，从而提高穿刺以及活检手术的准确率，减轻病人的痛苦。

尤其对应前列腺穿刺活检术，在前列腺穿刺活检术中，前列腺超声图像中的活检针具有以下特性：活检针一般由金属材料制成，在超声成像中属于高回声物体，因此在超声图像中会呈现出高亮、直线状这一特性。基于这个特性，国内外学者提出了以下几种活检针检测算法。第一种是Hough变换法，该方法是检测直线的最经典、有效的方法之一，常用的有标准Hough变换、随机Hough变换，这些方法将图像空间中较为困难的全局检测问题转化为参数空间中相对容易解决的局部峰值检测问题，可以识别和检测直线。第二种是基于投影的方法，例如积分投影算法，Radon变换法。第三种是基于迭代的方法，例如随机噪声一致性算法。第四种是基于统计学的方法，例如主成分分析法。上述方法中，基于Hough变换的方法存在计算复杂度高和运算量大的缺点，后两种算法的稳定性不足，容易得到错误检测结果且难以达到临床实时应用的要求。

发明内容

针对现有技术的以上缺陷或改进需求，本发明的目的在于提供一种超声图像中穿刺活检针的检测方法，其中通过对该检测方法整体流程以及各个步骤所基于的原理、相应处理所采用的具体计算过程等进行改进，与现有技术相比，由于该检测方法基于动态窗口、以及粗-精搜索的灰度投影积分，能够提高超声图像(尤其是前列腺超声图像)中穿刺活检针的检测正确率，同时减少处理时间以满足实际临床的实时需求。

为实现上述目的，按照本发明，提供了一种超声图像中穿刺活检针的检测方法，其特征在于，包括以下步骤：

(1)根据背景差分法利用穿刺活检针插入前后的超声图像从穿刺活检针插入后的超声图像中找出灰度变化满足预先设定要求的区域，该区域包含已插入的穿刺活检针的超声成像，记该区域为感兴趣区域；

所述感兴趣区域为矩形区域，在所述感兴趣区域上建立图像坐标系{F:Ouv}，其原点坐标O是该感兴趣区域的起始位置，水平向右为u轴，垂直向下为v轴；

(2)在所述感兴趣区域内选取一个矩形形状的窗口S，记该窗口S的中心为O_S，该在窗口S上重新建立一个新的坐标系{F_k:O_Su_kv_k},所述中心O_S即为该坐标系{F_k:O_Su_kv_k}的坐标原点，并且记该坐标原点O_S在{F:Ouv}坐标系下的坐标为O_S(u_S,v_S)；该坐标系{F_k:O_Su_kv_k}中，记坐标轴u_k为投影方向，坐标轴v_k为扫描方向，与坐标轴u_k垂直的平面Q记为投影平面；

将该窗口S及坐标系{F_k:O_Su_kv_k}一起以该窗口S的中心O_S为中心点按照[0°-180°]的范围每隔一定的角度θ旋转，每次旋转后将该窗口S内的超声图像在平行于所述投影方向上的像素灰度值进行累加得到一个灰度投影矩阵GPM，直到旋转全部完成后根据矩阵极值所在的行和列得出直线方程的参数；所述矩阵极大值所在的行即对应直线的斜率，所述矩阵极值所在的列即对应直线距所述窗口S中心O_S的距离；

所述窗口S的宽度为2b+1,高度为2d+1，b和d的取值为在保证该窗口S在旋转过程中不超出所述感兴趣区域情况下的任意整数，单位为像素；

(3)根据所述步骤(2)得到的所述直线方程参数，在该直线附近重新选取一个矩形形状的窗口S0，记该窗口S0的中心为O_S0，将该窗口S0以O_S0为中心点沿着所述直线的方向角度按照预先设置的某一范围每隔一定的角度旋转，相应得到一个新的灰度投影矩阵NPGM，根据该NPGM矩阵极值所在的行和列，再次计算出直线方程，得到再次直线方程参数；

(4)根据所述步骤(3)得到的所述再次直线方程参数，得出超声图像中穿刺活检针的位置，该穿刺活检针即位于与所述再次直线方程参数相对应的直线上。

作为本发明的进一步优选，记所述步骤(1)中所述感兴趣区域的灰度分布函数为f(u,v)；

则所述步骤(2)中，当第k次旋转完成后，所述灰度投影矩阵GPM满足：

GPM＝[p₀>1>2 ...>k]^T

其中，p_k＝[g_-d,k>-d+1,k>-d+2,k ...>d,k]；

若记该窗口S内的一个像素点在所述坐标系{F_k:O_Su_kv_k}下的坐标为(u_i^k,v_i^k)，该像素点在所述坐标系{F:Ouv}下的坐标为(u_i,v_i)，则两者的转换关系满足：

并且，该窗口S内经过该像素点、且沿平行于所述投影方向上的超声图像像素灰度值累加值满足：

v^k分别取-d,-d+1,-d+2,…,d；

直到旋转全部完成后根据矩阵极值所在的行和列得出直线方程的参数；

记所述矩阵极值位于矩阵的第m行，第n列，则：

其中，k'为直线的倾斜角；h为当所述窗口S的旋转角度为k'时，原点O_s距离该直线的距离。

作为本发明的进一步优选，所述步骤(3)中，所述窗口S0是根据所述步骤(2)得到的直线方程参数k'、h得出的，具体是将旋转角度为k'的所述窗口S及其坐标系{F_k:O_Su_kv_k}沿着垂直于该直线的方向向直线靠近h个像素平移所得，平移后窗口S的坐标原点O_S即对应窗口S0的坐标原点O_S0，然后以O_S0为该窗口S0的中心点重新选取窗口S0的宽度和高度，记窗口S0的宽度为w_S0＝2b_S0+1,高度为h_S0＝2d_S0+1，记所述坐标原点O_S0在{F:Ouv}坐标系下的坐标为O_S0(u_S0,v_S0)，则：

该窗口S0的旋转角度的取值范围为[k'-△k'+△]，其中△为预先设定的角度，b_S0和d_S0的取值为在保证该窗口S0在旋转过程中不超出所述感兴趣区域情况下的任意整数，单位为像素；

将该窗口S0以该窗口S0的中心O_S0为中心点从(k'-△)开始每隔一定的角度θ'旋转，相应得到所述新的灰度投影矩阵NPGM，记NPGM矩阵极值位于矩阵的第m₀行，第n₀列，则：

若记直线方程为y＝Kx+c，该直线上一点为P(i,j)，

如果k_s≠90°，则K＝k_s

c＝v_S0+h_ssin(k_s)-(u_S0-h_scos(k_s))×tan(k_s)

将K、c代入y＝Kx+c即得到直线方程；

如果k_s＝90°，则直线方程为：x＝u_S0；

该直线方程即为再次直线方程。

作为本发明的进一步优选，所述步骤(4)具体包括：

记与所述再次直线方程参数相对应的直线为目标直线，将所述感兴趣区域内被该目标直线经过的各个点的灰度值由其八邻域均值代替，并将所述各个点的灰度值它们的均值作为阈值T，然后对所述各个点进行判定，记该点的八邻域灰度均值为g₀，则：

若|g₀-T|小于预先设置的θ”，则该点被所述穿刺活检针经过；

若|g₀-T|不小于θ”，并且该点距离被所述穿刺活检针经过的其他点的最短距离不超过预先设定的长度值，则该点依然被所述穿刺活检针经过；

若|g₀-T|不小于θ”，并且该点距离被所述穿刺活检针经过的其他点的最短距离超过所述预先设定的长度值，则该点没有被所述穿刺活检针经过；

由此判断所述超声图像中穿刺活检针的针尖位置。

作为本发明的进一步优选，所述步骤(2)中，所述矩阵极值为所述矩阵GPM的最大值；

所述步骤(3)中，所述NPGM矩阵极值为所述矩阵NPGM的最大值。

作为本发明的进一步优选，所述步骤(3)中所述窗口S0的宽度与所述步骤(2)中所述窗口S的宽度保持相等，所述窗口S0的高度也与所述窗口S的高度保持相等。

总体而言，本发明提出的上述技术方案与现有技术相比，能够取得下列有益效果：

本发明超声图像中穿刺活检针的检测方法，是同时基于动态窗口和粗-精搜索的灰度投影积分检测穿刺活检针，通过对该检测方法整体流程以及各个步骤所基于的原理、相应处理所采用的具体计算过程等进行改进，利用各个步骤的整体配合，能够提高超声图像(尤其是前列腺超声图像)中穿刺活检针的检测正确率，同时减少处理时间以满足实际临床的实时需求。

本发明中超声图像中穿刺活检针的检测方法，是先根据背景差分法确定穿刺活检针的插入区域，即图像的感兴趣区域；然后，在感兴趣区域内选取一个合适的窗口，按照一定的角度将像素灰度值进行累加得到一个灰度投影矩阵，根据矩阵极值所在的行和列得出直线方程参数；其次，根据上一步得到的参数在直线附近重新选取一个窗口，重复上一个步骤再次计算出直线方程参数；最后，根据图像块之间的相关性得到穿刺活检针的针尖位置。

对于穿刺活检成像技术，以前列腺穿刺为例，前列腺穿刺手术成功与否的关键在于能否将活检穿刺针准确地放置到病人的病灶区域。而超声导引下的影像具有以下特点：第一，活检穿刺针在影像中会以一条高亮的细线出现，但是由于探头的移动会导致活检穿刺针的某些部分没有出现在影像中从而形成断续的影像；第二，背景图像中有可能出现类似活检穿刺针的细长、高亮的物体(比如组织边缘)，会对活检穿刺针的检测起到干扰作用；第三，由于影像中活检穿刺针前端端点与周围组织混在一起，这就导致前端端点的位置难以准确确定，这是临床医生比较关注的问题，也是判断检测结果准确性的一个重要指标；第四，在实际临床应用中要求算法处理速度足够快，一般要求达到实时应用，这样才具备实际临床应用价值。穿刺针一般是由金属材料制成，在超声图像中属于高回声物体，是以一条高亮的直线形式呈现，因此，前列腺穿刺手术中的活检穿刺针的检测也就是对图像中的直线进行检测，人们通常采用线状结构的分割、检测算法来处理，比较常用的直线检测算法有以下几种：1)Hough变换算法，是直线检测中最经典、最有效的方法之一，将图像空间中较为困难的全局检测问题转化为参数空间中相对容易解决的局部峰值检测问题，具有局部缺损不敏感、对随机噪声具有较好的鲁棒性等特性，但是它具有计算量大，内存消耗大的问题且一般是在二值化图像上进行，使基于Hough的直线检测算法在前列腺穿刺手术中的应用收到了限制；2)相位编组法，与Hough变换主要利用梯度的幅度来提取边缘不同的是：它充分利用了梯度的方向相位信息，根据梯度方向的一致性，形成边缘支持区域，相位编组法将灰度变化的相位作为对边缘特征做局部决策的首要考虑因素，消除了对比度以及相位变化对边缘造成的不清晰的影响，这就使得相位编组法相对其他直线检测算法更容易检测出弱对比度的直线；但是却受到图像噪声的影像，噪声像素点梯度方向相位发生突变，使得相位编组法在检测直线的过程中受直线中的间隙和噪声影响比较大。

本发明中基于灰度投影的穿刺针检测算法，直接利用图像的灰度信息，避免了Hough变换算法检测结果过度依赖于阈值的选取问题，并且降低了算法的复杂度，本发明通过粗-精搜索地方法，大大减少了处理所需的时间，使其满足实际临床的实时需求。本发明通过动态窗口的选取，提高了直线检测的准确率。本发明通过利用图像块之前的相关性，提高了直线检测结果的精度，降低了活检针针尖误差，使其符合临床医生对误差的要求。

综上，本发明中超声图像中穿刺活检针的检测方法是种基于动态窗口和粗-精搜索的灰度投影积分方法，尤其适用于前列腺穿刺活检针的快速精确检测。本方法利用动态窗口和粗-精搜索的方法不仅大大减少了灰度投影积分的处理时间，而且提高了活检针的检测精度。

附图说明

图1是本发明所基于的前列腺活检检测活检针检测算法其流程示意图。

图2是端点位置误差图。

图3是角度误差图。

图4是本发明与未改进的灰度投影方法比较时所使用的前列腺穿刺手术的临床超声序列图以及改进前后方法的检测效果图；其中，(a)、(d)、(g)为一组，(b)、(e)、(h)为一组，(c)、(f)、(i)为一组；(a)、(b)、(c)分别是序列图像中的第3、10、27幅，(d)、(e)、(f)对应改进前的灰度投影(GPI)检测结果图，(g)、(h)、(i)对应改进后的灰度投影(IGPI)检测图。

图5是感兴趣区域内基于旋转窗口S的投影方向及投影平面示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

如图1所示为本发明中超声图像中穿刺活检针的检测方法其流程示意图，该超声图像中穿刺活检针的检测方法是基于动态窗口和粗-精搜索的灰度投影积分进行活检针检测的，该检测方法如图1所示包括以下步骤：

步骤1根据背景差分法找出活检针的插入区域，即感兴趣区域。

例如，可根据超声图像(如前列腺超声图像等)，将活检针插入前后超声图像中灰度变化达到预先设定要求(例如，大于某一预先设置的灰度变化值等)的区域记为感兴趣区域(若灰度变化达到要求的区域为分离的多个区域，则可将连接这些区域的区域也列入感兴趣区域中，从而使最终得到的感兴趣区域覆盖所有灰度变化达到要求的区域)。

步骤2选取一个合适的窗口S，以窗口中心O为中心点按照[0°-180°]的范围每隔一定的角度θ(如每1°)将该方向上的像素灰度值进行累加得到一个灰度投影矩阵(GrayProjection Matrix简称GPM)，根据矩阵极值所在的行和列得出直线方程的参数；

步骤3根据步骤2计算出的直线方程参数在直线附近重新选取一个窗口S0，以新窗口O₀为中心点沿着步骤2得到的直线方向角度按照设置的某一范围每隔一定角度再次进行灰度投影积分，得到一个新的灰度投影矩阵NPGM，再次计算出直线方程。

步骤4根据图像块之间的相关性得到活检针的针尖位置。

上述基于动态窗口和多尺度的灰度投影积分的算法，步骤1中利用背景差分法找出感兴趣区域的计算方法可以为：

假设输入图像为活检针已经插入的图像I，另有一幅活检针未插入的图像I₀，两幅图像对应灰度值相减得到灰度值变化满足要求的区域(对于前列腺超声图像，由于所采用的活检针一般为金属材料，在超声图像中会呈现出高亮的特性，因此活检针的插入区域往往对应于灰度值变化最大的区域)，该区域为所选择的感兴趣区域(ROI)。

上述基于动态窗口和多尺度的灰度投影积分算法，步骤2中计算灰度投影积分矩阵的方法可以为：

如图5所示，假设步骤1中得到的感兴趣区域为如图中所示为F，其灰度分布函数为f(u,v)。在此区域建立图像坐标系{F:Ouv}，其原点坐标O是感兴趣区域的起始位置，假设原点在输入图像I中的坐标为O(m,n)。水平向右为u轴，长度为w；垂直向下为v轴，长度为h，单位为像素。

所选窗口为S，在窗口S上重新建立一个新的坐标{F_k:O_Su_kv_k},坐标原点O_S在{F:Ouv}坐标系下的坐标为O_S(u_S,v_S)。子窗口宽度为w_S＝2b+1,高度为h_S＝2d+1；b,d的取值在保证窗口S在旋转过程中不超过F区域的情况下为任意整数。u_k为投影方向，v_k为扫描方向。与u_k垂直的平面Q称其为投影平面。

k为旋转角度其取值范围为[0°-180°]，在此步骤中，k每隔θ取值(例如θ为1°)。

将子窗口S连同其上坐标系沿中心点O_S旋转到角度k，将子窗口S沿着投影轴方向投影到投影平面Q上得到灰度投影积分矩阵GPM：

GPM＝[p₀>1>2 ...>k]^T

因为图像处理将像素值离散化，因此每个角度下的投影向量p_k由v_k上的2d+1个投影积分值组成：

p_k＝[p_1k>2k>3k ...>2d+1,k]

假设子窗口S内的一个像素点(u_i^k,v_i^k)在感兴趣区域内F的坐标为(u_i,v_i),则两者的转换关系为：

子窗口中该点所在的那一行所对应的

v^k分别取-d,-d+1,-d+2,…,d-1,d(+1,+2……的单位均为像素)；

则可见：

p_k＝[p_1k>2k>3k ...>2d+1,k]＝[g_-d,k>-d+1,k>-d+2,k ...>d,k]；

在灰度投影积分矩阵GPM中搜索矩阵极值(即最大值)所在的行和列(若最大值为多个，则可取其中任意一个)，假设极值在第m行，第n列。即：

k'值对应直线的倾斜角，h值表示在旋转角度为k'的情况下原点O_S到直线的距离。在此步骤中计算出了直线方程参数。

上述基于动态窗口和多尺度的灰度投影积分检测活检针的算法，步骤3中计算直线方程可包括以下子步骤：

步骤3-1根据步骤2所计算出的k'，h值重新选取一个新的子窗口S0，该新窗口是由原来的窗口S沿着垂直于直线的方向向直线靠近h个像素所得，重新选取窗口的宽度和高度以及坐标原点O_S0(此处可保持窗口的宽度和高度不变)。O_S0在{F:Ouv}坐标系下的坐标为O_S0(u_S0,v_S0)：

新的旋转角度k₀的取值范围为[k'-5°，k'+5°](即将△预先设置为5°)，重复步骤2得到新的灰度投影矩阵NGPM(k₀可每隔θ'取值，θ'与θ的值可以相同)，查找出NGPM极值点所在的行列值(m₀,n₀)，得到直线方程的新参数k_s和h_s。

步骤3-2确定直线方程：

假设直线方程为y＝Kx+c，直线上一点为P(i,j)。

如果k_s≠90°，则K＝k_s

带入直线方程可得c＝v_S0+h_ssin(k_s)-(u_S0-h_scos(k_s))×tan(k_s)

如果k_s＝90°，则直线方程为：x＝u_S0(即垂直于x轴)。

上述基于动态窗口和多尺度的灰度投影积分的算法，步骤4中关于寻找针尖端点位置的计算方法可以为：

将在步骤3所确定的直线上的感兴趣区域内的点的灰度值由其八邻域均值代替，求出所有直线上点的灰度值的均值将其作为阈值T，

假设直线角度上的点的八邻域灰度均值为g₀；如果满足|g₀-T|＜θ”则认为该点在活检针上(θ”可预先设定)，如果不满足则暂时标记，寻找下一个点；如果该标记点距离下一个在活检针上的点的距离不超过step则认为该标记点依然是活检针上的点(step即为预先设定的长度)，如果不满足此条件就认为该点不是活检针上的点(即超声图像中的背景点)。直至遍历结束，不再有在活检针上的点就停止寻找。这样就确定了针尖的位置。

实施例1

为了体现本方法的优点，我们会将本发明改进后的灰度投影方法与未改进的方法在检测结果上面进行比较。

检测结果评价标准有四个：

端点位置误差

角度误差▽α＝|α-α^*|

单帧图像处理时间(毫秒)

检测准确率

其中p^*和α^*分别是金标准(医生手动划分)中的针尖端点和活检针角度，p和α分别是算法检测到的针尖端点和活检针角度。

本实例采用实际临床超声图像进行检测结果测试，本实例所用的超声图像来自于华中科技大学同济医学院附属协和医院，表1列出了两种算法所得端点位置误差、角度误差、单帧图像处理时间、检测准确率，结果表明本发明能在保证算法检测准确率的情况下大大减少处理时间，并且提高了检测精度。

表1各方法检测结果

为了更直观的看出本发明相对于改进前方法的优越性，我们提供了在上述实验中检测结果的视觉效果图。从图中我们可以直观的看出IGPI算法在减少处理时间的基础上提高检测结果的精度，使本算法能够满足临床需求，体现了本方法在实际应用中的优势。

本发明中针对超声图像中的穿刺活检针检测方法，是基于动态窗口和粗-精搜索的灰度投影积分，尤其适用于实现超声图像中前列腺活检穿刺针的检测。本发明窗口S、窗口S0旋转过程中，为求取灰度投影矩阵，各个像素点的灰度值可按某一固定的取整规则进行，如整体向上、整体向下、取中间值等。本发明未详细论述的技术细节，如背景差分法等，均可参考相关现有技术。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。