CN105616003B - 一种基于径向样条插值的软组织三维视觉跟踪方法 - Google Patents

Abstract

本发明公开了一种基于径向样条插值的软组织三维视觉跟踪方法，在基于软组织图像的三角形目标区下，构造了一种低复杂度的非线性空间的三维形变模型，在每一帧图像下通过快速求解出三维形变模型的最优参数，实现对软组织图像的目标区域实时、准确的三维视觉跟踪，具有简单、灵活等特点，符合当前临床医疗发展的需要。

Description

一种基于径向样条插值的软组织三维视觉跟踪方法

技术领域

本发明属于视觉跟踪技术领域，更为具体地讲，涉及一种基于径向样条插值的软组织三维视觉跟踪方法。

背景技术

近年来，以内窥镜手术为代表的微创手术越来越多的用于临床医疗。在这类手术中，医生通过观察探入内窥镜反馈回的实时影像，操控微创手术器械完成复杂的心脏外科手术。

以脱泵冠脉搭桥术(off pump CABG)为例，这种心脏不停跳的手术可使患者摆脱人工心肺机，降低体外循环对患者身体的损害，降低并发症的风险，缩短患者的康复时间。但是由于直接视觉的缺失，手术医生无法准确、直观的感知跳动心脏表面手术区域的三维运动，这给完成复杂精细的手术操作带来困难。

为了给手术医生提供实时、准确的软组织三维运动信息，国外学者尝试利用各种传感器跟踪软组织的三维运动。由于微创手术空间狭小，无法植入大型三维扫描设备，人体软组织结构复杂又含较多的非线性形变，软组织三维跟踪成为医学工程领域中的难点。在不需要额外测量设备，基于立体内窥镜的软组织三维视觉跟踪逐渐受到科研和临床界的重视。W.Lau等在文献《Stereo-based endoscopic tracking of cardiac surfacedeformation》中利用B样条模型基于立体内窥镜影像对心脏表面形变进行三维跟踪，D.Stoyanov等在文献《A practical approach towards accurate dense 3-D depthrecovery for robotic laparoscopic surgery》中则使用分片双线性映射模型在机器人辅助内窥镜手术中对心脏表面进行三维重构和跟踪，R.Richa等在文献《Three-dimensional motion tracking for beating heart surgery using a thin-platespline deformable model》中使用薄板样条模型跟踪跳动心脏表面的运动。

然而，上述方法中所采用的形变模型通常较为复杂，很难实时应用；同时由于模型较复杂，在每一帧模型参数求解时，往往无法获得最优的参数解，甚至导致迭代求优算法无法收敛，最终导致跟踪失败。而现有的简单形变模型，如线性仿射变换模型，在用于跟踪时通常可以满足实时性要求，但由于无法准确描述软组织表面的非线性形变，因而无法获得精确的三维跟踪结果。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于径向样条插值的软组织三维视觉跟踪方法，通过构造了一种低复杂度的非线性空间形变模型，可以在每一帧快速的求解出模型参数，实现对软组织目标区域实时、准确的三维视觉跟踪。

为实现上述发明目的，本发明一种基于径向样条插值的软组织三维视觉跟踪方法，其特征在于，包括以下步骤：

(1)、将立体内窥镜获取的首帧k＝1左图像记为I_L1；

(2)、确定三角形目标区

在左图像I_L1中，确定需要跟踪的三角形目标区，其三个顶点的像素坐标分别用二维列向量表示为：m_a，m_b和m_c；在三角形目标区内包含的N个像素点，其坐标分别记为：m₁,m₂,...,m_N；

(3)、保存跟踪模板

从左图像I_L1中提取目标区内N个点的像素值，按顺序排列并保存为N维行向量，记为：

T＝[I_L1(m₁) I_L1(m₂) ... I_L1(m_N)]

其中，I_L1(m_n)为左图像I_L1中像素点m_n(n＝1,2,…,N)的像素值；

(4)、根据步骤(2)确定的三角形目标区，构造径向样条插值三维形变模型；

(4.1)、构造模型的设计矩阵

