CN105434046A - 基于消除红外滤光片折射效应的手术导航仪定位方法 - Google Patents

Abstract

本发明提供了一种基于消除红外滤光片折射效应的手术导航仪定位方法，包括步骤1：标定无红外滤光片时左右两侧相机的内参数矩阵、镜头畸变参数以及两个相机之间的相对位置关系；步骤2：测量红外滤光片的厚度、折射率；步骤3：对左右相机采集的图像分别进行镜头畸变校正；步骤4：提取所述镜头畸变校正后的图像中每个红外标志点的中心，得到红外标志点对应的图像点对；步骤5：根据相机成像模型、左右两侧相机的相对位置关系以及红外标志点对应的图像点对，计算红外标志点的空间位置。本发明消除了红外滤光片的折射带来的系统定位误差，提高了手术导航仪的定位精度。

Description

基于消除红外滤光片折射效应的手术导航仪定位方法

技术领域

本发明涉及医学手术导航，具体地，涉及基于消除红外滤光片折射效应的手术导航仪定位方法。

背景技术

光学导航是目前手术导航中的主流方法，具有精度高、应用灵活方便、无侵入伤害等优点，因而在神经外科、脊柱外科(骨科)、耳鼻喉科、整形外科等领域有着广泛的应用。目前光学导航产品几乎被国外垄断，国内应用的产品有加拿大NDI公司的Polaris系统，美国Medtronic公司的Stealthstation系统，德国Stryker公司的StrykeNavigationSystem系统和BrainLab公司的3DVectorVision导航系统等，为了打破垄断，降低医疗成本，国内一些高校和公司竞相开展了对光学手术导航系统的研究，但距离国外得最高水平还有很长的一段距离要走。

光学手术导航系统其核心为基于近红外光的立体视觉定位仪，其由一对灰度相机和镜头前加装的红外滤光片构成。红外滤光片能够剔除环境光的干扰，快速准确地分割出手术器械上的标志点，使用红外滤光片提取红外标志点的位置具有设计简单、成本低廉、易于实现等优点，因此获得了广泛的应用。根据手术器械标志点是否发光，立体视觉定位仪分为主动式和被动式两种：主动式立体视觉定位仪是在手术器械上安装数个红外发光二极管；被动式立体视觉定位仪是在手术器械上安装数个红外反光小球，在相机周围安装红外光源，然后使用相机透过红外滤光片采集红外反光小球的位置。

在立体视觉定位仪定位过程中，其两个相机分别透过红外滤光片对红外标志点进行成像，在两视图中分别提取到红外标志点的中心后，使用传统的三角测量法计算红外标志点的空间位置，在此过程中并没有考虑红外滤光片对定位精度的影响。通过分析发现，当相机通过红外滤光片对目标进行成像时，光线从目标点出发穿过红外滤光片，由于光的折射作用，光线在红外滤光片的表面发生两次折射，改变光线的传播方向。由于红外滤光片折射的影响，传统的针孔相机成像模型以及相关的理论方法将不再适用，此时使用传统的三角测量法对空间点进行定位，将引进系统定位误差。

目前，由于红外滤光片的折射对成像影响的复杂性，国内所研发的光学手术导航仪均没有考虑红外滤光片的折射对空间定位的影响。由于光学手术导航仪对精度的要求非常高，可谓精度越高越好，因此有必要考虑红外滤光片的折射效应对空间定位的影响，消除红外滤光片折射带来的定位误差对于提高手术导航仪的空间定位精度具有重要的实用价值。

发明内容

针对现有技术中的缺陷，本发明的目的是提供基于消除红外滤光片折射效应的手术导航仪定位方法。

根据本发明提供的基于消除红外滤光片折射效应的手术导航仪定位方法，包括如下步骤：

步骤1：标定在无红外滤光片时左右两侧相机的内参数矩阵、镜头畸变参数以及左右两侧相机之间的相对位置关系；

步骤2：测量红外滤光片的厚度、折射率，然后将红外滤光片分别安装在左右两侧相机的镜头前；

步骤3：使用光学导航仪采集手术器械图像，利用相机内参矩阵和镜头畸变参数对手术器械图像进行镜头畸变校正，得到镜头畸变校正后的图像；

步骤4：将镜头畸变校正后的图像进行二值处理，提取所述镜头畸变校正后的图像中每个红外标志点的成像区域，得到每个红外标志点在图像中对应的中心点图像坐标，然后左右两幅图像匹配得到红外标志点对应的图像点对；

步骤5：根据相机的内参矩阵、红外滤光片的厚度和折射率得到相机透过红外滤光片成像的相机成像模型，然后根据该相机成像模型、左右两侧相机之间的相对位置关系以及红外标志点对应的图像点对，计算红外标志点的空间位置。

优选地，所述步骤1包括：

