CN101467887B - 一种手术导航系统中x射线透视图像标定方法 - Google Patents

Abstract

本发明属于医学图像处理及应用领域，涉及一种手术导航系统中X射线透视图像标定方法。本发明获取含有标记点信息的术中X射线透视图像，经滤波减影后，得只含标记点信息的图像后，利用模板匹配、聚类分析以及信息熵统计的方法获得标记点的中心坐标及排列方向，用B样条进行标记点坐标的校正后，进行标定参数的计算。根据得到的成像模型参数进行获得标记点的投影坐标，根据此位置与图像上实际坐标的距离，对图像的所有像素进行校正，获得变形校正后的X射线透视图像，供导航系统使用。该方法实施简单，算法可靠，便于临床应用，能提高导航系统精度。

Description

一种手术导航系统中X射线透视图像标定方法

技术领域

本发明属医学图像处理及应用领域，涉及一种用于X射线透视图像标定方法，具体涉及一种脊柱外科手术导航系统中X射线透视图像标定方法。

背景技术

基于C型臂X射线透视图像的手术导航是脊柱外科手术导航的一种主要类型，多用于胸腰段的脊柱手术。此类手术导航，必须要根据的X射线透视图像计算出成像系统的内部参数，如焦距、图像尺度因子、径向失真等，以及外部参数，如焦点空间位置、投影方向等。上述过程称之为成像系统的标定。同时由于X射线机的影像增强器球面效应以及重力导致的磁场偏移等原因，得到的透视图像存在径向失真及局部不规则的变形，而且变形的程度和方式随着增强器的位置变化而不同。该变形会导致标定计算的误差，最终影响到导航的精度，因此在标定的同时还需要解决图像的变形问题。现有技术解决这一问题的途径有两种，第一种是在标定前进行图像校正，如利用B样条，薄板样条拟和变形情况，对进行图像的恢复，标定采用恢复后图像进行计算。由于图像的变形校正与标定是两个独立的过程，导致计算误差较大。第二种是将图像几何变形作为成像系统的内部参数，直接进行标定计算，如Tsai RY在1987年的IEEE Transaction OnRobotics and Automation杂志上的A versatile camera calibration technique for high-accuracy 3d machine vision metrology using off-the-shelf TV cameras and lenses一文中，在建立标定方程时将径向失真作为内部参数计算。该方法不能解决图像局部变形的问题。

标定过程按照C型臂X射线机与患者的相对位置关系分成离线标定与在线标定两种方式。离线标定，就是将C型臂X射线机推离手术区域进行标定。由于此时所拍摄的X射线图像仅包含标定特征的信息，没有患者等其它物体的干扰，因此标定算法较简单。然而此类方法存在以下缺点：首先手术导航仅能在成像系统内部参数固定不变的透视体位下进行，限制了一些C型臂X射线机的功能应用，影响了手术导航的灵活性。其次，手术导航时，随着C型臂转动，影像增强器的位置会发生变化，因此X射线透视图像的变形会产生变化，因此按照离线标定的计算结果再进行手术导航，精度会受到影响。在线标定，就是在术中确定好导航透视体位的时候，直接在获取病人透视图像的时候完成标定，这样做解决了离线标定的缺点。但由于透视图像中除了标记点信息外还有患者身体、手术器械等影像的干扰，标定过程复杂，稳定性有待提高。

发明内容

本发明的目的是为临床应用提供一种操作方便，实施简单，定位精度高的一种实时在线的X射线透视图像标定方法，具体涉及一种脊柱外科手术导航系统中X射线透视图像标定方法。

为了获取图像标定特征，本发明利用FluroTrax标定靶(图1)实时获得术中的含有标定特征的图像。所述的标定靶由两个平行的靶面组成，使用时将其固定在C型臂X射线机的影像增强器上，靠近影像增强器一端的靶面称作近端靶面，远离增强器的另一个靶面称作远端靶面。在两个靶面上分别镶嵌有规则排列的钢球(图2)，其中在近端靶面上还有数个体积较大的钢球排列成U型，用来确定钢球的排列规则与位置。钢球的空间坐标通过红外线定位仪追踪标定靶侧面的发光二极管得到。通过本发明对标定特征的自动识别以及与图像变形校正紧密结合的标定计算，可以实时地获得病人当前透视体位下的光学成像参数参数，并进行手术导航，能解决图像变形与离线标定对手术导航精度产生的不利影响。