其中，三维列向量λ_n(n＝1,2,...,N)为m_n在三角形目标区中的重心坐标，由下式获得：

λ_o为三角形目标区外接圆圆心的重心坐标，由下式获得：

s_n(n＝1,2,…,N)为像素点m_n的径向样条插值系数，由下式计算：

其中，m_o＝[m_a m_b m_c]λ_o为三角形目标区外接圆圆心的像素坐标；

(4.2)、构造径向样条插值三维形变模型

P＝[p_a p_b p_c p_o]Q

其中，模型参数p_a、p_b、p_c、p_o分别为三角形目标区的三个顶点m_a、m_b、m_c，和外接圆圆心m_o在真实三维空间中的三维坐标列向量；3×N维的输出矩阵P包含了三角形目标区N个像素点的三维坐标，即P＝[p₁ p₂ … p_N]，其中p_n为像素点m_n对应的三维坐标；

(5)、获得第k帧的左右图像

将立体内窥镜获取的第k帧左、右图像分别记为I_Lk、I_Rk；(6)、基于步骤(4)构造的三维形变模型，求解最优的模型参数，使左、右图像I_Lk、I_Rk与模板T之间最匹配，可用以下方程表示为：

其中，为最优的模型参数，它们使目标函数||I_L-T||²+||I_R-T||²最小；

I_L＝[I_Lk(m′₁) I_Lk(m′₂) ... I_Lk(m′_N)]、I_R＝[I_Rk(m″₁) I_Rk(m″₂) ... I_Rk(m″_N)]分别为三角形目标区内N个点在第k帧左、右图像中的像素值排列而成的行向量，其中m′_n和m″_n为m_n在第k帧左、右图像中的新坐标；

(7)、输出目标区三维形态

将步骤(6)获得的模型参数代入模型方程，得：

输出矩阵P^*包含了目标区N个点在第k帧时刻的三维坐标，从而获得目标区三维形态；

(8)、k＝k+1：等待下一帧图像到来，帧序号加1，然后返回步骤(5)。

进一步的，所述的步骤(6)中，像素点m_n在第k帧左、右图像中的新坐标m′_n和m″_n的计算方法为：

其中，p_n为像素点m_n对应的三维坐标，C_L和C_R分别为立体内窥镜左、右相机的投影矩阵，Φ是从三维齐次坐标到二维像素坐标的投影函数，具体定义如下：

其中，代表任意的三维向量。

本发明的发明目的是这样实现的：

本发明基于径向样条插值的软组织三维视觉跟踪方法，在基于软组织图像的三角形目标区下，构造了一种低复杂度的非线性空间的三维形变模型，在每一帧图像下通过快速求解出三维形变模型的最优参数，实现对软组织图像的目标区域实时、准确的三维视觉跟踪，具有简单、灵活等特点，符合当前临床医疗发展的需要。

同时，本发明基于径向样条插值的软组织三维视觉跟踪方法还具有以下有益效果：

(1)、运算量低、实时性好：由于所构造的径向样条插值模型仅有四个控制点，即三个顶点加外接圆圆心，共3×4个模型参数，求解简单快速，因此可以在每一帧实时求解模型参数，实现对目标区的实时三维跟踪。

(2)、跟踪准确性好：尽管所采用的形变模型较简单，但其包含线性仿射变换和非线性径向三次样条插值两部分；在第二部分使用的三次样条曲线可以较好的拟合软组织表面的非线性特性；因此可以获得较好的跟踪结果。

(3)、模型参数具有实际物理意义，可在每一帧直接获得目标区三个顶点和外接圆圆心的三维坐标，方便对跟踪结果的直接使用。

(4)、采用三角形目标区，目标区划定更灵活，容易避开易对跟踪产生不利影响的无特征或动态特征区域，以及容易被手术器械遮挡的区域。

附图说明

图1是本发明基于径向样条插值的软组织三维视觉跟踪方法流程示意图；

图2是图1中构造径向样条插值模型的原理示意图；

图3是本发明实施例的三维跟踪结果示意图。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

实施例