步骤1.1：使用单摄像机标定算法对左右两侧相机分别进行摄像机标定，得到摄像机内参矩阵和镜头畸变参数，得到相机内参矩阵(A_l,A_r)和镜头畸变参数(S_l,S_r)，计算公式如下：

$A_l=\left[\begin{array}{ccc}{f}_{lx}& 0& c_{lx}\\ 0& f_{ly}& c_{ly}\\ 0& 0& 1\end{array}\right];$

$A_r=\left[\begin{array}{ccc}{f}_{rx}& 0& c_{rx}\\ 0& f_{ry}& c_{ry}\\ 0& 0& 1\end{array}\right];$

S_l＝[k_l1k_l2p_l1p_l2]；

S_r＝[k_r1k_r2p_r1p_r2]；

式中，其中下标l、r分别表示左右相机对应的相关参数，A_l表示左侧相机内参矩阵，A_r表示右侧相机内参矩阵，S_l表示左侧相机畸变参数，S_r表示右侧相机畸变参数，f_lx表示左侧相机在像素量纲下关于x轴方向上的焦距，f_rx表示右侧相机在像素量纲下关于x轴方向上的焦距，f_ly表示左侧相机在像素量纲下关于y轴方向上的焦距，f_ry表示右侧相机在像素量纲下关于y轴方向上的焦距，c_lx表示左侧相机主点关于x轴的值，c_ly表示左侧相机主点关于y轴的值，c_rx表示右侧相机主点关于x轴的值，c_ry表示右侧相机主点关于y轴的值，k_l1表示左侧相机镜头2阶径向畸变参数，k_l2表示左侧相机镜头4阶径向畸变参数，k_r1表示右侧相机镜头2阶径向畸变参数，k_r2表示右侧相机镜头4阶径向畸变参数，p_l1表示左侧相机镜头离心畸变参数，p_l2表示左侧相机镜头离心畸变参数，p_r1表示右侧相机镜头离心畸变参数，p_r2表示右侧相机镜头离心畸变参数；

步骤1.2：使用双目相机标定算法得到左右两侧相机之间的相对旋转矩阵和平移向量(R,t)，R表示右侧相机坐标系相对于左侧相机坐标系的旋转矩阵，t表示右侧相机坐标系相对于左侧相机坐标系的平移向量。

优选地，所述步骤2包括：

步骤2.1：使用游标卡尺或螺旋测微器测量红外滤光片的厚度，多次测量后取均值，记所述左右两侧相机所用的红外滤光片的厚度分别为d_l、d_r；

步骤2.2：使用折射率测量仪测量红外滤光片的折射率为n；

步骤2.3：将红外滤光片安装在左右两侧相机的镜头前。

优选地，所述步骤3包括：

步骤3.1：使用光学导航仪采集手术器械图像，即安装红外滤光片后相机采集到的包含红外标志点的图像，其中，将左侧相机采集到的图像记为I_l；将右侧相机采集到的图像记为I_r；

步骤3.2：利用相机内参矩阵(A_l,A_r)和镜头畸变参数(S_l,S_r)对图像I_l、I_r进行镜头畸变矫正，在此只考率径向畸变和离心畸变，对左右两幅图像分别采用如下畸变校正方法进行校正，校正公式如下：

$\left[\begin{array}{c}U\\ V\end{array}\right]=\left[\begin{array}{c}{U}_d\\ V_d\end{array}\right]+\left[\begin{array}{c}{\mathrm{\δ}}_{Udr}\\ {\mathrm{\δ}}_{Vdr}\end{array}\right]+\left[\begin{array}{c}{\mathrm{\δ}}_{Udd}\\ {\mathrm{\δ}}_{Vdd}\end{array}\right];$

式中，δ_Udr表示观测到的图像点(U_d，V_d)在x方向上的径向畸变的畸变量，δ_Vdr分别表示观测到的图像点(U_d，V_d)在y方向上的径向畸变的畸变量，δ_Udd表示观测到的图像点(U_d，V_d)在x方向上的离心畸变的畸变量，δ_Vdd分别表示观测到的图像点(U_d，V_d)在y方向上的离心畸变的畸变量，U_d为观测到的图像点在x方向上的图像坐标，V_d为观测到的图像点在y方向上的图像坐标，U表示观测到的图像坐标(U_d，V_d)在经过畸变校正后在x方向上的图像坐标，V分别表示观测到的图像坐标(U_d，V_d)在经过畸变校正后在y方向上的图像坐标，其中：

$\left[\begin{array}{c}{\mathrm{\δ}}_{Udr}\\ {\mathrm{\δ}}_{Vdr}\end{array}\right]=(k_1{r_d}^2+k_2{r_d}^4)\left[\begin{array}{c}{U}_d\\ V_d\end{array}\right],$

r_d ²＝U_d ²+V_d ²

式中，r_d是从图像主点至畸变图像点(U_d，V_d)的距离，k₁表示二阶径向畸变参数，k₂表示四阶径向畸变参数。

优选地，所述步骤4包括：