本发明的技术方案是：首先在手术中采集含有标定靶信息的透视图像，经过滤波减影预处理后得到只包含标记点信息的图像；然后利用模板匹配算法、聚类分析以及信息熵统计方法对图像中的标记点进行自动识别，获得标记点的中心坐标及排列方向；采用B样条拟合的方法对标记点坐标进行校正；然后根据红外线定位仪获得的标记点的空间及其在图像上校正后的坐标，计算透视系统的成像参数；最后再根据上述参数将标定靶上标记点投影到图像中，根据投影坐标与实际图像上的坐标距离，利用B样条拟合的方法对图像中的全部像素进行校正，形成手术导航的参考图像。

本发明的技术方案通过下述步骤实现：

(1)X射线透视图像采集及预处理；

(2)采用模板匹配、聚类分析及信息熵统计方法对透视图像上的标记点进行标记点的识别；

(3)利用B样条对透视图像中的标记点进行坐标校正；

(4)采用红外线定位仪获得标记点的三维空间坐标，结合其在图像上的二维坐标进行标定计算，获得成像模型参数；

(5)根据成像模型参数，计算标记点的投影坐标，根据投影坐标与实际的坐标间的误差，利用B样条对透视图像进行校正，生成手术导航的参考图像。

本发明将标定与图像变形校正紧密结合，既保证了成像模型预测精度，又解决了变形难以估计这一难题，从而可以在临床上实施，大幅度提高手术导航系统的精度。本发明方法实施简单，算法可靠，操作灵活，无须更改现有的C型臂X射线透视设备，能提高导航系统精度，便于临床应用。

附图说明

图1是FluroTrax标定靶。

图2是靶面标记点，其中a：远端靶面标记点，b：近端靶面标记点。

图3是标记点横向位移校正示意图。

图4照相机成像模型。

具体实施方式

实施例1

1、X射线透视图像采集及预处理

使用视频采集技术将X射线透视透视图像数字化，采用多帧平均滤波及小波变换的方法进行噪声综合抑制。首先对图像进行多幅叠加平均，然后将平均后的图像作为输入图像进行小波变换滤波。小波的类型采用由Daubechies提出的对称小波symN，N＝2，3，…8，采用软阈值化，固定阈值准则方式进行去噪处理。经过上述系列处理后，图像噪声被有效滤除并且较好地保留了原始图像中的细节信息。

2、标记点识别与校正

首先对背景信息进行过滤。利用数字减影的原理：设原图像为I_Source，利用以标记点直径为大小的核函数对I_Source进行中值滤波处理，得到不含标记点的背景图像I_Background。将I_Source与I_Background相减，获得只含标记点信息的图像。

然后利用模板匹配方法对标记点进行识别，模板采用事先根据标准图像分割好的标记点图像。设被搜索图像S的大小为N×N，模板T的大小为M×M。把模板T在搜索图像S上移动，在模型覆盖下的搜索图称子图S^i，j，i，j为该块子图的左上角像素在S图像中的坐标，成为参考点。比较T和S^i，j的内容，若两者一致，则两者之差为零。本发明采用如下测度衡量T和S^i，j的相似程度：

$R(i,j)=\frac{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{M}{\mathrm{\Σ}}}[S^{i,j}(m,n)\×T(m,n)]}{\sqrt{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{M}{\mathrm{\Σ}}}{\left[S^{i,j}(m,n)\right]}^2}\sqrt{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{M}{\mathrm{\Σ}}}{\left[T(m,n)\right]}^2}}$

