1.本发明涉及医学图像处理的技术领域，尤其涉及一种超声探头自动空间标定方法，以及超声探头自动空间标定装置。


背景技术：

2.徒手跟踪超声成像是一种廉价，安全，无创的技术，跟踪允许将超声波与其他方式进行配准以进行计划和导航。尽管存在3d超声探头，但b型或2d探头由于其相对较高的图像质量和较低的成本而在临床实践中仍被广泛使用。
3.为了进行有效的可视化和导航，使用手持式探头获取的2d us切片必须与患者的基础3d解剖结构相关，包括其他相关信息，例如诊断3d mr或ct扫描。通常，这需要将跟踪器牢固地连接到常规超声探头上，通过外部跟踪系统(例如，光学或电磁跟踪器)连续地监视跟踪器，从而实时提供us探头的位置和方向。这种超声引导的导航系统的关键组成部分是标定2d us图像平面与探头跟踪器之间的刚性转换矩阵。
4.由于标定原理简单，基于n线模型的标定方法仍然是应用最为广泛的方法。在经典的n线体模中，每层n线由两条平行线和一条斜线组成，只有中间斜线的基准点才是有效基准点。对于3个n线的标定体模，每张图像上的有效基准点只有3个，需要采集大量不同位置的图像才得到稳定的标定结果。理论上，可以增加单张图像特征点的数量，从而提高标定过程的效率。pagoulatos首先使用了多层n线标定体模来提高标定效率。但是该方法中的n基准点识别是通过手动选点完成，每张图需要大约2分钟。lindseth提出了使用多层金字塔n线体模实现超声探头标定，证明了多层n线体模标定精度优于单点体模和交叉线体模。hsu设计了一种多层n线体模，其体模需要使用一个薄膜作为地标去辅助分割所有基准点。但是该方法不能够识别位于平行线上的斑点噪声。此外，该方法要求检测到超声图像中的全部n线基准点，缺少任意一个基准点都无法标定。因此，多层n线结构增加了特征的识别难度，尤其是对于图像伪影较多的低频超声探头，因为相似尺度的图像伪影和噪声有时会被误认为是图像特征。此外，在标定过程中还需要考虑斑点识别的准确性。在超声图像中，不同位置斑点的形态结构都不一样，在靠近探头位置斑点成像质量较好，远离探头的位置，图像伪影变大，影响斑点识别的准确性，从而影响标定结果准确性。


技术实现要素：

