-
导航系统一般由软件系统及主动或被动跟踪定位系统组成，而光学跟踪定位系统因其可靠、非接触等特点在导航系统中广泛应用。光学跟踪定位系统一般应用在手术导航、精准焊接、工件定位、打磨抛光等领域。例如在手术导航方面，光学跟踪定位系统可以通过实时跟踪手术工具和患者实体，然后将手术工具与患者实体之间的相对姿态映射到虚拟导航界面，从而使手术更加便捷，而且也减少了创伤面积[1-3]。为了减少可见光的干扰，手术导航一般使用近红外跟踪定位系统，由两台高速相机搭建而成，通过跟踪被动标记球而实时重建空间坐标，达到精准导航的目的。
跟踪定位技术是衡量手术导航系统的关键因素。高精度和实时性是手术导航成功的先决条件。目前，在手术导航系统中应用的主要跟踪方法是电磁跟踪和光学跟踪。电磁导航系统虽然成本低，使用灵活，但对金属物体特别是铁磁仪器特别敏感，其精度对电磁波的干扰很敏感。光学跟踪系统（OTS）不仅定位精度高，而且使用方便[4]。它迅速发展成为当前手术导航的主要产品，目前，主流商用产品如NDI的Polaris VEGA 系列，最多可跟踪25个工具，跟踪精度仅为0.12 mm。
北京理工大学首先通过并行算法在FPGA中实现了三维点重建定位附着在仪器上的标记，然后根据这些标记的几何关系估计姿势[5]。极线约束通常用于立体匹配，但它无法确定与两个光学中心共面的标记的正确对应，从而导致重影标记在重建3D点后生成[6]。为了解决重影标记的形成，学者研究了其他改进的方法，如标记点与极线之间的几何距离约束的方法[7]；将最大轮廓重投影一致性和极线约束结合到目标函数中的方法[8]；根据标记提供的几何信息查找标记投影之间的对应关系，以及基于投影的大小在帧间快速跟踪2D图像中的标记方法[9]；基于标记点几何排列的仪器识别方法和基于标记点运动矢量的仪器跟踪方法[10]。
对于新手外科医生来说，计算机辅助手术的术前培训必不可少，旨在提高外科医生的动手实践能力。然而，Polaris VEGA 相关产品不仅价格昂贵，而且其所依赖的技术受专利保护，不利于二次开发和定制化，因此，通常用于临床手术，不能广泛用于培训新手外科医生。
文中基于双目视觉原理，设计了一种低成本、高精度的光学跟踪系统。在设计的光学跟踪系统中，6个被动标记球用作安装在器械上的标记。该系统使用两个高速相机来捕捉标记的图像，两个红外灯用作主动光源，为了增强标记亮度减弱可见光的干扰，在每个相机镜头前添加了近红外滤光片。为提取标记中心的像素坐标，采用轮廓过滤方法消除干扰轮廓，然后用最小二乘椭圆拟合算法获取中心像素坐标。这种图像处理方法大大增强了抗干扰能力，提高了系统的精度。接下来，器械识别算法根据器械特有几何特征对其正确识别并有序储存器械上标记球的像素坐标序列，用于匹配左右摄像机图像中的标记。最后，通过标记的3D坐标重建，得到器械尖端的坐标，并通过实验验证了系统的准确性和稳定性，可以满足手术导航的需要。
-
一般来说，可见光相机在可见光条件下工作，但它也可以检测近红外光。量子效率 (QE) 是器件对各个波长灵敏度的重要表征，该系统使用的CMOS USB3.0相机，它在不同波长下QE的分布如图1（a）所示。由图1（a）可知，可见光波段的QE大于近红外波段的QE。如果相机直接用于跟踪器械，手术室内环境光的干扰是不可避免的，而且标记球的亮度会很暗，如图1（b）所示，这将严重影响器械尖端定位的准确性。
图 1 传感器测试
Figure 1. Sensor test
因此，在工业定焦镜头前添加了一个中心波长为850 nm的窄带通滤波器。由图1（c）可知，添加近红外滤光片可以有效阻止可见光等进入传感器。
整体系统如图2所示，硬件部分使用两台同型号高频黑白相机（MV-CA013-21UM，HIKROBOT，China），相机分辨率为1280×1024，最大帧率为210 fps，二者间距420 mm。相机内侧为两个中心波长为850 nm的近红外照明灯。镜头采用定焦镜头，焦距为8 mm，镜头前添加有近红外滤光片。软件部分可以实时采集处理得出标记中心像素坐标及器械尖点在世界坐标系中的坐标。
图 2 光学跟踪系统
Figure 2. Sensor acquisition
系统原理框图如图3所示，当系统启动时，近红外光源照射固定在器械上的被动标记球的表面，双目视觉系统收集标记球反射的光线在相机的传感器上成像，并将图像数据通过USB3.0接口发送到计算机，软件可以实时显示标记及器械尖点的世界坐标。软件算法处理主要分为三步：首先，应用标记轮廓过滤算法得到有效轮廓并利用最小二乘椭圆拟合算法获取标记投影中心的像素坐标，结合双目视觉系统的内外参数，对捕获的中心坐标进行畸变校正；其次，根据器械的几何排列对左右视图的器械逐一识别并有序储存器械上标记球的像素坐标序列，为立体匹配做好基础；最后，通过已标定的双目相机内外参数重建标记中心的3D坐标。更重要的是，结合器械尖端在器械坐标系中的坐标，可以计算出尖端的世界坐标。从而实现对器械的实时定位和跟踪。
图 3 系统原理框架
Figure 3. System principle framework
-
如图4所示，文中设计的器械安装六个NDI公司生产的被动标记球，这种标记球对近红外光具有高反射率，其组成一个等腰直角三角形，三角形每条边都有三个标记球，但是处于每条边中间的标记球的位置会不同，比如有些位于中点，而有些位于三等分点。尖点位于图4最右侧，通过尖点旋转验证系统抗干扰性。D16、D14、D24这三组距离的准确性也会在动态跟踪实验中验证 。
图 4 标记器械
Figure 4. Marking instrument
-
Zhang’s标定法[11]近年来在相机标定中得到广泛应用，并且已经成熟。它是传统相机标定和自标定的折衷方法，这种方法是稳健的，无需开发高精度定位参考。更重要的是，标定精度可以满足要求，因此，文中采用此方法来标定相机，标定后将得到左右相机的旋转矩阵、平移矩阵以及每个相机的内参矩阵和畸变参数。
相机光学系统中透镜形状引起的畸变称为径向畸变，安装过程中由于不能使透镜和成像面严格平行引起的畸变称为切向畸变。以上两种畸变会使相机坐标系中的归一化坐标
$({({X_c}/{Z_c})_{{\rm{distorted}}}},{({Y_c}/{Z_c})_{{\rm{distorted}}}})$ 与真实归一化坐标$(({X_c}/{Z_c}),({Y_c}/{Z_c}))$ 存在偏差。二者的对应关系为：$$ \begin{split} & X:{{({X_c}} / {{Z_c}}}{)_{{\rm{distorted}}}} = {{({X_c}} / {{Z_c}}})(1 + {k_{11}}{r^2} + {k_{12}}{r^4} + {k_{13}}{r^6}) + \hfill \\ & 2{p_1}{{({X_c}} / {{Z_c}}}){{({Y_c}} / {{Z_c}}}) + {p_2}{({r^2} + 2{{({X_c}} / {{Z_c}}})^2}) \\ \end{split} $$ (1) $$ \begin{split} & Y:{{({Y_c}} / {{Z_c}}}{)_{{\rm{distorted}}}} = {{({Y_c}} / {{Z_c}}})(1 + {k_{11}}{r^2} + {k_{12}}{r^4} + {k_{13}}{r^6}) + \hfill \\ &2{p_2}{{({X_c}} / {{Z_c}}}){{({Y_c}} / {{Z_c}}}) + {p_1}{({r^2} + 2{{({Y_c}} / {{Z_c}}})^2}) \hfill \\ \end{split} $$ (2) 式中：
$ {r^2} = {({{{X_c}} \mathord{\left/ {\vphantom {{{X_c}} {{Z_c}}}} \right. } {{Z_c}}})^2} + {({{{Y_c}} \mathord{\left/ {\vphantom {{{Y_c}} {{Z_c}}}} \right. } {{Z_c}}})^2} $ ；${k_{11}},{k_{12}},{k_{13}}$ 为径向畸变系数，${p_1},{p_2}$ 为切向畸变系数。5个参数均可由相机标定得到。则投影在像素坐标系中的位置变为：
$$ \begin{split} &{u}_{d}={f}_{x}\cdot ({X}_{c}/{Z}_{c}{)}_{{\rm{distorted}}}+{u}_{0}\\ &{v}_{d}={f}_{y}\cdot ({Y}_{c}/{Z}_{c}{)}_{{\rm{distorted}}}+{v}_{0} \end{split} $$ (3) 式中：
${f_x}、{f_y}、{u_0}、{v_0}$ 均可由相机标定的内参矩阵得到。因此，系统首先应通过畸变校正，得到畸变校正后的图像。畸变校正一般通过迭代法解算，当前后两次迭代的差值小于预设阈值，则终止迭代。
-
在光学跟踪系统中，标记中心的像素坐标提取将直接影响器械尖端世界坐标的精度。在捕获的图像中，标记是高亮区域，标记中心的像素坐标是高亮区域的中心，用于计算器械尖端的世界坐标。
然而，因为工作环境中其他物体的高强度反光，形成高亮区域，这种区域会对标记球轮廓的提取造成干扰，因此，文中采用轮廓过滤算法提取标记轮廓。在进行标记轮廓提取之前，必须使用图像分割的方法提取标记边缘。虽然阈值分割可以快速实现边缘提取，但使用单阈值或多阈值的图像分割方法会造成结果不稳定。而Canny边缘检测器作为提取标记边缘的常用方法[12]，算法稳定性已经被多位学者验证，因此，文中使用Canny算法进行边缘检测。边缘提取后使用OpenCV中轮廓检测函数获得轮廓集
${C^d} = \{ c_i^d\} $ ，$c_i^d$ 表示第i个轮廓，$c_i^d{\text{ = \{ }}p_j^d{\text{\} }}$ ，表示第i个轮廓中的单个像素$p_j^d$ 。但这些像素值都是未经校正的，需要利用畸变参数进行校正。然而校正整个原始图像的畸变是很耗时的，因此，只需要校正$c_i^d$ 。最后，得到校正后的轮廓${c_i} = \{ {p_j}\} $ ，以及校正后的最终获得轮廓集$C = \{ {c_i}\} $ 。标记轮廓提取的第一步：判断边缘是否近似于圆。如果
${c_i}$ 是一个圆，那么轮廓上的像素${p_j}$ 到质心$o$ 的距离${d_j}$ 与平均距离$ {o_i} $ 的均方根误差很小。假设${M_i}$ 为${c_i}$ 中的像素点个数，那么$\left| {{d_j} - {o_i}} \right|$ 的均方根误差（RMSE）为：$$ RMS{E_i} = \sqrt {\dfrac{{\displaystyle\sum\limits_{j = 1}^{{M_i}} {{{\left( {{\raise0.7ex\hbox{${{d_j} - {o_i}}$} \mathord{\left/ {\vphantom {{{d_j} - {o_i}} {{o_i}}}}\right.} \lower0.7ex\hbox{${{o_i}}$}}} \right)}^2}} }}{{{M_i}}}} $$ (4) 如果
$ RSM{E_i} $ 小于阈值${t^r}$ ，${c_i}$ 可以被视为一个标记轮廓，如表1所示，分别选取0.4000、0.3000、0.226、0.225、0.224、0.200、0.100等7个阈值进行实验，发现在0.100~0.225区间中，识别数量满足要求，文中${t^r}$ 取值0.224。表 1 tr阈值选取
Table 1. tr threshold selection
Threshold Number 0.400 8 0.300 7 0.226 7 0.225 6 0.224 6 0.200 6 0.100 6 标记轮廓提取的第二步：判断轮廓中的像素数量是否满足一个取值范围，在1.2 m处采集到标记球的像素数量最大值为250个像素，2 m处采集到标记球的像素数量为60个像素，因此，如果轮廓中的像素数量满足此范围，则该轮廓正式被视为一个标记轮廓，如若未执行此步骤，轮廓提取结果将会如图5（a）所示。
图 5 (a)轮廓错误提取; (b)轮廓正确提取
Figure 5. (a) Contour error extraction; (b) Correct contour extraction
最后使用最小二乘椭圆拟合算法[13]过滤后的轮廓
${c_i}$ 拟合得到投影标记中心像素坐标，最终提取结果如图5（b）所示。首先给出椭圆的一般表达式：
$$ f(\vec m,\vec n) = a{x^2} + bxy + c{y^2} + dx + ey + f $$ (5) 式中：矢量
$\vec m = [a\;\;b\;\;c\;\;d\;\;e\;\;f]$ ；$\vec n = [ {{x^2}}\;\;{xy}\;\;{{y^2}} \;\;x \;\;y\;\;1]$ 。令${c_i} = \{ {p_j}\} $ 中像素点的坐标为$({x_j},{y_j})(j = 1,2, \cdots {M_i})$ ，${M_i}$ 为椭圆边界上像素点的个数。引入约束条件
${|| {\vec m} ||^2} = 1$ ，同时，建立如下的目标函数：$$ F(\vec m) = \displaystyle\sum\limits_{j = 1}^{{M_i}} {{f^2}} (\vec m,{\vec n_j}) + M{({|| {\vec m} ||^2} - 1)^2} $$ (6) 式中：
$M$ 为惩罚系数；$ {\vec n_j} =[{{x_j}^2}\;\;\;{{x_j}{y_j}}\;\;\;{{y_j}^2}\;\;\; {{x_j}}\;\;\;{{y_j}}\;\;\;1 ] $ 。最小化目标函数
$F(\vec m)$ ，求得矢量$\vec m$ 的最优解。这是一个非线性最小二乘问题，可以利用牛顿-高斯法或Levenberg-Marquardt算法求解。根据矢量，计算椭圆中心的像素坐标，公式如下：
$$ x{}_c = \dfrac{{2cd - be}}{{{b^2} - 4ac}} $$ (7) $$ y{}_c = \dfrac{{2ae - bd}}{{{b^2} - 4ac}} $$ (8) -
器械识别影响立体匹配的精确性，通过器械识别，可以有序储存器械中每个标记的坐标信息，为立体匹配做好准备。
首先，对左视图进行标记中心提取，得出一组二维像素坐标
$\{ {u_i},{v_i}\} $ ，然后，将像素坐标系转换成图像坐标系：$$ \begin{split} &{x}_{i}=({u}_{i}-{u}_{0})\cdot {d}_{x}\\ &{y}_{i}=({v}_{i}-{v}_{0})\cdot {d}_{y} \end{split} $$ (9) 式中：
${d_x} = {d_y}$ ，表示像元尺寸，文中二者相等。此时，图像坐标系中点集
$\{ {x_i},{y_i}\} $ 与相机坐标系中三维点集具有下述转换关系：$$ \{ {x_i},{y_i}\} \to \{ {x_i},{y_i},f\} $$ (10) 式中：
$\{ {x_i},{y_i},f\} $ 表示图像坐标系中的二维点集与光心之间距离为焦距$f$ 时的相机坐标系空间坐标。假设空间点集的数量为N，子集元素数目为3的集合一共可以找出
$C_N^3$ 个组合。现找出其中的一个子集$\{ {x_i},{y_i},f|i = 0,1,2\} $ ，对每个点列出空间直线方程：$$ x/{x_i} = y/{y_i} = {\textit{z}}/f = t$$ (11) $\{ {x_i},{y_i},f|i = 0,1,2\} $ 对应的$t$ 为1，位于图6中${Z_c} = f$ 平面。在${Z_c} = f$ 平面前后建立${Z_c} = 0.5f$ 与${Z_c} = 4f$ 两个平面，因为每个点的空间直线方程一定经过相机光心，因此，与${Z_c} = 0.5f$ 和${Z_c} = 4f$ 这两个平面存在交点，如图6所示。
图 6 平面拟合
Figure 6. Plane fitting
如图6所示，子集
$\{ {x_i},{y_i},f|i = 0,1,2\} $ 与光心连线的空间直线与${Z_c} = 0.5f$ 和${Z_c} = 4f$ 这两个平面共有6个交点，然后通过这6个交点拟合平面，计算${Z_c} = f$ 平面上$\{ {x_i},{y_i},f|i = 0,1,2\} $ 到该平面的距离平方和${S_{i = 012}}$ ，并判断${S_{i = 012}}$ 是否小于阈值${t^s}$ ，如若满足，则认为$\{ {x_i},{y_i}, f|i = 0,1,2\} $ 在拟合平面上，即该子集属于器械的一条边。最后将其存储到二维向量output中。接下来，依次取出其他子集重复上述过程。如表2所示，分别选取0.012、0.011、0.010、0.009、0.008等5个阈值分别实验，发现在0.008~0.010区间中，识别数量满足要求，文中${t^s}$ 取值为0.009。表 2 ts阈值选取
Table 2. ts threshold selection
Threshold Number 0.012 4 0.011 4 0.010 3 0.009 3 0.008 3 当所有的子集经过上述过程后，不满足的子集被过滤，二维向量output中保存的子集均为目标子集，但是子集之间并不能有效区分，因次，器械设计时对每条边添加了对称性属性，使子集通过对称性区分，分为对称子集与非对称子集。子集对称性区分即：找出output中的一个子集，假设子集是
$\{ {x_i},{y_i},f|i = 0,1,2\} $ 。下面对该子集执行对称性区分操作，首先找出$\{ {x_i},{y_i},f|i = 0,1,2\} $ 中位于边缘两点中间的点，然后计算边缘两点与中间点二者空间向量的弧度值，并用较小的弧度值与较大弧度值相比，该比值记为${\sigma _{i = 012}}$ 。如果${\sigma _{i = 012}}$ 满足公式（12）：$$ {\varsigma _{\min }} < {\sigma _{i = 012}} < {\varsigma _{\max }} $$ (12) ${\varsigma _{\min }} = 0.9$ ，${\varsigma _{\max }}{\text{ = }}1.0$ ，则认为子集$\{ {x_i},{y_i},f|i = 0,1,2\} $ 为对称子集，如图7-a所示；
图 7 (a)对称子集; (b)非对称子集
Figure 7. (a) Symmetric subset; (b) Asymmetric subset
如果满足公式（13）：
$$ {\upsilon _{\min }} < {\sigma _{i = 012}} < {\upsilon _{\max }} $$ (13) ${\upsilon _{\min }}{\text{ = }}0.35$ ,${\upsilon _{\max }}{\text{ = }}0.58$ ，则认为子集$\{ {x_i},{y_i},f|i = 0,1,2\} $ 为非对称子集，如图7（b)所示。因为系统工作距离远大于器械标记球之间的距离，所以
${\sigma _{i = 012}}$ 在对称子集中非常接近1.0；在非对称子集中，${\sigma _{i = 012}}$ 接近器械最初设计的比例，即中间点位于三等分点。依次计算output中的子集，将对称子集按照中间点、边缘点的顺序存储到二维向量Sym中，非对称子集也按照中间点、边缘点的顺序存储到二维向量Asym中。
output中的子集全部处理完之后，子集被分为对称子集与非对称子集，分别存储在二维向量Sym与Asym中。从二维向量Sym中找出一个子集，然后从二维向量Asym中找出两组子集，如果这两组非对称子集满足——存在边缘点与对称子集边缘点相等，则将这三组子集存放到同一个二维向量Tool中。
目前，只是将属于同一器械中的被动标记球中心坐标存入到二维向量Tool中，但并没有明确标记球之间的关系。因此，需要定义一种可以辨别标记球相互关系的方法。
文中定义的方法是以对称子集中间点为起点，对器械中每个标记球按照逆时针顺序编号
$\{ {A_1},{A_2},{A_3}, {A_4},{A_5},{A_6}\} $ ，起点为${A_1}$ ，两个非对称子集的重合点为${A_4}$ ，这两点可以快速确定。但是因空间旋转的不确定性，其他点并不能确定。因此，利用确定的${A_1}$ 与${A_4}$ 的位置关系来确定${A_2}$ ，相互关系如图8（a）所示。
图 8 (a) 相互位置; (b) 器械标记顺序
Figure 8. (a) Mutual position; (b) Instrument marking sequence
由图8（a）可知，当
${A_1}$ 与${A_4}$ 的位置确定之后，${A_2}$ 的位置将唯一确定。例如${A_4}$ 的横纵坐标均大于${A_1}$ 的横、纵坐标，则${A_2}$ 应满足：横坐标小于${A_1}$ 的横坐标而纵坐标大于${A_1}$ 的纵坐标。确定${A_2}$ 之后，对称子集已经确定了两个点。${A_3}$ 则为属于${A_2}$ 与${A_4}$ 那组非对称点集中的中间点，${A_4}$ 为公共点，之前已经确认。${A_5}$ 则为另外一个非对称子集的中间点，${A_6}$ 则为对称子集的另外一个端点。至此，Tool中的6个标记球中心坐标已经有序排列，最后将$\{ {A_1}、{A_2}、{A_3}、{A_4}、{A_5}、{A_6}\} $ 按照下标由小到大的顺序储存，为下一步立体匹配做准备，器械上标记球的编号如图8（b）所示。经过上述过程，器械中的标记球被有序存储。对于每个器械，分别找出
${A_2}、{A_3}、{A_4}、{A_5}、{A_6}$ 计算图8（b）中${{{L_1}} \mathord{\left/ {\vphantom {{{L_1}} {{L_2}}}} \right. } {{L_2}}}$ 与${{{L_3}} \mathord{\left/ {\vphantom {{{L_3}} {{L_4}}}} \right. } {{L_4}}}$ 的比值确定器械种类，文中根据${{{L_1}} \mathord{\left/ {\vphantom {{{L_1}} {{L_2}}}} \right. } {{L_2}}}$ 与${{{L_3}} \mathord{\left/ {\vphantom {{{L_3}} {{L_4}}}} \right. } {{L_4}}}$ 的比值设计了3种器械，如表3所示。表 3 器械种类
Table 3. Types of instruments
${{{L_1}} \mathord{\left/ {\vphantom {{{L_1}} {{L_2}}}} \right. } {{L_2}}}$ ${{{L_3}} \mathord{\left/ {\vphantom {{{L_3}} {{L_4}}}} \right. } {{L_4}}}$ Types of devices > 1 > 1 A < 1 < 1 B > 1 < 1 C 最后，依次从二维向量Sym、Asym取出子集，重复上述过程，确定器械种类。当左视图所有子集处理完成后，再对右视图进行类似处理，最终得到左右视图相同种类器械有序存储的标记球中心坐标，为后续匹配做好准备。
-
器械正确识别后，可以对左右视图中的相同器械按照存储的顺序匹配重建。对于左相机，假设其坐标系的原点与世界坐标系原点重合且无旋转，那么左相机坐标系可以表示为
$o - xy{\textit{z}}$ ，图像坐标系为${O_l} - {X_l}{Y_l}$ ，焦距为${f_l}$ ；右相机坐标系可以表示为$o - {x_r}{y_r}{{\textit{z}}_r}$ ，图像坐标系为${O_r} - {X_r}{Y_r}$ ，焦距为${f_r}$ 。因此，通过相机的透视投影关系可以得到[14]：$$ \left\{ \begin{gathered} {s_l}\left[ {\begin{array}{*{20}{c}} {{X_l}} \\ {{Y_l}} \\ 1 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{f_l}}&0&0 \\ 0&{{f_l}}&0 \\ 0&0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} x \\ y \\ {\textit{z}} \end{array}} \right] \hfill \\ {s_r}\left[ {\begin{array}{*{20}{c}} {{X_r}} \\ {{Y_r}} \\ 1 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{f_r}}&0&0 \\ 0&{{f_r}}&0 \\ 0&0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{x_r}} \\ {{y_r}} \\ {{{\textit{z}}_r}} \end{array}} \right] \hfill \\ \end{gathered} \right. $$ (14) 而两个相机之间的位置关系可以通过空间转换矩阵
${M_{lr}}$ 表示：$$ \left[ {\begin{array}{*{20}{c}} {{x_r}} \\ {{y_r}} \\ {{{\textit{z}}_r}} \end{array}} \right] = {M_{lr}}\left[ {\begin{array}{*{20}{c}} x \\ y \\ {\textit{z}} \\ 1 \end{array}} \right]{\text{ = }}\left[ {\begin{array}{*{20}{c}} {{r_1}}&{{r_2}}&{{r_3}}&{{t_x}} \\ {{r_4}}&{{r_5}}&{{r_6}}&{{t_y}} \\ {{r_7}}&{{r_8}}&{{r_9}}&{{t_{\textit{z}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} x \\ y \\ {\textit{z}} \\ 1 \end{array}} \right] $$ (15) 对于空间中的一点P，成像于两个相机的点的对应关系可以用下面的等式表示：
$$ {\rho _l}\left[ {\begin{array}{*{20}{c}} {{X_r}} \\ {{Y_r}} \\ 1 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{f_r}{r_1}}&{{f_r}{r_2}}&{{f_r}{r_3}}&{{f_r}{t_x}} \\ {{f_r}{r_4}}&{{f_r}{r_5}}&{{f_r}{r_6}}&{{f_r}{t_y}} \\ {{r_7}}&{{r_8}}&{{r_9}}&{{t_{\textit{z}}}} \end{array}} \right]{\text{ = }}\left[ {\begin{array}{*{20}{c}} {{{{\textit{z}}{X_l}} / {{f_l}}}} \\ {{{{\textit{z}}{Y_l}} / {{f_l}}}} \\ {\textit{z}} \\ 1 \end{array}} \right] $$ (16) 于是，空间中的三维点可以通过上式解出：
$$\begin{split} &x = {\textit{z}}{X_l}/{f_l}\\ &y = {\textit{z}}{Y_l}/{f_l}\\ &{\textit{z}} = \dfrac{{{f_l}({f_r}{t_x} - {X_r}{t_{\textit{z}}})}}{{{X_r}({r_7}{X_l} + {r_8}{Y_l} + {f_l}{r_9}) - {f_r}({r_1}{X_l} + {r_2}{Y_l} + {f_l}{r_3})}} =\\ & \dfrac{{{f_l}({f_r}{t_y} - {Y_r}{t_{\textit{z}}})}}{{{Y_r}({r_7}{X_l} + {r_8}{Y_l} + {f_l}{r_9}) - {f_r}({r_4}{X_l} + {r_5}{Y_l} + {f_l}{r_6})}} \end{split}$$ (17) -
文中使用的器械经过精密加工，器械尖点在器械坐标系中的坐标
${P_{{\rm{ins}}}}$ $(0,214.14, - 13.77)$ 已知。如图8（b）所示，在器械坐标系中，${A_1}$ 为原点，${A_1}{A_2}$ 方向为x轴正方向，${A_1}{A_4}$ 方向为y轴正方向，z轴通过右手法则得出。经过立体匹配，${A_1}$ 的空间坐标${P_{{A_1}}}$ 则为世界坐标系与器械坐标系之间的平移矩阵$T$ ，旋转矩阵可通过公式（18）、（19）、（20）得出：$$ {n_x} = \dfrac{{{P_{{A_2}}} - {P_{{A_1}}}}}{{\left\| {{P_{{A_2}}} - {P_{{A_1}}}} \right\|}} $$ (18) $$ {n_y} = \dfrac{{{P_{{A_4}}} - {P_{{A_1}}}}}{{\left\| {{P_{{A_4}}} - {P_{{A_1}}}} \right\|}} $$ (19) $$ {n_{\textit{z}}} = {n_x} \times {n_y} $$ (20) 世界坐标系与器械坐标系的旋转矩阵
$R$ ，可以表示为：$ R = \left[ {{n_x},{n_y},{n_{\textit{z}}}} \right] $ 。则器械尖点在世界坐标系中的坐标
${P_{{\rm{cl}}}}$ 可以表示为：$$ {P}_{{\rm{cl}}}=R\cdot {P}_{{\rm{ins}}}+T $$ (21) -
为了能够更好的测试该系统，笔者进行了4组实验，测试器械识别算法的准确性和鲁棒性，并对器械进行了稳定性、定位精度、动态跟踪等测试。
首先，器械识别实验用于验证算法的准确性和鲁棒性，通过无规则摆放多组器械并且设置干扰进行识别得出结论；其次，稳定性测试用于验证系统的稳定性，这是通过以尖点为圆心旋转而重复测量尖端坐标来完成的；接下来，应用三维坐标仪来测量系统的静态定位精度，这是通过测量由尖端坐标跟随坐标仪移动的前后距离差值计算得出的；最后，通过不断改变器械空间位置，并实时计算两个标记之间的距离，与真值进行比较，用于评估该系统动态跟踪精度。文中中，平均绝对误差（MAE）表示程度测量值与真实值的偏差，而均方根误差（RMSE）用于评估定位精度。它们的表达式如下：
$$ RMSE = \sqrt {\dfrac{{\displaystyle\sum\limits_{i = 1}^n {{{\left( {{X_i} - {X_{ref}}} \right)}^2}} }}{n}} $$ (22) $$ MAE = \dfrac{1}{n}\displaystyle\sum\limits_{i = 1}^n {\left| {{X_i} - {X_{ref}}} \right|} $$ (23) -
文中使用Matlab相机标定箱对相机进行标定，为了提高标定的精确性，首先，分别对左右相机进行单目标定，然后再使用单目标定参数作为初始参数进行双目标定，最终标定出左右相机的内参矩阵
${A_{\rm{l}}}、{A_{\rm{r}}}$ ，畸变系数${D_{\rm{l}}}、{D_{\rm{r}}}$ ，右相机相对左相机的旋转矩阵${R_{{\rm{rl}}}}$ ，平移矩阵${T_{{\rm{rl}}}}$ ，结果如下：$$ {A_{\rm{l}}} = {A_{\rm{r}}} = \left[ {\begin{array}{*{20}{c}} {1751.783}&{0.3731}&{641.047} \\ 0&{1749.971}&{495.806} \\ 0&0&{1.0} \end{array}} \right] $$ (24) $$ {R_{{\rm{rl}}}} = \left[ {\begin{array}{*{20}{c}} {0.998}&{ - 0.016}&{ - 0.056} \\ {0.013}&{0.998}&{ - 0.043} \\ {0.057}&{0.042}&{0.997} \end{array}} \right] $$ (25) $$ {T_{{\rm{rl}}}} = \left[ {\begin{array}{*{20}{c}} {429.585}&{4.763}&{2.711} \end{array}} \right] $$ (26) $$ {D_{\rm{l}}} = \left[ {\begin{array}{*{20}{l}} { - 0.094} & {0.34} & { - 5.412e - 4} & { - 0.007} & { - 0.05} \end{array}} \right] $$ (27) $$ {D_{\rm{r}}} = \left[ {\begin{array}{*{20}{l}} { - 0.054} & {1.197} & { - 7.05e - 4} & { - 7.7e - 4} & { - 10.305} \end{array}} \right] $$ (28) -
首先，系统对表3中的任意两种器械识别，验证算法是否可以识别不同种类的器械。实验时将两个器械无规则摆放，但没有相互遮挡，共采集了100张图片。结果表明：只有4次没有识别成功，成功率为96%。未识别成功的主要原因是两个器械重合度过高，实际应用不会出现这种情况。识别成功情况如图9（a）所示。
图 9 (a)多器械识别; (b) 添加干扰
Figure 9. (a) Multiple instruments identification; (b) Add interference
为了验证器械识别算法的鲁棒性，还设计了干扰实验，对需识别器械添加干扰，该实验干扰设置为4个无规则摆放的被动标记球，如图9（b）所示。在存在干扰的情况下识别器械，共采集了100张图片。结果表明：只有5次没有成功识别，成功率为95% 。
图9（a）中红色连线组成的器械代表A类器械，蓝色代表C类器械。图9（b）中绿色连线组成的器械代表B类器械。
此次实验还测量了在有干扰和没有干扰标记的器械识别所需的时间，如表4 所示。结果表明：识别器械时间花费很少，并且不会随着器械的数量增加而显著增加。
表 4 多个器械识别所需时间（单位：毫秒）
Table 4. Time required for multiple instruments identification (Unit:ms)
Number of instruments Four interfering light sources No interfering light sources 1 4 3 2 4.2 3.5 3 4.5 4 -
光学跟踪系统工作时，系统必须具有足够的稳定性。如果器械尖端保持静止，由于捕获环境的变化和算法的稳定性，计算的尖端世界坐标可能会发生变化。
为了测试系统的稳定性，随机选取8个不同的位置，使器械尖端在每个位置都保持静止。对于每个位置， 采集100 组图像。然后分别计算 X、Y 和 Z 方向的RMSE，用于评估系统的稳定性，结果如图10所示。
图 10 尖点旋转测试
Figure 10. Tip rotation test
由图10可知，8个位置的X、Y和Z方向的平均RMSE分别为0.039625、0.024625、0.1284 mm。Z方向的RMSE误差大于其他两个方向，但是误差处于可接受范围，可能由于器械加工误差所致。
-
在定位精度验证实验中，将器械固定到三维坐标仪上，在距离光心1.2、1.4、2 m三个位置及1.6、1.8 m水平偏离光轴、1.6 m垂直偏离光轴另外三组位置进行测试。实验中器械将在三维坐标仪上沿一个方向移动5、10、15、20 mm，记录尖点移动的距离如表5所示。
表 5 定位精度测试（单位：毫米）
Table 5. Positioning accuracy test(Unit:mm)
Position 5 mm 10 mm 15 mm 20 mm 1 5.1087 10.3104 15.5100 20.5739 2 6.3136 10.6348 15.7608 20.8322 3 5.1060 10.0926 15.1268 20.3387 4 6.4282 12.4944 18.3975 23.2675 5 5.2923 10.4994 15.6380 20.8859 6 6.0823 11.6198 16.8552 22.4211 RMSE 0.6195 0.9252 1.2160 1.1768 MAE 0.5528 0.7435 0.9411 0.9718 由表5可知，5、10、15和20 mm的RMSE分别为0.6195、0.9252、1.2160、1.1768 mm，MAE分别为0.5528、0.7435、0.9411、0.9718 mm。表5中第4个位置误差较大，对应的是距离光心2 m处的实验结果，由第1、3、5组位置数据可知，系统在1.2 ~1.6 m区间范围内静态定位精度良好，但是随着距离的增加，误差会逐渐增大。实际使用时应合理设置工作距离，尽量设置在1.2 ~1.6 m的范围区间且在垂直方向不严重偏离光轴。由此，可知：系统在1.2~1.6 m的范围区间工作时，定位精度可达0.373 mm。
-
在动态跟踪实验中，笔者在系统视场内距离左右相机中心1.2~1.8 m的范围内随机采集100组图像，利用图4所示器械中的3组距离D16、D14、D24验证，结果如表6所示。
表 6 动态跟踪验证(单位：毫米)
Table 6. Dynamic tracking verification (Unit:mm)
Distance True value Max error RMSE MAE D16 70.5 0.4072 0.6184 0.5254 D14 70.7 0.1720 0.5205 0.4097 D24 100 1.0654 0.5765 0.4619 由表6可知，3组距离在1.2~1.8 m的范围内RMSE值较为接近。所有3种情况下的MAE
均小于0.53 mm，其中D24这组距离的最大误差超过1 mm，产生此结果的原因可能是当器械在系统的工作空间中移动时，标记球位于不同的视角，标记投影的变化可能会影响定位精度。在实际应用中，应避免出现严重偏离光轴的情况，这可以通过术前合理设置该系统中心光轴与器械的角度改善。 -
文中提出了一种基于双目视觉原理的低成本、高精度的光学跟踪系统。结合轮廓过滤算法，系统可以快速消除场景中非目标轮廓干扰。自主设计的器械识别算法根据器械特有几何特征对其正确识别并有序储存器械上标记球的像素坐标序列。实验结果表明：该系统具有较高的稳定性和准确性，可为缺少经验的外科医生的术前培训带来福音。但在实际应用过程中，该系统仍有改进的余地，首先，可以为器械选择一种对近红外光反射率比较低的材料，反射较高将影响标记中心提取的精确度；其次，可以为系统匹配处理速度更高的并行处理模块，这将大大提高系统的实时性。目前，随着对手术导航系统越来越精确的运行要求，光学跟踪系统的精度和稳定性需要不断提高。
A near-infrared binocular system for optical instrument tracking
-
摘要: 近红外光学跟踪系统以其高精度、便捷的特点迅速发展成为手术导航中的重要组成部分。基于双目视觉设计了一套高精度、低成本的光学跟踪系统，该系统使用安装在器械上的被动标记球（NDI Passive Sphere）作为标记，并在双目相机镜头前添加近红外滤光片，以消除环境光的干扰。首先，采用轮廓过滤算法，用于标记轮廓的提取，使用最小二乘椭圆拟合算法获取标记投影中心的像素坐标；其次，设计了器械识别算法，使每个器械可以单独有序储存标记球的中心像素坐标，用于多个器械的分辨；最后，通过左右视图对应标记中心匹配，重建其空间坐标，进而推导出器械尖端在世界坐标系中的坐标。实验使用该系统跟踪器械，利用器械的无规则摆放实验验证了器械识别算法的准确性和鲁棒性，准确率达95%，平均识别时间仅需4 ms。并对系统进行了稳定性、静态定位精度、动态跟踪等测试。结果表明：稳定性误差可达0.13 mm，静态定位精度可达0.373 mm，所提出的近红外光学跟踪系统具有较高的精度和稳定性，可以满足手术导航的要求。Abstract: The near-infrared optical tracking system has rapidly developed into an important part of surgical navigation because of its high precision and convenience. A high-precision and low-cost optical tracking system was designed based on binocular vision. The system used the passive marker ball (NDI passive sphere) installed on the instrument as the marker, and added a near-infrared filter in front of the binocular camera lens to eliminate the interference of ambient light. Firstly, the contour filtering algorithm was used to extract the marker contour, and the least square ellipse fitting algorithm was used to obtain the pixel coordinates of the marker projection center; Secondly, an instrument recognition algorithm was designed so that each instrument could store the central pixel coordinates of the marker ball separately and orderly for the resolution of multiple instruments; Finally, by matching the mark centers corresponding to the left and right views, the spatial coordinates were reconstructed, and then the coordinates of the instrument tip in the world coordinate system were derived. The system was used to track the instruments in the experiment, and the irregular placement of instruments was used to verify the accuracy and robustness of the instrument recognition algorithm. The accuracy rate was 95%, and the average recognition time was only 4 ms. The stability, static positioning accuracy and dynamic tracking of the system were tested. The results showed that the stability error could reach 0.13 mm, the static positioning accuracy could reach 0.373 mm. The proposed near-infrared optical tracking system has high accuracy and stability, and can meet the requirements of surgical navigation.
-
Key words:
- instrument recognition /
- surgical navigation /
- optical positioning /
- binocular vision
-
图 5 (a)轮廓错误提取; (b)轮廓正确提取
Figure 5. (a) Contour error extraction; (b) Correct contour extraction
图 8 (a) 相互位置; (b) 器械标记顺序
Figure 8. (a) Mutual position; (b) Instrument marking sequence
图 9 (a)多器械识别; (b) 添加干扰
Figure 9. (a) Multiple instruments identification; (b) Add interference
表 1 tr阈值选取
Table 1. tr threshold selection
Threshold Number 0.400 8 0.300 7 0.226 7 0.225 6 0.224 6 0.200 6 0.100 6 
下载: 导出CSV
表 3 器械种类
Table 3. Types of instruments
${{{L_1}} \mathord{\left/ {\vphantom {{{L_1}} {{L_2}}}} \right. } {{L_2}}}$ ${{{L_3}} \mathord{\left/ {\vphantom {{{L_3}} {{L_4}}}} \right. } {{L_4}}}$ Types of devices > 1 > 1 A < 1 < 1 B > 1 < 1 C
下载: 导出CSV
表 4 多个器械识别所需时间（单位：毫秒）
Table 4. Time required for multiple instruments identification (Unit:ms)
Number of instruments Four interfering light sources No interfering light sources 1 4 3 2 4.2 3.5 3 4.5 4
下载: 导出CSV
表 5 定位精度测试（单位：毫米）
Table 5. Positioning accuracy test(Unit:mm)
Position 5 mm 10 mm 15 mm 20 mm 1 5.1087 10.3104 15.5100 20.5739 2 6.3136 10.6348 15.7608 20.8322 3 5.1060 10.0926 15.1268 20.3387 4 6.4282 12.4944 18.3975 23.2675 5 5.2923 10.4994 15.6380 20.8859 6 6.0823 11.6198 16.8552 22.4211 RMSE 0.6195 0.9252 1.2160 1.1768 MAE 0.5528 0.7435 0.9411 0.9718
下载: 导出CSV
表 6 动态跟踪验证(单位：毫米)
Table 6. Dynamic tracking verification (Unit:mm)
Distance True value Max error RMSE MAE D16 70.5 0.4072 0.6184 0.5254 D14 70.7 0.1720 0.5205 0.4097 D24 100 1.0654 0.5765 0.4619
下载: 导出CSV
-
[1] Ma S, Zhao Z. A new method of surgical tracking system based on fiducial marker[C]//Annual Conference on Medical Image Understanding and Analysis, 2017: 886-896. [2] An H W, Hou D, Zheng R, et al. A near-infrared peptide probe with tumor-specific excretion-retarded effect for image-guided surgery of renal cell carcinoma [J]. ACS Nano, 2020, 14(1): 927-936. doi: 10.1021/acsnano.9b08209 [3] Jolesz Ferenc A. Intraoperative Imaging and Image-guided Therapy Preoperative [M]. US:Springer, 2014, : 677-683. [4] Zhou Z, Wu B, Duan J, et al. Optical surgical instrument tracking system based on the principle of stereo vision [J]. J Biomed Opt, 2017, 22(6): 65005. doi: 10.1117/1.JBO.22.6.065005 [5] Zhou P, Liu Y, Wang Y. Multiple infrared markers based real-time stereo vision positioning system for surgical navigation[C]//2009 IEEE Instrumentation and Measurement Technology Conference, 2009: 692-696. [6] Steinicke F, Jansen C P, Hinrichs K H, et al. Generating optimized marker-based rigid bodies for optical tracking systems[C]//VISAPP, 2007,2: 387-395. [7] Yan G, Li J, Huang Y, et al. Ghost marker detection and elimination in marker‐based optical tracking systems for real‐time tracking in stereotactic body radiotherapy [J]. Medical Physics, 2014, 41(10): 101713. doi: 10.1118/1.4896126 [8] Lin Q, Yang R, Zhang Z, et al. Robust stereo-match algorithm for infrared markers in image-guided optical tracking system [J]. IEEE Access, 2018, 6: 52421-52433. doi: 10.1109/ACCESS.2018.2869433 [9] Wu H, Lin Q, Yang R, et al. An accurate recognition of infrared retro-reflective markers in surgical navigation [J]. Journal of Medical Systems, 2019, 43(6): 1-13. [10] Zhang M, Wu B, Ye C, et al. Multiple instruments motion trajectory tracking in optical surgical navigation [J]. Optics Express, 2019, 27(11): 15827-15845. doi: 10.1364/OE.27.015827 [11] Zhang Z. A flexible new technique for camera calibration [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2000, 22(11): 1330-1334. doi: 10.1109/34.888718 [12] Canny J. A computational approach to edge detection [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1986(6): 679-698. [13] Fitzgibbon A W. Direct least square fitting of ellipses[C]//Pattern Analysis and Machine Intelligence, 1999: 253-257. [14] 李云雷. 基于近景工业摄影的光学三维形貌测量关键技术研究[D]. 上海: 上海大学, 2020. Li Yunlei. Research on key technologies of optical three-dimensional shape measurement based on close-range industry photogrammetry[D]. Shanghai: Shanghai University, 2020. (in Chinese) -
点击查看大图
计量
- 文章访问数: 383
- HTML全文浏览量: 116
- PDF下载量: 97
- 被引次数: 0