本发明的基本思想是通过构造一种简单的空间形变模型快速匹配软组织目标区在当前帧的三维形态，实现实时的三维视觉跟踪。所构造的形变模型包含两部分：第一部分为线性仿射变换，第二部分为非线性的三次样条插值。尽管模型简单、复杂度低，仅需求解四个(三个顶点加外接圆圆心)控制点的三维坐标(共12个模型参数)，但是由于利用了三次样条函数匹配软组织表面的非线性形变，因此可以实现对软组织目标区的实时、精确三维跟踪。

下面结合附图对本发明进行详细说明，具体如下：

图1是本发明基于径向样条插值的软组织三维视觉跟踪方法流程示意图。

在本实施例中，如图1所示，本发明一种基于径向样条插值的软组织三维视觉跟踪方法具体包括以下流程：

S1、获取k＝1帧左图像：将立体内窥镜获取的首帧(k＝1)左图像记为I_L1。

S2、确定三角形目标区：在左图像I_L1中，确定需要跟踪的三角形目标区，其三个顶点的像素坐标分别用二维列向量表示为：m_a，m_b和m_c；在三角形目标区内包含的N个像素点，其坐标分别记为：m₁,m₂,...,m_N。

S3保存跟踪模板：从图像I_L1中提取目标区N个像素点的像素值，按其位置顺序(如从左至右或由上而下)排列并保存为N维行向量，记为：

T＝[I_L1(m₁) I_L1(m₂) ... I_L1(m_N)]

其中，I_L1(m_n)为左图像I_L1中像素点m_n(n＝1,2,…,N)的像素值。

S4、构造径向样条插值模型：根据步骤S2确定的目标区，构造径向样条插值三维形变模型；

首先，构造模型的设计矩阵

其中，三维列向量λ_n(n＝1,2,...,N)为m_n在三角形目标区中的重心坐标，由下式获得：

λ_o为三角形目标区外接圆圆心的重心坐标，由下式获得：

s_n(n＝1,2,…,N)为像素点m_n的径向样条插值系数，由下式计算：

其中，m_o＝[m_a m_b m_c]λ_o为三角形目标区外接圆圆心的像素坐标；

其次，构造径向样条插值三维形变模型

P＝[p_a p_b p_c p_o]Q

其中，模型参数p_a、p_b、p_c、p_o分别为三角形目标区的三个顶点m_a、m_b、m_c，和外接圆圆心m_o在真实三维空间中的三维坐标列向量；3×N维的输出矩阵P包含了三角形目标区N个像素点的三维坐标，即P＝[p₁ p₂ … p_N]，其中p_n为像素点m_n对应的三维坐标；

由步骤S4中的模型方程可知，一旦模型参数p_a、p_b、p_c、p_o确定，即目标区三个顶点和外接圆圆心的三维坐标确定，则目标区内其它像素点的三维坐标，便由三维形变模型唯一确定。因此，在后续图像帧中，对目标区的三维跟踪，实际就简化为在每一帧中求解最优的模型参数。一旦模型参数确定，目标区的三维形态就得以确定。

S5、获得第k帧左右图像：将立体内窥镜获取的第k帧左、右图像分别记为I_Lk、I_Rk。

S6、求解形变模型参数：基于步骤S4构造的三维形变模型，求解最优的模型参数，使左、右图像I_Lk、I_Rk与模板T之间最匹配，可用以下方程表示为：

其中，为最优的模型参数，它们使目标函数||I_L-T||²+||I_R-T||²最小；

I_L＝[I_Lk(m′₁) I_Lk(m′₂) ... I_Lk(m′_N)]、I_R＝[I_Rk(m″₁) I_Rk(m″₂) ... I_Rk(m″_N)]分别为三角形目标区内N个点在第k帧左、右图像中的像素值排列而成的行向量，其中m′_n和m″_n为m_n在第k帧左、右图像中的新坐标，通过下式计算：

其中，p_n为像素点m_n对应的三维坐标，C_L和C_R分别为立体内窥镜左、右相机的投影矩阵，Φ是从三维齐次坐标到二维像素坐标的投影函数，具体定义如下：