5.为克服现有技术的缺陷，本发明要解决的技术问题是提供了一种超声探头自动空间标定方法，其准确、高效且易于使用，即使在图像成像缺失和噪声干扰的情况下也可以正确识别，即使没有经验的用户也可以执行精确的探头标定。
6.本发明的技术方案是：这种超声探头自动空间标定方法，其包括以下步骤：
7.(1)设计一种扩展n线标定体模，通过在横向和纵向上排布多个n线结构，增加超声图像中的有效基准点数量；
8.(2)基于扩展n线标定体模超声图像进行快速标定：通过斑点分割、平行线限制和
基准点识别的步骤准确识别所有斑点编号，并通过非线性迭代优化方法计算出超声探头标定矩阵；
9.所述步骤(1)中，扩展n线标定体模的导线孔为三角形，设计多层并且每层上有两个以上的n线结构，一种线为n的斜线，另一种线为n的平行线，将三维n线连接方式投影在前面板上，根据前后面板交点位置定义每条线的编号；在设计连线方式上，上下相邻两层的结构不重复，斜线两端点的水平方向距离小于30mm，相邻层之间的宽度为10mm，共有12层n线，模型中共有51根平行线和33根斜线，采用集合模型中共有51根平行线和33根斜线，采用集合表示体模的l层线集，其中第i层线集中包含：个平行线w
p
、个斜线wo。
10.本发明设计了一种扩展n线标定体模，通过在横向和纵向上排布多个n线结构，增加了超声图像中的有效基准点数量，仅使用几张图像就可以得到收敛的标定结果；设计了基于扩展n线体模超声图像的快速标定方法，该方法通过斑点分割、平行线限制和基准点识别的步骤准确识别所有斑点编号，并通过非线性迭代优化方法计算出超声探头标定矩阵；因此该方法准确、高效且易于使用，即使在图像成像缺失和噪声干扰的情况下也可以正确识别，即使没有经验的用户也可以执行精确的探头标定。
11.还提供了超声探头自动空间标定装置，其包括：
12.扩展n线标定体模，其配置来通过在横向和纵向上排布多个n线结构，增加超声图像中的有效基准点数量；扩展n线标定体模的导线孔为三角形，设计多层并且每层上有两个以上的n线结构，一种线为n的斜线，另一种线为n的平行线，将三维n线连接方式投影在前面板上，根据前后面板交点位置定义每条线的编号；在设计连线方式上，上下相邻两层的结构不重复，斜线两端点的水平方向距离小于30mm，相邻层之间的宽度为10mm，共有12层n线，模型中共有51根平行线和33根斜线，采用集合型中共有51根平行线和33根斜线，采用集合表示体模的l层线集，其中第i层线集中包含：个平行线w
p
、个斜线wo；
13.快速标定模块，其配置来基于扩展n线标定体模超声图像进行快速标定：通过斑点分割、平行线限制和基准点识别的步骤准确识别所有斑点编号，并通过非线性迭代优化方法计算出超声探头标定矩阵。
附图说明
14.图1示出了探头标定过程涉及的组件、连接和坐标系转换关系。标签是指：em(electromagnetic),emt(电磁跟踪器),us(超声图像),es
pr
(固定在超声探头上的传感器),ph(标定体模),es
ph
(固定在体模上的传感器)。类似符号表示从a到b的转换。实线转换表示实时跟踪或者已知的转换矩阵，虚线转换代表着未知的转换矩阵。
15.图2示出了将三维n线连接方式投影在前面板上的图。
16.图3示出了基于扩展n线体模超声图像的标定方法流程图。
17.图4是根据本发明的超声探头自动空间标定方法的流程图。
具体实施方式
18.如图4所示，这种超声探头自动空间标定方法，其包括以下步骤：
19.(1)设计一种扩展n线标定体模，通过在横向和纵向上排布多个n线结构，增加超声图像中的有效基准点数量；
20.(2)基于扩展n线标定体模超声图像进行快速标定：通过斑点分割、平行线限制和基准点识别的步骤准确识别所有斑点编号，并通过非线性迭代优化方法计算出超声探头标定矩阵；
21.所述步骤(1)中，扩展n线标定体模的导线孔为三角形，设计多层并且每层上有两个以上的n线结构，一种线为n的斜线，另一种线为n的平行线，将三维n线连接方式投影在前面板上，根据前后面板交点位置定义每条线的编号；在设计连线方式上，上下相邻两层的结构不重复，斜线两端点的水平方向距离小于30mm，相邻层之间的宽度为10mm，共有12层n线，模型中共有51根平行线和33根斜线，采用集合模型中共有51根平行线和33根斜线，采用集合表示体模的l层线集，其中第i层线集中包含：个平行线w
p
、个斜线wo。
22.本发明设计了一种扩展n线标定体模，通过在横向和纵向上排布多个n线结构，增加了超声图像中的有效基准点数量，仅使用几张图像就可以得到收敛的标定结果；设计了基于扩展n线体模超声图像的快速标定方法，该方法通过斑点分割、平行线限制和基准点识别的步骤准确识别所有斑点编号，并通过非线性迭代优化方法计算出超声探头标定矩阵；因此该方法准确、高效且易于使用，即使在图像成像缺失和噪声干扰的情况下也可以正确识别，即使没有经验的用户也可以执行精确的探头标定。
23.优选地，在所述步骤(1)中，使用直径为0.234mm的尼龙线安装所有n线结构。
24.优选地，在所述步骤(1)中，基于以下三个几何特征，通过形态学操作和像素聚类的方式提取斑点的特征，包括位置、形态和强度
25.a.不同层n线交点的所在直线平行；
26.b.每层n线的交点中，至少有三个点共线；
27.c.平行线和任意位置超声平面的交点集与平行线的二维投影点集之间都存在等比变换。
28.优选地，在所述步骤(2)中，为了全面描述斑点的形态，使用二元高斯函数拟合斑点形状，为公式(2)：
[0029][0030]
其中(xi,yi)是散斑中第i个像素的位置，ii是第i个像素的强度，考虑到散斑的强度分布，该方法中μ和δ分别设置为255和10；之后，对所有聚类像素的强度呈线性总和的散斑进行评分。
[0031]
优选地，在所述步骤(2)的平行线限制中，统计两个候选点连线与水平方向夹角在-20度到20度之间的直线簇，在这些直线簇中找出包含最多直线的角度区间，定义为主方向角度，对角度排序之后，用一个小窗口搜索，找到角度最集中的位置为平行线主方向，选
择角度窗口大小为2度；
[0032]
统计所有保留直线簇的截距，使用具有噪声的基于密度的聚类方法得到截距的聚类集合，保留斑点数量大于3的集合，利用每个集合中的斑点位置和评分计算当前集合的直线方程，并根据直线截距进行排序，从而识别到不同层n基准点所在直线，用d＝{d
ij
,i∈[1,m],j∈[1,ni]}表示图中检测出m层斑点集合，第i层斑点集合包含有ni个斑点d。
[0033]
优选地，在所述步骤(2)的基准点识别中，对斑点集合d构建检索kd树，计算体模坐标系到图像坐标系的初始旋转矩阵
[0034][0035]
其中，α为第一层斑点所在直线方程的角度；s为图像间隔，通过前两层斑点的直线方程和体模配置计算；穷举斑点集合d中前两层斑点和集合w中前两层平行线之间的所有匹配《w
p
,d》；对于每一组匹配，计算平移向量t＝d-r*p，p为平行线w
p
在前面板上的端点，最终得到初始配准矩阵t
2d
＝[r|t]；
[0036]
从线集w第一层开始，使用配准矩阵t计算出每条平行线w
p
在超声图像中投影点p，在kd树中搜索所有可能匹配的斑点d；定义以下匹配相似度
[0037][0038]
其中，v(d)为斑点d的形态评分，‖p-d‖为投影点p和斑点d之间的欧氏距离，s为配准后图像间距；将平行线w
p
与相似度最大的斑点组成映射关系《w
p
,d
p
》，存入平行线映射集合中m
p
；当计算到第三行之后，根据m
p
中的映射关系，重新计算配准矩阵t
2d
，然后重新计算所有平行线的最佳映射关系，存入到m
p
中；完成所有平行线映射之后，计算斜线的映射关系，使用配准矩阵t
2d
计算出每条斜线wo在超声图像中的投影线段l，在kd树中搜索所有位于线段l的斑点{do}m，组成一对多的映射关系《wo,{do}m》，加入候选斜线映射集合在多次循环之后，保留集合m
p
和集合中总映射关系最多的集合；根据n线几何结构，筛选中每条斜线与斑点的唯一映射《wo,do》，记集合m
p
中第i层存在个平行线映射集合集合中第i层任意斜线wo对应m种映射通过组合任意两条平行线和一条斜线wo，构成一组n线结构，并且根据n线标定原理计算出斜线斑点在体模中的空间位
[0039][0040]
对于由同一条斜线和不同组合的平行线构成的n线结构，得到不同的空间位置设计加权平均计算公式
[0041][0042]
其中，d(
·
)为斑点的弥散度因子，为两个斑点之间的欧氏距离，对所有斜线斑点do计算对应体模空间中的坐标点do，组成集合s(do)和s(do)，通过随机采样一致性ransac方法，计算超声平面到标定体模的初始转换矩阵
[0043]
对于集合中每条斜线wo对应m种映射使用转换矩阵计算斜线wo与超声平面的交点定义斜线wo与标定体模的前面板和后面板的交点分别为c
of
和c
ob
，使用的逆矩阵计算出c
of
和c
ob
在超声图像坐标系中的位置u
of
和u
ob
，进而计算出与超声平面的交点，在所有候选斑点{do}m中，选择与交点po欧氏距离最小的斑点do，定义为唯一映射
[0044][0045]
将所有的斜线映射存入最终的斜线映射集合mo。
[0046]
优选地，在所述步骤(2)中，对于多个位置的超声图像和跟踪数据，通过定义平面内误差ipe的损失函数优化标定矩阵并通过非线性迭代优化方法进行全局优化，公式为
[0047][0048]
其中，‖p
i-di‖为连线在超声平面的交点pi和对应斑点di之间的欧氏距离；n为标定图片数量；nw为图像中线映射数量，计算交点需要使用优化过程中的逆矩阵，图像的跟踪数据和的逆矩阵，以及体模标定矩阵非线性优化的目标被定义为最小化所有超声图像的平面误差之和，公式为
[0049][0050]
使用powell迭代的非线性优化方法计算最小化损失函数，采用所有超声图像的闭式结果作为初始值。
[0051]
优选地，在所述步骤(2)中优化时，使用8个参数进行优化，参数包括3个旋转参数、3个平移参数、2个缩放参数。
[0052]
本领域普通技术人员可以理解，实现上述实施例方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，所述的程序可以存储于一计算机可读取存储介质中，该程序在执行时，包括上述实施例方法的各步骤，而所述的存储介质可以是：rom/ram、磁碟、光盘、存储卡等。因此，与本发明的方法相对应的，本发明还同时包括一种超声探头自动空间标定装置，该装置通常以与方法各步骤相对应的功能模块的形式表示。该装置包括：
[0053]
扩展n线标定体模，其配置来通过在横向和纵向上排布多个n线结构，增加超声图像中的有效基准点数量；扩展n线标定体模的导线孔为三角形，设计多层并且每层上有两个以上的n线结构，一种线为n的斜线，另一种线为n的平行线，将三维n线连接方式投影在前面板上，根据前后面板交点位置定义每条线的编号；在设计连线方式上，上下相邻两层的结构不重复，斜线两端点的水平方向距离小于30mm，相邻层之间的宽度为10mm，共有12层n线，模型中共有51根平行线和33根斜线，采用集合型中共有51根平行线和33根斜线，采用集合表示体模的l层线集，其中第i层线集中包含：个平行线w
p
、个斜线wo；
[0054]
快速标定模块，其配置来基于扩展n线标定体模超声图像进行快速标定：通过斑点分割、平行线限制和基准点识别的步骤准确识别所有斑点编号，并通过非线性迭代优化方法计算出超声探头标定矩阵。
[0055]
优选地，所述扩展n线标定体模使用直径为0.234mm的尼龙线安装所有n线结构。
[0056]
以下更详细地说明本发明。
[0057]
探头标定过程涉及的组件、连接和坐标系转换关系如图1所示。使用标定体模进行闭式求解过程可以表示为
[0058][0059]
为了获得超声探头矩阵(红色虚线转换)，首先需要计算出转换矩阵使得超声图像中的像素位置[u,v,0,1]
t
能够完美匹配到体模中的空间位置[x,y,z,1]
t
。对于和矩阵是可以通过电磁跟踪设备可以直接获得，矩阵需要使用额外的探针标定。
[0060]
标定体模由前后面板和侧板组成，前后面板上等间距排布相同位置的孔洞。使用细线连接前后板上的不同孔洞，组成“n”线结构。其中，一种线为“n”的斜线，另一种线为“n”的平行线。为了直观显示连接关系内部结构，将三维n线连接方式投影在前面板上，如图2所示。根据前后面板交点位置，可以定义每条线的编号，例如，第一行最左侧的平行线(红线)编号分别为“c1-c1”，第七行最左侧的斜线(黄线)编号分别为“d7-g7”，以此类推。
[0061]
在设计连线方式时候，尽可能保证上下相邻两层的结构不重复，以方便识别图像中每个斑点的编号。此外，斜线两端点的水平方向距离小于30mm，使得斜线斑点成像尽可能为圆形。同时，为了保证在垂直方向上n线斑点的成像质量，相邻层之间的宽度为10mm，共有12层n线。模型中共有51根平行线和33根斜线。采用集合表示体模的l层线集，其中第i层线集中包含个平行线w
p
和个斜线wo。最后我们使用直径为0.234mm的尼龙线安装所有n线结构。
[0062]
在使用超声探头扫描上述扩展n线体模的时候，会得到包含大量n-fiducial斑点的超声图像。一个主要的困难是如何在存在散斑噪声的情况下准确而稳健地识别n-fiducial斑点，因为他们具有相似的图像强度和形状。类似于双n线体模，可以发现扩展n线体模的超声图像中存在三个独特几何特征：
[0063]
·
不同层n线交点的所在直线几乎平行；
[0064]
·
每层n线的交点中，至少有三个点共线；
[0065]
·
平行线(红线)和任意位置超声平面的交点集与平行线的二维投影点集(红点)之间都存在等比变换。
[0066]
基于上述的三个几何特征，提出了基于扩展n线体模超声图像的快速标定方法，主要分为四个部分，如图3所示。第一部分，通过形态学操作和像素聚类的方式提取斑点的特征，包括位置、形态和强度。在第二部分，通过平行线的限制将斑点分类，得到每一层上的斑点集合。之后，提出了一种基于点云增量配准的方法，将斑点和体模中的连线关联，得到平行线映射集合和斜线映射集合。最后，结合跟踪数据，使用非线性优化方式计算标定矩阵。
[0067]
斑点分割
[0068]
为了分割出有效斑点，首先使用一种形态结构算子对原始超声图像进行开运算操作，去除超声图像中的大斑点，再利用另一种形态结构算子去除小斑点。注意，经过两次形态学操作之后，此处图像中保留的斑点不仅包括所需的n基准点，还包括类似形状和大小的噪点。之后通过像素聚类的方法得到全部的斑点区域。
[0069]
根据探头的成像原理，图像横向分辨率由超声横向波束宽度决定，并随轴向深度而变化。点状小物体被错误表示并在图像中显示为短线(除非在扫描平面焦深处)。为了全面描述斑点的形态，使用点扩散函数(psf)，即二元高斯函数拟合斑点形状，如下所示：
[0070][0071]
其中(xi,yi)是散斑中第i个像素的位置，ii是第i个像素的强度。考虑到散斑的强度分布，该方法中μ和δ分别设置为255和10。之后，对所有聚类像素的强度呈线性总和的散斑进行评分。这种方法有效地驱动了一个点的大小和图像强度来影响它的分数，这个度量将在线条拟合和散斑匹配的后续阶段找到它的用途。
[0072]
平行线约束
[0073]
完成斑点分割和评分后，需要利用n基准点之间的特殊几何关系：不同层n线交点的所在直线几乎平行。由于超声探头在采集图像过程中，在水平方向上不会有非常大的倾斜，所以只统计两个候选点连线与水平方向夹角在-20度到20度之间的直线簇。需要在这些直线簇中找出包含最多直线的角度区间，定义为主方向角度。最简单的方法是对角度排序之后，用一个小窗口搜索，找到角度最集中的位置为平行线主方向。在实验中，选择角度窗口大小为2度，这已被证明足以找到准确的主方向。
[0074]
之后，统计所有保留直线簇的截距，使用具有噪声的基于密度的聚类方法dbscan得到截距的聚类集合。考虑到n基准点的第二个独特关系，只保留斑点数量大于3的集合。利用每个集合中的斑点位置和评分计算当前集合的直线方程，并根据直线截距进行排序，从而识别到不同层n基准点所在直线。用d＝{d
ij
,i∈[1,m],j∈[1,ni]}表示图中检测出m层斑点集合，第i层斑点集合包含有ni个斑点d。值得注意的是，在每层n线斑点集合中仍然存在位于直线上的斑点噪声。
[0075]
斑点识别
[0076]
考虑到n基准点的第三个独特关系，在平行线和任意位置超声平面的交点集合与
平行线的二维投影点集合之间都存在等比变换。理论上，可以直接通过随机采样一致性ransac方法将斑点集合d和体模连线集合w进行配准，但是这种方式迭代时间长且随机性比较大。为了提高配准效率和精度，我们提出了一种基于点云增量配准的斑点识别算法，伪代码见表1。
[0077]
表1
[0078][0079][0080]
首先，为了快速检索超声图像中斑点，对斑点集合d构建了检索kd树。之后，计算体模坐标系到图像坐标系的初始旋转矩阵
[0081][0082]
其中，α为第一层斑点所在直线方程的角度；s为图像间隔，可以通过前两层斑点的
直线方程和体模配置计算。穷举了斑点集合d中前两层斑点和集合w中前两层平行线之间的所有匹配《w
p
,d》。对于每一组匹配，计算平移向量t＝d-r*p，p为平行线w
p
在前面板上的端点，最终得到初始配准矩阵t
2d
＝[r|t]。之后，从线集w第一层开始，使用配准矩阵t计算出每条平行线w
p
在超声图像中投影点p，在kd树中搜索所有可能匹配的斑点d。本文定义以下匹配相似度
[0083][0084]
其中，v(d)为斑点d的形态评分；‖p-d‖为投影点p和斑点d之间的欧氏距离；s为配准后图像间距。该测度结合了斑点的结构特征以及欧式距离，减少微小斑点的干扰。将平行线w
p
与相似度最大的斑点组成映射关系《w
p
,d
p
》，存入平行线映射集合中m
p
。当计算到第三行之后，需要根据m
p
中的映射关系，重新计算配准矩阵t
2d
，然后重新计算所有平行线的最佳映射关系，存入到m
p
中。完成所有平行线映射之后，计算斜线的映射关系。使用配准矩阵t
2d
计算出每条斜线wo在超声图像中的投影线段l。在kd树中搜索所有位于线段l的斑点{do}m，组成一对多的映射关系《wo,{do}m》，加入候选斜线映射集合不同于一个平行线只能对应一个斑点，一个斜线可能对应多个斑点，所以集合并不是最终的斜线映射集合。在多次循环之后，保留集合m
p
和集合中总映射关系最多的集合。
[0085]
之后，还需要根据n线几何结构，筛选中每条斜线与斑点的唯一映射《wo,do》。记集合m
p
中第i层存在个平行线映射集合集合中第i层任意斜线wo对应m种映射通过组合任意两条平行线和一条斜线wo，都可以构成一组n线结构，并且可以根据n线标定原理(carbajal et al.2013)计算出斜线斑点在体模中的空间位置
[0086][0087]
对于由同一条斜线和不同组合的平行线构成的n线结构，可以得到不同的空间位置根据n线标定原理，两条平行线的空间距离越大，抗干扰性会越强。同时考虑到斑点形态，设计加权平均计算公式
[0088][0089]
其中，d(
·
)为斑点的弥散度因子；为两个斑点之间的欧氏距离。对所有斜线斑点do计算对应体模空间中的坐标点do，组成集合s(do)和s(do)。通过随机采样一致性ransac方法，计算超声平面到标定体模的初始转换矩阵
[0090]
对于集合中每条斜线wo对应m种映射使用转换矩阵
计算斜线wo与超声平面的交点定义斜线wo与标定体模的前面板和后面板的交点分别为c
of
和c
ob
。值得注意的是，标定体模坐标系中所有连线与标定体模壁的交点都是已知的，因为标定体模是精确制造的。使用的逆矩阵计算出c
of
和u
ob
在超声图像坐标系中的位置u
of
和u
ob
，进而计算出与超声平面的交点。在所有候选斑点{do}m中，选择与交点po欧氏距离最小的斑点do，定义为唯一映射
[0091][0092]
最后，将所有的斜线映射存入最终的斜线映射集合mo。
[0093]
通过提出的点云增量配准方式，可以准确地匹配超声图像中的斑点与体模中的连线，识别出每个有效斑点的编号。
[0094]
全局优化策略
[0095]
在前面的部分中，已经估计标定体模与超声图像之间的平行线映射集合m
p
和斜线映射集合mo以及转换矩阵根据公式(1)，可以通过闭式解方程得到当前图像到探头传感器的转换通常地，对于多个位置的超声图像和跟踪数据，可以通过定义平面内误差(in-plane error,ipe)的损失函数优化标定矩阵并通过非线性迭代优化方法进行全局优化。
[0096][0097]
其中，‖p
i-di‖为连线在超声平面的交点pi和对应斑点di之间的欧氏距离；n为标定图片数量；nw为图像中线映射数量。计算交点需要使用优化过程中的逆矩阵，图像的跟踪数据和的逆矩阵，以及体模标定矩阵因此，非线性优化的目标被定义为最小化所有超声图像的平面误差之和，公式为
[0098][0099]
使用powell迭代的非线性优化算法最小化损失函数，采用所有超声图像的闭式结果作为初始值。考虑到超声图像的各向异性缩放(图像平面中的水平和垂直缩放不同)，使用8个参数(3旋转+3平移+2缩放)进行优化。
[0100]
以上所述，仅是本发明的较佳实施例，并非对本发明作任何形式上的限制，凡是依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属本发明技术方案的保护范围。