步骤4.1：对镜头畸变校正后的图像I_l1、I_r1进行二值处理，得到左侧相机的二值图像I_l2和右侧相机的二值图像I_r2；

步骤4.2：利用区域提取算法对左侧相机的二值图像I_l2和右侧相机的二值图像I_r2进行处理，得到左侧相机和右侧相机采集到的图像中每个红外标志点的成像区域；

步骤4.3：使用重心法或者任何区域中心提取算法，计算红外标志点对应区域的中心，并把该中心作为红外标志点的中心在图像中对应的位置；

步骤4.4：使用对极约束法快速匹配左侧相机和右侧相机采集到的图像中的红外标志点，得到红外标志点对应的图像点对(x_l,x_r)，其中x_l表示红外标志点的中心在左图像中对应的图像点，x_r表示同一个红外标志点的中心在右图像中对应的图像点。

优选地，所述步骤5包括：

步骤5.1：根据相机的内参矩阵、红外滤光片的厚度和折射率得到相机透过红外滤光片的相机成像模型；

当以左侧相机坐标系作为世界坐标系时，即令相机坐标系的原点位于相机中心，则对于每一个图像点对(x_l,x_r)，其中 $x_l=\left(\begin{array}{c}{u}_l\\ v_l\end{array}\right)$ 表示左图像点的图像坐标， $x_r=\left(\begin{array}{c}{u}_r\\ v_r\end{array}\right)$ 表示右图像点的图像坐标，u_l表示左图像中的图像点在x方向上的坐标值，v_l表示左图像中的图像点在y方向上的坐标值，u_r表示右图像中的图像点在x方向上的坐标值，v_r表示左图像中的图像点在y方向上的坐标值；

左相机投影矩阵P_l为P_l＝A_l[It_dl]，右相机投影矩阵P_r为P_r＝A_r[Rt-t_dr]，其中：左侧相机虚拟视点在左相机坐标系中的位置为t_dl＝(00d_Ol)^T；右侧相机虚拟视点在右相机坐标系中的位置为t_dr＝(00d_Or)^T；d_Ol表示左相机对应于图像点对(x_l,x_r)的相机虚拟视点在左相机光轴上的平移距离；d_Or表示右两侧相机对应于图像点对(x_l,x_r)的相机虚拟视点在右相机光轴上的平移距离；d_Ol和d_Or可以使用如下公式进行计算：

$d_{Ol}=d_l(1-\frac{\sqrt{1-{\sin }^2{\mathrm{\θ}}_l}}{\sqrt{n^2-{\sin }^2{\mathrm{\θ}}_l}}),d_{Or}=d_r(1-\frac{\sqrt{1-{\sin }^2{\mathrm{\θ}}_r}}{\sqrt{n^2-{\sin }^2{\mathrm{\θ}}_r}}),$

式中，d_l表示左相机镜头前红外滤光片的厚度；d_r表示右相机镜头前红外滤光片的厚度，通过游标卡尺或是千分尺进行测量；n表示红外滤光片的折射率；其中，

${\mathrm{sin\θ}}_l=\frac{r_l}{\sqrt{{r_l}^2+{f_l}^2}};$

f_l＝(f_lx+f_ly)/2；

$r_l=\sqrt{{u_l}^2+v_l^2};$

${\mathrm{sin\θ}}_r=\frac{r_r}{\sqrt{{r_r}^2+f_r^2}};$

f_r＝(f_rx+f_ry)/2；

$r_r=\sqrt{{u_r}^2+v_r^2};$

式中：θ_l表示左相机中图像点对应的光线相对于红外滤光片的入射角，θ_r表示右相机中图像点对应的光线相对于红外滤光片的入射角，r_l表示左图像点到左图像主点的距离，r_r表示右图像点到右图像主点的距离，f_l表示左相机的焦距，f_r表示右相机的焦距；

步骤5.2：根据相机成像模型、左侧和右侧相机的相对位置关系以及红外标志点对应的图像点对计算红外标志点的空间位置；

对于每一对对应的图像点对(x_l,x_r)，可以使用任何传统的三角测量法进行计算，此处以传统的线性三角测量法为例进行计算，具体计算如下：

$s_l{u}_l=p_{l1}^T\hat{X},s_l{v}_l=p_{l2}^T\hat{X},s_l=p_{l3}^T\hat{X};$

$s_r{u}_r=p_{r1}^T\hat{X},s_r{v}_r=p_{r2}^T\hat{X},s_r=p_{r3}^T\hat{X};$

其中： $\begin{array}{c}{u}_l{p}_{l3}^T\hat{X}-p_{l1}^T\hat{X}=0\\ v_l{p}_{l3}^T\hat{X}-p_{l2}^T\hat{X}=0\\ u_r{p}_{r3}^T\hat{X}-p_{r1}^T\hat{X}=0\\ v_r{p}_{r3}^T\hat{X}-p_{r2}^T\hat{X}=0\end{array};$