根据许瓦尔兹不等式，0≤R(i，j)≤1，当且仅当S^i，j(m，n)＝kT(m，n)时，R(i，j)＝1，其中，k是标量常数。

由于图像与标准模板间的差异，对R(i，j)设定一个阈值，凡是大于该阈值的则认为可能是标记点。根据所述标记点之间的距离采用聚类方法将其合并，并得到每个类的类中心，在类中心周围小区域内采用如下搜索算法确定标记点的中心位置：

$P_c=\frac{\underset{p\∈\mathrm{\Ω}}{\mathrm{\Σ}}(I\left(p\right)-I_T)P}{\underset{p\∈\mathrm{\Ω}}{\mathrm{\Σ}}(I\left(p\right)-I_T)}$

其中P_c是标记点中心坐标，I(p)为p点的密度，I_T为目标标记的分割阈值，Ω是p点的小邻域。通上述方法确定的标记点中心精度可以达到到亚像素水平。

为了获得标记点排列方向，本发明将图像旋转到其标记点排列方向与图像坐标轴的方向一致。本发明利用信息熵H(θ_j)来确定旋转角度，该函数定义如下：

$H\left({\mathrm{\θ}}_j\right)=-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}p\left(x_i\left({\mathrm{\θ}}_j\right)\right){\log }_{2}p\left(x_i\left({\mathrm{\θ}}_j\right)\right)$

其中θ是旋转角度，x_i(θ_j)是图像向X轴的投影值。由熵的定义我们得到结论：当图像中标记点的列方向与Y轴一致时，每一列的标记点在X轴上的投影将重合成一点，因此这时的信息量最少，其H(θ_j)值也必然最小。

3、标记点的位置校正

标记点的位置校正通过下述步骤：

1)估算在未变形情况下近端靶面标记点在图像上的投影坐标，

利用正方形网格T表示在未变形情况下近端靶面标记点在图像上的投影位置。T通过以下三个参数确定：网格中心，排列方向，网格间距。其中网格中心由近端靶面标记点投影坐标的平均值得到，排列方向已经在步骤2中获得。网格间距根据近端靶面的实际间隔与图像缩放因子相乘得到，或者根据对图像上标记点间距的测量值直接得到；

2)根据上述估算的坐标与实际投影坐标的距离，对远端靶面标记点的实际投影坐标进行校正，

首先，以T的中心为旋转中心，旋转图像至网格排列方向与坐标轴一致，其旋转的角度在步骤2中已经给出，然后利用基于B样条的坐标校正算法对图像中的远端标记点位置进行校正，最后将校正后的图像旋转回原来的方向。

所述的基于B样条的坐标校正算法是：初始的数据有参考网格T以及源网格S。S由实际图像上(变形后的)的近端标记点中心构成。S与T的大小一致，且标记点是一一对应的。根据S和T构造新的网格I，I中结点的x坐标来自于T，y坐标来自于S。