其中，代表任意的三维向量。

本步骤中的最优参数求解，为常规的无约束最优化问题，可通过最速下降法、拟Newton法或共轭梯度法等经典的最优化方法求解。

S7、输出目标区三维形态：将步骤(6)获得的模型参数代入模型方程，得：

输出矩阵P^*包含了目标区N个点在第k帧时刻的三维坐标，从而获得目标区三维形态。

S8、k＝k+1：等待下一帧图像到来，帧序号加1，然后返回步骤(5)。

图2是图1中构造径向样条插值模型的原理示意图。

在本实施例中，如图2所示，径向样条插值三维形变模型分解成两部分：第一部分是线性的仿射变换，如图2(a)所示；第二部分是非线性的径向样条插值，如图2(b)所示。以下模型方程的推导中，在步骤S4的叙述中已定义或说明过的变量在此不再赘述。

第一部分，线性仿射变换模型用方程表示为：

其中，3×N维的输出矩阵包含了目标区N个像素点仿射变换后的三维坐标，为像素点m_n(n＝1,2,...,N)仿射变换后的三维坐标；λ_n(n＝1,2,...,N)为像素点m_n在三角形目标区中的重心坐标。由仿射变换原理可知：仿射变换后的目标区仍为一平面三角形，其形态完全由三个顶点m_a、m_b、m_c确定。

第二部分，非线性形变插值，用于匹配目标区的非线性形变。设目标区外接圆圆心由点形变至p_o，形变量为其中，和p_o均代表三维坐标值；目标区内其它像素点的形变量可据其到外接圆圆心的像素距离通过三次样条插值获得，具体方程为：

其中，3×N维的输出矩阵包含了目标区N个像素点的非线性形变量；为目标区外接圆圆心经第一部分仿射变换后的三维坐标(即未发生非线性形变前的坐标)，可表示：

由插值系数s_n的计算方程(见步骤S4)可知，在外接圆圆心处(即m_n＝m_o时)，插值系数为1，故外接圆圆心样条插值后的坐标为p_o；而在三个顶点处(即m_n＝m_a，m_b和m_c时)，插值系数为0，故三个顶点在样条插值后的三维坐标仍为p_a，p_b和p_c。由于插值系数是通过像素点m_n到外接圆圆心m_o之间像素距离的三次样条函数计算获得，三次样条函数的连续性保证输出矩阵构成一光滑的三角形曲面，可以精确贴合软组织表面。

将上述两部分相加，可以得到步骤S4中构造的径向样条插值模型：

图3是本发明实施例的三维跟踪结果示意图。

如图3所示，在本实施例中，在第k＝1帧左图像中，划定三角形目标区，其三个顶点的像素坐标分别为：

目标区内共包含N＝3281个像素点，外接圆圆心的像素坐标和重心坐标分别为：

列出部分帧的跟踪结果如下：

在第1帧，求解出形变模型参数(即三个顶点和外接圆圆心的三维坐标)分别为：

目标区在左、右图像中的位置和在三维空间中的形态如图第一行所示。

在第11帧，求解出形变模型参数分别为：

当前时刻(帧)目标区在左、右图像中的位置和在三维空间中的形态如图第二行所示。

在第21帧，求解出形变模型参数分别为：

当前时刻(帧)目标区在左右图像中的位置和在三维空间中的形态如图第三行所示。

在第31帧，求解出形变模型参数分别为：

当前时刻(帧)目标区在左右图像中的位置和在三维空间中的形态如图第四行所示。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (3)

1.一种基于径向样条插值的软组织三维视觉跟踪方法，其特征在于，包括以下步骤：

(1)、将立体内窥镜获取的首帧k＝1左图像记为I_L1；

(2)、确定三角形目标区

在左图像I_L1中，确定需要跟踪的三角形目标区，其三个顶点的像素坐标分别用二维列向量表示为：m_a，m_b和m_c；在三角形目标区内包含的N个像素点，其坐标分别记为：m₁,m₂,...,m_N；