则改写为 $A\hat{X}=0$ 时， $A=\left[\begin{array}{c}{u}_l{p}_{l3}^T-p_{l1}^T\\ v_l{p}_{l3}^T-p_{l2}^T\\ u_r{p}_{r3}^T-p_{r1}^T\\ v_r{p}_{r3}^T-p_{r2}^T\end{array}\right];$

式中：表示左侧相机投影矩阵的第一行，表示待求空间点的坐标，表示左侧相机投影矩阵的第二行，表示左侧相机投影矩阵的第三行，s_l表示比例因子，s_r表示比例因子，表示右侧相机投影矩阵的第一行，表示右侧相机投影矩阵的第二行，表示右侧相机投影矩阵的第一行。

优选地，所述步骤5中所应用的相机模型能够逐一对应图像的像素变化，即每一像素点相当于相机平移到了新的位置，相机成像模型也变为逐像素点相机中心变化的针孔相机成像模型。

与现有技术相比，本发明具有如下的有益效果：

1、本发明提出的基于消除红外滤光片折射效应的手术导航仪定位方法，简单可靠容易实现，消除了红外滤光片的折射带来的系统定位误差。

2、本发明提出的基于消除红外滤光片折射效应的手术导航仪定位方法不仅能够应用在光学手术导航仪的精确定位中，还应用于相机通过透明玻璃或其他透明平板对空间目标进行高精度定位测量以及三维重建等工作，具有重要的实用价。

附图说明

通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：

图1为本发明的手术导航仪的成像示意图；

图2为本发明的手术导航仪的成像原理示意图；

图3为本发明的手术导航仪的空间定位原理示意图。

图中：

1-目标；

2-红外滤光片；

3-左侧相机；

4-右侧相机。

具体实施方式

下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进。这些都属于本发明的保护范围。