第一步进行横向位移的校正：以图像中的标记点为控制点，分别对S与I的每一列构造纵向B样条X_S、X_I。对于所有远端靶面标记点的投影(变形后)位置(x_fiducial，y_fiducial)，插入水平线Y＝y_fiducial。该线与X_S、X_I分别相交于(x_S，y_fiducial)及(x_I，y_fiducial)，这样在直线Y上x_s点处的横向校正位移为(x_I-x_S)。{(x_S，，(x_I-x_S)}是直线Y与每条X_S的交点以及该点的横向校正位移组成的点集。根据该点集再构造一个横向的B样条。这样，标记点(x_fiducial，y_fiducial)的横向校正位移是此样条在x_fiducial坐标的插值计算结果(图3)。这样依次对远端标记点进行校正，确定其最终的x坐标；

第二步进行纵向位移的校正，以图像中的标记点为控制点，分别对I与T的每一行构造横向B样条Y_I、Y_T。对每个远端靶面变形的标记点(x_fiducial，y_fiducial)，插入垂直线X＝x_fiducial。该线与Y_I、Y_T分别相交于(x_fiducial，y_I)及(x_fiducial，y_T)，这样在该垂直线的y_S点处的纵向校正位移为(y_T-y_S)。{(y_S，，(y_T-y_S)}是直线X与每条Y_I的交点以及该点的纵向校正位移。根据该点集再一个构造竖直的B样条。这样，标记点(x_fiducial，y_fiducial)的纵向校正位移是此样条在y_fiducial坐标的插值计算结果。这样依次对远端标记点进行校正，确定其最终的y坐标。

至此，得到全部标记点的在图像上的变形校正后的投影坐标。

4、标定参数计算

X射线成像模型同一般小孔相机成像模型(图4)。照相机C具有焦点c，即X射线的焦点，以及成像平面plane(C)＝(o，h，v)。

其中o是图像坐标系的原点，h，v是图像平面的坐标轴方向。点p是三维空间中的一点，其在2D图像平面的投影点为p′。焦点c在像平面上的投影点为c′。与像平面垂直的轴为a＝h×v。s_x，s_y分别为图像沿h，v方向近端靶面相邻两点间隔像素数目。p_y是p点在过c点与像平面垂直的平面上的投影点。p_y′是p_y在像平面上的投影点。p_a是p_y在直线c c′上的投影点。

根据上述模型建立3D空间到2D图像空间的坐标转换方程。

三角形Δcpp_y与Δcp′p_y′是相似三角形，因此

$\frac{\left|p^{\′}{p}_y\right|}{\left|{\mathrm{pp}}_y\right|}=\frac{\left|{{\mathrm{cp}}_y}^{\′}\right|}{\left|{\mathrm{cp}}_y\right|}$ ...公式1

三角形Δcp_ap_y与Δcp′p_y′也是相似三角形，有

$\frac{\left|{{\mathrm{cp}}_y}^{\′}\right|}{\left|{\mathrm{cp}}_y\right|}=\frac{\left|{\mathrm{cc}}^{\′}\right|}{\left|{\mathrm{cp}}_a\right|}$ ...公式2

因此，由公式1、2得到：

$\frac{\left|p^{\′}{p_y}^{\′}\right|}{\left|{\mathrm{pp}}_y\right|}=\frac{\left|{\mathrm{cc}}^{\′}\right|}{\left|{\mathrm{cp}}_a\right|}$ ...公式3

设f为相机焦距，(x′，y′)是p′的坐标，(x₀，y₀)是c′的坐标，公式3变换为

$\frac{(x^{\′}-x_0)s_x}{{(p-c)}^{T}h}=\frac{f}{{(p-c)}^{T}a}$ ...公式4

定义f_x＝f/s_x，将f_x代入公式4得：

$x^{\′}=\frac{{(p-c)}^T(f_{x}h+x_{0}a)}{{(p-c)}^{T}a}$ ...公式5

定义

$\tilde{h}=f_{x}h+x_{0}a$

$\tilde{p}=p-c$

代入公式5得到点p在图像上投影点的x坐标：

$x^{\′}=\frac{{\tilde{p}}^T\tilde{h}}{{\tilde{p}}^{T}a}$ ...公式6

同理，定义

f_y＝f/s_y

v′＝f_yv+y₀a

则点p在图像上投影点的y坐标为：

$y^{\′}=\frac{{\tilde{p}}^T\tilde{v}}{{\tilde{p}}^{T}a}$ ...公式7

由此得到三维空间坐标到二维图像坐标的转换关系。

标定参数计算即已知标记点的三维空间坐标与二维图像坐标的位置，求解成像模型参数的过程。首先根据已知靶面上标记点的空间位置及图像上其对应的位置，采用视觉标定中的直接线性变换(DLT)算法求出图像坐标到靶面坐标的的转换矩阵。

在透视投影的情况下，像点与物点的成像几何关系在齐次坐标下可以表示成：

P＝HQ

其中P是靶面标记点的三维空间坐标[x，y，z]^T，Q是其在图像上的二维空间坐标[x′，y′，1]^T，H是要求的3×3空间转换矩阵。利用最小二乘法原理，将靶面上所有标记点代入方程，用奇异值分解法求得H：

H＝P·Q^T·(Q·Q^T)-¹

然后根据该矩阵，可以得到图像上任意一点在靶面上的位置。据此计算出图像上若干个点在近端及远端靶面上的空间位置，构成多条虚拟的射线R＝{r₀，r₁，…}，用最小二乘法计算出这些射线的交汇处，这就是X射线的焦点位置。然后根据公式6与公式7计算投影方向、投影中心点等其它参数。

5、图像的校正

标定完成后要将透视图像中的所有像素进行校正，生成一张每个像素在原始图像上的坐标查找表。根据该表对原始图像重新插值采样，生成变形校正后图像，作为手术导航的参考图像。

该校正算法步骤3中标记点的位置校正类似，不同之处在此步骤是要对所有像素进行校正，并且是寻找校正后图像的每个整数像素点在原始图像上的位置，以便在原始图像上插值采样，因此校正方向与步骤3的算法相反。具体描述如下：

校正算法初始的数据分别是由源网格S、目标网格T以及透视图像组成。所述的S的结点由原始图像中(变形后)的近端及远端靶面的所有标记点的中心位置构成；T的结点是根据标定参数计算出的近端及远端靶面的所有标记点的投影位置组成。根据S与T构造中间网格I，其结点的x坐标来自S，y坐标来自T。首先按照上述步骤3中的方法将S、T、I及透视图像及旋转至网格排列方向与坐标轴一致，然后进行校正：

第一步，计算每个像素在校正后与校正前x方向上的位移。设T与I共有W_grid列与H_grid行，透视图像的像素有W_image列以及H_image行。以T与I的每一列的点为控制点构造竖直的B样条，然后在行方向上对每个B样条进行采样，分别建立表T_T与T_I，其大小为W_grid×H_image。这样在所有坐标为T_T(i，j)点的横向校正位移是T_i-T_T。对属于j的每一行，0＜j＜h_image，以下列点为控制点构造横向B样条：

{(T_T(i，j)，T_I(i，j)-T_T(i，j))|0≤i ≤w_grid}

那么像素(x_pixel，j)在此行的横向校正位移是此B样条在x_pixel处的插值结果。

第二步，计算每个像素在校正后与校正前y方向上的位移。设S与I共有W_grid列与H_grid行，透视图像的像素有W_image列以及H_image行。以S与I的每一行的点为控制点构造横向的B样条，然后在列方向上对每个B样条进行采样，分别建立表S_S与S_I，其大小为W_image×H_grid。这样在所有坐标为S_S(i，j)点的横向校正位移是S_I-S_T。对属于i的每一列，0＜i＜W_image，以下列点为控制点构造纵向B样条：

{(S_S(i，j)，S_I(i，j)-S_S(i，j))|0≤j≤H_grid}

那么像素(i，y_pixel)在此行的纵向校正位移是此B样条在y_pixel处的插值结果。

根据每个像素的x，y方向的校正位移构造一张像素坐标查找表，根据此表，对原始图像进行重采样，得到变形校正后的X射线透视图像，将此做为手术导航的参考图像。

本发明将标定与图像变形校正紧密结合，既保证了成像模型预测精度，又解决了变形难以估计这一难题，从而可以在临床上实施，大幅度提高手术导航系统的精度。

实施例2临床透视图像的采集、标定试验

将标定靶固定在C型壁X射线机的影像增强器上，近端靶面贴近增强器输入屏，进行患者透视图像的采集。考虑到采样时间及患者身体的移动，本发明将透视踏板松开后前0.5秒内采集的图像进行叠加平均，进行小波滤波；

对标记点进行自动识别，中值滤波器的核函数为2.5mm；本发明利用滤波以及减影方法去除图像中非标记点的信息。在模板匹配算法中，以0.85为阈值，凡R(i，j)＞0.85的点均列为可能的标记点中心，进行聚类处理，寻找类中心，然后通过求类中心的重心法获得标记点的具体位置，经过信息熵统计的方法确定图像旋转角度为-16°；

将图像旋转-16°后进行标记点校正，然后旋转16°，回到原来的方向；

利用红外线定位仪跟踪标定器件上的发光二极管，获得标记点的三维空间位置；根据标记点的三维坐标以及之前得到的其在图像上的二维坐标进行标定计算，获得射线焦点位置，焦距，投影平面的坐标，投影中心点坐标等参数信息；

按照标定计算的参数将标记点全部投影到图像上，与透视图像一起旋转-16°，使标记点排列方向与网格一致，然后对所有像素进行校正，生成一张校正像素的坐标查找表；。根据此表对原始图像重新插值采样，生成新的变形校正图像，旋转16°后作为手术导航的参考图像；

将成像模型参数与参考图像输出，供手术导航系统使用，在奔IV2.6G、1G内存、WinXP操作系统下运行，需要15秒左右。

Claims (6)

1.一种手术导航系统中X射线透视图像标定方法，其特征是包括下述步骤：

(1)X射线透视图像采集及预处理；

(2)采用模板匹配、聚类分析及信息熵统计方法对透视图像上的标记点进行标记点的识别；

(3)利用B样条对透视图像中的标记点进行坐标校正；

(4)采用红外线定位仪获得标记点的三维空间坐标，结合其在图像上的二维坐标进行标定计算，获得成像模型参数，

所述的成像模型参数包括X射线的焦点位置、投影方向、投影中心点和成像平面；

(5)根据成像模型参数，计算标记点的投影坐标，根据投影坐标与实际的坐标间的误差，利用B样条对透视图像进行校正，生成手术导航的参考图像。

2.根据权利要求1所述的方法，其特征是，所述的步骤(1)的X射线透视图像是目标图像。

3.根据权利要求1所述的方法，其特征是，其中步骤(2)的标记点的识别采用下述方法：利用中值滤波获得背景图像，应用数字减影以及模板匹配的方法，得到标记点中心，用聚类分析与图像搜索算法得到标记点中心坐标，用信息熵统计得到网格方向。

4.根据权利要求1所述的方法，其特征是，其中步骤(3)的B样条校正标记点通过下述算法：根据S与T上对应点的距离之差，分别进行横向位移的校正与竖直位移的校正，所述的T是由估算的未变形情况下近端靶面上标记点在图像上的投影位置构成；S由实际图像上的近端标记点中心位置构成。

5.根据权利要求1所述的方法，其特征是，其中步骤(4)所述的X射线焦点的求解算法为：

第一步：D算法求出图像点到靶面点的的坐标转换矩阵，在透视投影的情况下，

像点与物点的成像几何关系在齐次坐标下表示成：

P_prox＝H_proxQ_prox

其中，P_prox是近端靶面标记点的三维空间坐标[x，y，z]^T，Q_Prox是其在图像上的二维空间坐标[x′，y′，1]^T，H_Prox是要求的3×3空间转换矩阵，利用最小二乘法原理，将靶面上标记点代入方程，用奇异值分解法求得H_Prox：

H_prox＝P_prox·Q_prox ^T·(Q_prox·Q_prox ^T)^-1

同理，求得远端靶面标记点的坐标转换矩阵H_dist；

第二步：根据转换矩阵确定图像上一组点在两个靶面上对应点的坐标，然后将每两个对应点相连，构成一组射线，用最小二乘法计算出这些射线的交汇处，得X射线的焦点。

6.根据权利要求1所述的方法，其特征是，其中步骤(5)所述的透视图像进行校正算法为，根据S与目标网格T对应点的距离之差，分别进行横向位移的校正与竖直位移的校正，在每一项校正中，构造两类B样条，最后生成图像中每一个像素在原始图像上的坐标查找表，根据所得的表，对原始图像重构，得到变形校正后的X射线透视图像；

所述的S由原始图像中所有标记点的中心位置构成；目标网格T是根据标定参数计算得的近端及远端靶面的所有标记点的投影位置组成。

Images