(3)、保存跟踪模板

从左图像I_L1中提取目标区内N个点的像素值，按顺序排列并保存为N维行向量，记为：

T＝[I_L1(m₁) I_L1(m₂) ... I_L1(m_N)]

其中，I_L1(m_n)为左图像I_L1中像素点m_n(n＝1,2,…,N)的像素值；

(4)、根据步骤(2)确定的三角形目标区，构造径向样条插值三维形变模型；

(4.1)、构造模型的设计矩阵

其中，三维列向量λ_n(n＝1,2,...,N)为m_n在三角形目标区中的重心坐标，由下式获得：

λ_o为三角形目标区外接圆圆心的重心坐标，由下式获得：

s_n(n＝1,2,…,N)为像素点m_n的径向样条插值系数，由下式计算：

其中，m_o＝[m_a m_b m_c]λ_o为三角形目标区外接圆圆心的像素坐标；

(4.2)、构造径向样条插值三维形变模型

P＝[p_a p_b p_c p_o]Q

其中，模型参数p_a、p_b、p_c、p_o分别为三角形目标区的三个顶点m_a、m_b、m_c，和外接圆圆心m_o在真实三维空间中的三维坐标列向量；3×N维的输出矩阵P包含了三角形目标区N个像素点的三维坐标，即P＝[p₁ p₂ … p_N]，其中p_n为像素点m_n对应的三维坐标；

(5)、将立体内窥镜获取的第k帧左、右图像分别记为I_Lk、I_Rk；

(6)、基于步骤(4)构造的三维形变模型，求解最优的模型参数，使左、右图像I_Lk、I_Rk与模板T之间最匹配，可用以下方程表示为：

其中，为最优的模型参数，它们使目标函数||I_L-T||²+||I_R-T||²最小；

I_L＝[I_Lk(m₁′) I_Lk(m₂′) ... I_Lk(m_N′)]、I_R＝[I_Rk(m₁″) I_Rk(m₂″) ... I_Rk(m_N″)]分别为三角形目标区内N个点在第k帧左、右图像中的像素值排列而成的行向量，其中m_n′和m_n″为m_n在第k帧左、右图像中的新坐标；

(7)、输出目标区三维形态

将步骤(6)获得的模型参数代入模型方程，得：

输出矩阵P^*包含了目标区N个点在第k帧时刻的三维坐标，从而获得目标区三维形态；

(8)、k＝k+1：等待下一帧图像到来，帧序号加1，然后返回步骤(5)。

2.根据权利要求1所述的基于径向样条插值的软组织三维视觉跟踪方法，其特征在于，所述的步骤(6)中，像素点m_n在第k帧左、右图像中的新坐标m_n′和m_n″的计算方法为：

其中，p_n为像素点m_n对应的三维坐标，C_L和C_R分别为立体内窥镜左、右相机的投影矩阵，Φ是从三维齐次坐标到二维像素坐标的投影函数，具体定义如下：

其中，代表任意的三维向量。

3.根据权利要求1所述的基于径向样条插值的软组织三维视觉跟踪方法，其特征在于，所述的径向样条插值三维形变模型包括线性仿射变换模型和非线性的径向样条插值；

所述的线性仿射变换模型为：

其中，3×N维的输出矩阵包含了目标区N个像素点仿射变换后的三维坐标，为像素点m_n(n＝1,2,...,N)仿射变换后的三维坐标；

所述的非线性的径向样条插值为：设三角形目标区的外接圆圆心由点形变至p_o，形变量为目标区内其它像素点的形变量可据其到外接圆圆心的像素距离通过三次样条插值获得，具体方程为：

其中，s_n(n＝1,2,…,N)为像素点m_n的径向样条插值系数，3×N维的输出矩阵包含了三角形目标区N个像素点的非线性形变量；为目标区外接圆圆心经第一部分仿射变换后的三维坐标可表示：

其中，λ_o为三角形目标区外接圆圆心的重心坐标；

将线性仿射变换模型与非线性的径向样条插值相加，得到径向样条插值三维形变模型P：