根据本发明提供的基于消除红外滤光片折射效应的手术导航仪定位方法，利用单相机标定方法和立体相机标定方法标定摄像机的内参数、畸变参数以及左右两侧相机的相对位置关系(参考文献：J.Y.Bouguet.Cameracalibrationtoolboxformatlab.http://www.vision.caltech.edu/bouguetj/calib_doc/index.html.)；然后测量红外滤光片的厚度和折射率；再次通过红外滤光片采集红外标志点的图像；最后通过计算得到红外滤光片的折射率。

使用单相机标定算法和立体相机标定方法进行相机标定，得到相机内参矩阵和畸变参数，使用单相机标定方法和立体相机标定方法进行相机标定，得到相机内参矩阵(A_l,A_r)和镜头畸变参数(S_l,S_r)的计算公式如下：使用单摄像机标定算法(参考文献：Z.Zhang.Aflexiblenewtechniqueforcameracalibration.IEEETransactionsonPatternAnalysisandMachineIntelligence,22(11):1330-1334,2000.)对左右两侧相机分别进行摄像机标定，得到摄像机内参矩阵和镜头畸变参数，其中左侧相机内参矩阵为 $A_l=\left[\begin{array}{ccc}{f}_{lx}& 0& c_{lx}\\ 0& f_{ly}& c_{ly}\\ 0& 0& 1\end{array}\right];$ 右侧相机内参矩阵为 $A_r=\left[\begin{array}{ccc}{f}_{rx}& 0& c_{rx}\\ 0& f_{ry}& c_{ry}\\ 0& 0& 1\end{array}\right];$ 左侧相机畸变参数为S_l＝[k_l1k_l2p_l1p_l2]；右侧相机畸变参数为S_r＝[k_r1k_r2p_r1p_r2]；其中：下标l和r分别表示左右两侧相机的相应参数，f_lx表示左侧相机在像素量纲下关于x轴方向上的焦距，f_rx表示右侧相机在像素量纲下关于x轴方向上的焦距，f_ly表示左侧相机在像素量纲下关于y轴方向上的焦距，f_ry表示右侧相机在像素量纲下关于y轴方向上的焦距；c_lx表示左侧相机主点关于x轴的值，c_ly表示左侧相机主点关于y轴的值，c_rx表示右侧相机主点关于x轴的值，c_ry表示右侧相机主点关于y轴的值；k_l1表示左侧相机镜头2阶径向畸变参数，k_l2表示左侧相机镜头4阶径向畸变参数，k_r1表示右侧相机镜头2阶径向畸变参数，k_r2表示右侧相机镜头4阶径向畸变参数，p_l1表示左侧相机镜头离心畸变参数，p_l2表示左侧相机镜头离心畸变参数，p_r1表示右侧相机镜头离心畸变参数，p_r2表示右侧相机镜头离心畸变参数；A_l表示左侧摄像机矩阵，A_r表示右侧摄像机矩阵，S_l表示左侧相机镜头畸变参数，S_r表示右侧相机镜头畸变参数；使用双目相机标定算法(参考文献：J.-Y.Bouguet.Cameracalibrationtoolboxformatlab.Available:http://www.vision.caltech.edu/bouguetj/calibdoc/index.html)得到左右相机之间的相对旋转矩阵和平移向量(R,t)，R表示右相机坐标系相对于左相机坐标系的旋转矩阵，t表示右相机坐标系相对于左相机坐标系的平移向量；

使用游标卡尺或螺旋测微器测量红外滤光片的厚度d，然后使用折射率测量仪测量红外滤光片的折射率为n。

对加装红外滤光片后采集的红外标志点的图像I_l、I_r，进行畸变矫正，在此只考率径向畸变和离心畸变(参考文献：J.Weng.P.Cohen,andM.Herniou,CameraCalibrationwithDistortionModelsandAccuracyEvaluation,IEEETrans.PatternAnalysisandMachineIntelligence,1992,14(10):965-980.)。

具体地，

$\left[\begin{array}{c}U\\ V\end{array}\right]=\left[\begin{array}{c}{U}_d\\ V_d\end{array}\right]+\left[\begin{array}{c}{\mathrm{\δ}}_{Udr}\\ {\mathrm{\δ}}_{Vdr}\end{array}\right]+\left[\begin{array}{c}{\mathrm{\δ}}_{Udd}\\ {\mathrm{\δ}}_{Vdd}\end{array}\right]$

式中，δ_Udr和δ_Vdr为径向畸变，δ_Udd和δ_Vdd为离心畸变，(U_d，V_d)为畸变图像点，(U，V)为校正后图像点，

$\left[\begin{array}{c}{\mathrm{\δ}}_{Udr}\\ {\mathrm{\δ}}_{Vdr}\end{array}\right]=(k_1{r_d}^2+k_2{r_d}^4)\left[\begin{array}{c}{U}_d\\ V_d\end{array}\right],$ r_d ²＝U_d ²+V_d ²

式中，r_d是从畸变中心至畸变图像点(U_d，V_d)的距离，

$\left[\begin{array}{c}{\mathrm{\δ}}_{Udd}\\ {\mathrm{\δ}}_{Vdd}\end{array}\right]=\left[\begin{array}{c}{p}_1(3{U_d}^2+{V_d}^2)+2p_2\·U_d\·V_d\\ p_2(3{V_d}^2+{U_d}^2)+2p_1\·U_d\·V_d\end{array}\right]$

畸变校正后的图像，记为I_l1、I_r1。

对图像I_l1、I_r1进行二值处理，得到二值图像I_l2、I_r2，在图像I_l2、I_r2上应用区域提取算法，在左右两幅图像中提取每个红外标志点的成像区域，使用重心法或者使用任何区域中心提取算法计算每个成像区域的中心作为红外标志点的中心在图像中的成像点，使用对极约束法(参考文献：R.HartleyandA.Zisserman,Multipleviewgeometryincomputervision.CambridgeUniversityPress,2000,vol.2.)快速匹配左右两幅图像中的红外标志点，得到红外标志点对应的图像点对(x_l,x_r)。

具体地，如图3所示，红外滤光片折射改变光线传播方向，此处使用逐像点相机平移的方法构建相机成像模型，然后根据该成像模型使用线性三角测量法(此处可以使用任何传统的三角测量法)计算空间点X的位置。

构建逐像素点相机成像模型，具体实施如下：

假设以左相机坐标系作为世界坐标系，相机坐标系的原点位于相机中心，对于每一个图像点对(x_l,x_r)，对应于对应图像点对的投影矩阵分别为：

P_l＝A_l[It_dl]

P_r＝A_r[Rt-t_dr]

其中t_dl＝(00d_Ol)^T和t_dr＝(00d_Or)^T是两个虚拟相机对应与真实相机的平移向量，红外滤光片的折射对成像的影响相当于相机沿着光轴发生了平移，把折射成像问题转变为逐像素点的针孔成像问题。根据光线穿过平面平行板的相关理论，可得：

$d_{Ol}=d_l(1-\frac{\sqrt{1-{\sin }^2{\mathrm{\θ}}_l}}{\sqrt{n^2-{\sin }^2{\mathrm{\θ}}_l}})$

$d_{Or}=d_r(1-\frac{\sqrt{1-{\sin }^2{\mathrm{\θ}}_r}}{\sqrt{n^2-{\sin }^2{\mathrm{\θ}}_r}})$

其中，

${\mathrm{sin\θ}}_l=\frac{r_l}{\sqrt{r_l^2+f_l^2}},$ f_l＝(f_lx+f_ly)/2， $r_l=\sqrt{{u_l}^2+v_l^2},x_l=\left(\begin{array}{c}{u}_l\\ v_l\end{array}\right)$

${\mathrm{sin\θ}}_r=\frac{r_r}{\sqrt{{r_r}^2+f_r^2}},$ f_r＝(f_rx+f_ry)/2， $r_r=\sqrt{{u_r}^2+v_r^2},x_r=\left(\begin{array}{c}{u}_r\\ v_r\end{array}\right)$

式中：

对于每一像素点相当于相机平移到了新的位置，成像模型也变为逐像素点虚拟相机中心变化的针孔相机成像模型，计算得到每一像素点对的摄像机矩阵P_l和P_r，对于每一对对应的图像点对(x_l,x_r)，可以使用任何传统的三角测量法进行计算，此处以传统的线性三角测量法为例进行计算，具体计算如下：

上述两视图中，空间三维点以及图像点可以写作如下的形式，

$s_l{u}_l=p_{l1}^T\hat{X},s_l{v}_l=p_{l2}^T\hat{X},s_l=p_{l3}^T\hat{X}$

$s_r{u}_r=p_{r1}^T\hat{X},s_r{v}_r=p_{r2}^T\hat{X},s_r=p_{r3}^T\hat{X}$

其中和分别表示对应相机投影矩阵的第i行，s_l和s_r为比例因子。

消去s_l和s_r，从每一图像点取两个方程组成一个具有四个方程的线性方程组：

$\begin{array}{c}{u}_l{p}_{l3}^T\hat{X}-p_{l1}^T\hat{X}=0\\ v_l{p}_{l3}^T\hat{X}-p_{l2}^T\hat{X}=0\\ u_r{p}_{r3}^T\hat{X}-p_{r1}^T\hat{X}=0\\ v_r{p}_{r3}^T\hat{X}-p_{r2}^T\hat{X}=0\end{array}$

这个方程组可以写为的形式，它是关于的线性方程组

$A=\left[\begin{array}{c}{u}_l{p}_{l3}^T-p_{l1}^T\\ v_l{p}_{l3}^T-p_{l2}^T\\ u_r{p}_{r3}^T-p_{r1}^T\\ v_r{p}_{r3}^T-p_{r2}^T\end{array}\right]$

使用非齐次方法求解上述线性方程组即可以得到空间坐标的值，即为空间点X的空间点测量值。

以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变形或修改，这并不影响本发明的实质内容。

Claims (7)

1.一种基于消除红外滤光片折射效应的手术导航仪定位方法，其特征在于，包括如下步骤：

步骤1：标定在无红外滤光片时左右两侧相机的内参数矩阵、镜头畸变参数以及左右两侧相机之间的相对位置关系；

步骤2：测量红外滤光片的厚度、折射率，然后将红外滤光片分别安装在左右两侧相机的镜头前；

步骤3：使用光学导航仪采集手术器械图像，利用相机内参矩阵和镜头畸变参数对手术器械图像进行镜头畸变校正，得到镜头畸变校正后的图像；

步骤4：将镜头畸变校正后的图像进行二值处理，提取所述镜头畸变校正后的图像中每个红外标志点的成像区域，得到每个红外标志点在图像中对应的中心点图像坐标，然后左右两幅图像匹配得到红外标志点对应的图像点对；

步骤5：根据相机的内参矩阵、红外滤光片的厚度和折射率得到相机透过红外滤光片成像的相机成像模型，然后根据该相机成像模型、左右两侧相机之间的相对位置关系以及红外标志点对应的图像点对，计算红外标志点的空间位置。

2.根据权利要求1所述的基于消除红外滤光片折射效应的手术导航仪定位方法，其特征在于，所述步骤1包括：

步骤1.1：使用单摄像机标定算法对左右两侧相机分别进行摄像机标定，得到摄像机内参矩阵和镜头畸变参数，得到相机内参矩阵(A_l,A_r)和镜头畸变参数(S_l,S_r)，计算公式如下：

$A_l=\left[\begin{array}{ccc}{f}_{lx}& 0& c_{lx}\\ 0& f_{ly}& c_{ly}\\ 0& 0& 1\end{array}\right];$

$A_r=\left[\begin{array}{ccc}{f}_{rx}& 0& c_{rx}\\ 0& f_{ry}& c_{ry}\\ 0& 0& 1\end{array}\right];$

S_l＝[k_l1k_l2p_l1p_l2]；

S_r＝[k_r1k_r2p_r1p_r2]；

式中，其中下标l、r分别表示左右相机对应的相关参数，A_l表示左侧相机内参矩阵，A_r表示右侧相机内参矩阵，S_l表示左侧相机畸变参数，S_r表示右侧相机畸变参数，f_lx表示左侧相机在像素量纲下关于x轴方向上的焦距，f_rx表示右侧相机在像素量纲下关于x轴方向上的焦距，f_ly表示左侧相机在像素量纲下关于y轴方向上的焦距，f_ry表示右侧相机在像素量纲下关于y轴方向上的焦距，c_lx表示左侧相机主点关于x轴的值，c_ly表示左侧相机主点关于y轴的值，c_rx表示右侧相机主点关于x轴的值，c_ry表示右侧相机主点关于y轴的值，k_l1表示左侧相机镜头2阶径向畸变参数，k_l2表示左侧相机镜头4阶径向畸变参数，k_r1表示右侧相机镜头2阶径向畸变参数，k_r2表示右侧相机镜头4阶径向畸变参数，p_l1表示左侧相机镜头离心畸变参数，p_l2表示左侧相机镜头离心畸变参数，p_r1表示右侧相机镜头离心畸变参数，p_r2表示右侧相机镜头离心畸变参数；

步骤1.2：使用双目相机标定算法得到左右两侧相机之间的相对旋转矩阵和平移向量(R,t)，R表示右侧相机坐标系相对于左侧相机坐标系的旋转矩阵，t表示右侧相机坐标系相对于左侧相机坐标系的平移向量。

3.根据权利要求2所述的基于消除红外滤光片折射效应的手术导航仪定位方法，其特征在于，所述步骤2包括：

步骤2.1：使用游标卡尺或螺旋测微器测量红外滤光片的厚度，多次测量后取均值，记所述左右两侧相机所用的红外滤光片的厚度分别为d_l、d_r；

步骤2.2：使用折射率测量仪测量红外滤光片的折射率为n；

步骤2.3：将红外滤光片安装在左右两侧相机的镜头前。

4.根据权利要求3所述的基于消除红外滤光片折射效应的手术导航仪定位方法，其特征在于，所述步骤3包括：

步骤3.1：使用光学导航仪采集手术器械图像，即安装红外滤光片后相机采集到的包含红外标志点的图像，其中，将左侧相机采集到的图像记为I_l；将右侧相机采集到的图像记为I_r；

步骤3.2：利用相机内参矩阵(A_l,A_r)和镜头畸变参数(S_l,S_r)对图像I_l、I_r进行镜头畸变矫正，在此只考率径向畸变和离心畸变，对左右两幅图像分别采用如下畸变校正方法进行校正，校正公式如下：

$\left[\begin{array}{c}U\\ V\end{array}\right]=\left[\begin{array}{c}{U}_d\\ V_d\end{array}\right]+\left[\begin{array}{c}{\mathrm{\δ}}_{Udr}\\ {\mathrm{\δ}}_{Vdr}\end{array}\right]+\left[\begin{array}{c}{\mathrm{\δ}}_{Udd}\\ {\mathrm{\δ}}_{Vdd}\end{array}\right];$

式中，δ_Udr表示观测到的图像点(U_d，V_d)在x方向上的径向畸变的畸变量，δ_Vdr分别表示观测到的图像点(U_d，V_d)在y方向上的径向畸变的畸变量，δ_Udd表示观测到的图像点(U_d，V_d)在x方向上的离心畸变的畸变量，δ_Vdd分别表示观测到的图像点(U_d，V_d)在y方向上的离心畸变的畸变量，U_d为观测到的图像点在x方向上的图像坐标，V_d为观测到的图像点在y方向上的图像坐标，U表示观测到的图像坐标(U_d，V_d)在经过畸变校正后在x方向上的图像坐标，V分别表示观测到的图像坐标(U_d，V_d)在经过畸变校正后在y方向上的图像坐标，其中：

$\left[\begin{array}{c}{\mathrm{\δ}}_{Udr}\\ {\mathrm{\δ}}_{Vdr}\end{array}\right]=(k_1{r_d}^2+k_2{r_d}^4)\left[\begin{array}{c}{U}_d\\ V_d\end{array}\right],$

r_d ²＝U_d ²+V_d ²

式中，r_d是从图像主点至畸变图像点(U_d，V_d)的距离，k₁表示二阶径向畸变参数，k₂表示四阶径向畸变参数。

5.根据权利要求4所述的基于消除红外滤光片折射效应的手术导航仪定位方法，其特征在于，所述步骤4包括：

步骤4.1：对镜头畸变校正后的图像I_l1、I_r1进行二值处理，得到左侧相机的二值图像I_l2和右侧相机的二值图像I_r2；

步骤4.2：利用区域提取算法对左侧相机的二值图像I_l2和右侧相机的二值图像I_r2进行处理，得到左侧相机和右侧相机采集到的图像中每个红外标志点的成像区域；

步骤4.3：使用重心法或者任何区域中心提取算法，计算红外标志点对应区域的中心，并把该中心作为红外标志点的中心在图像中对应的位置；

步骤4.4：使用对极约束法快速匹配左侧相机和右侧相机采集到的图像中的红外标志点，得到红外标志点对应的图像点对(x_l,x_r)，其中x_l表示红外标志点的中心在左图像中对应的图像点，x_r表示同一个红外标志点的中心在右图像中对应的图像点。

6.根据权利要求5所述的基于消除红外滤光片折射效应的手术导航仪定位方法，其特征在于，所述步骤5包括：

步骤5.1：根据相机的内参矩阵、红外滤光片的厚度和折射率得到相机透过红外滤光片的相机成像模型；

当以左侧相机坐标系作为世界坐标系时，即令相机坐标系的原点位于相机中心，则对于每一个图像点对(x_l,x_r)，其中 $x_l=\left(\begin{array}{c}{u}_l\\ v_l\end{array}\right)$ 表示左图像点的图像坐标， $x_r=\left(\begin{array}{c}{u}_r\\ v_r\end{array}\right)$ 表示右图像点的图像坐标，u_l表示左图像中的图像点在x方向上的坐标值，v_l表示左图像中的图像点在y方向上的坐标值，u_r表示右图像中的图像点在x方向上的坐标值，v_r表示左图像中的图像点在y方向上的坐标值；

左相机投影矩阵P_l为P_l＝A_l[It_dl]，右相机投影矩阵P_r为P_r＝A_r[Rt-t_dr]，其中：左侧相机虚拟视点在左相机坐标系中的位置为t_dl＝(00d_Ol)^T；右侧相机虚拟视点在右相机坐标系中的位置为t_dr＝(00d_Or)^T；d_Ol表示左相机对应于图像点对(x_l,x_r)的相机虚拟视点在左相机光轴上的平移距离；d_Or表示右两侧相机对应于图像点对(x_l,x_r)的相机虚拟视点在右相机光轴上的平移距离；d_Ol和d_Or可以使用如下公式进行计算：

$d_{Ol}=d_l(1-\frac{\sqrt{1-{\sin }^2{\mathrm{\θ}}_l}}{\sqrt{n^2-{\sin }^2{\mathrm{\θ}}_l}}),d_{Or}=d_r(1-\frac{\sqrt{1-{\sin }^2{\mathrm{\θ}}_r}}{\sqrt{n^2-{\sin }^2{\mathrm{\θ}}_r}}),$

式中，d_l表示左相机镜头前红外滤光片的厚度；d_r表示右相机镜头前红外滤光片的厚度，通过游标卡尺或是千分尺进行测量；n表示红外滤光片的折射率；其中，

${\mathrm{sin\θ}}_l=\frac{r_l}{\sqrt{r_l^2+f_l^2}};$

f_l＝(f_lx+f_ly)/2；

$r_l=\sqrt{{u_l}^2+v_l^2};$

${\mathrm{sin\θ}}_r=\frac{r_r}{\sqrt{{r_r}^2+f_r^2}};$

f_r＝(f_rx+f_ry)/2；

$r_r=\sqrt{{u_r}^2+v_r^2};$

式中：θ_l表示左相机中图像点对应的光线相对于红外滤光片的入射角，θ_r表示右相机中图像点对应的光线相对于红外滤光片的入射角，r_l表示左图像点到左图像主点的距离，r_r表示右图像点到右图像主点的距离，f_l表示左相机的焦距，f_r表示右相机的焦距；

步骤5.2：根据相机成像模型、左侧和右侧相机的相对位置关系以及红外标志点对应的图像点对计算红外标志点的空间位置；

对于每一对对应的图像点对(x_l,x_r)，可以使用任何传统的三角测量法进行计算，此处以传统的线性三角测量法为例进行计算，具体计算如下：

$s_l{u}_l=p_{l1}^T\hat{X},s_l{v}_l=p_{l2}^T\hat{X},s_l=p_{l3}^T\hat{X};$

$s_r{u}_r=p_{r1}^T\hat{X},s_r{v}_r=p_{r2}^T\hat{X},s_r=p_{r3}^T\hat{X};$

其中： $\begin{array}{c}{u}_l{p}_{l3}^T\hat{X}-p_{l1}^T\hat{X}=0\\ v_l{p}_{l3}^T\hat{X}-p_{l2}^T\hat{X}=0\\ u_r{p}_{r3}^T\hat{X}-p_{r1}^T\hat{X}=0\\ v_r{p}_{r3}^T\hat{X}-p_{r2}^T\hat{X}=0\end{array};$

则改写为 $A\hat{X}=0$ 时， $A=\left[\begin{array}{c}{u}_l{p}_{l3}^T-p_{l1}^T\\ v_l{p}_{l3}^T-p_{l2}^T\\ u_r{p}_{r3}^T-p_{r1}^T\\ v_r{p}_{r3}^T-p_{r2}^T\end{array}\right];$

式中：表示左侧相机投影矩阵的第一行，表示待求空间点的坐标，表示左侧相机投影矩阵的第二行，表示左侧相机投影矩阵的第三行，s_l表示比例因子，s_r表示比例因子，表示右侧相机投影矩阵的第一行，表示右侧相机投影矩阵的第二行，表示右侧相机投影矩阵的第一行。

7.根据权利要求6所述的基于消除红外滤光片折射效应的手术导航仪定位方法，其特征在于，所述步骤5中所应用的相机模型能够逐一对应图像的像素变化，即每一像素点相当于相机平移到了新的位置，相机成像模型也变为逐像素点相机中心变化的针孔相机成像模型。