CN101582161A - 一种基于透视成像模型标定的c型臂图像校正方法 - Google Patents

Abstract

一种基于透视成像模型标定的C型臂图像校正方法，它有四大步骤；步骤一：设计C型臂系统标定的标定模板；步骤二：获取标记点的失真坐标，在失真图像上通过图像处理获得；步骤三：C型臂系统标定；包括标定模板上标记点的三维空间坐标、透视成像模型标定、畸变模型和失真坐标通过系统标定获得透视成像参数；步骤四：图像校正；由标定的透视成像模型，对失真图像进行校正。本发明在线将C型臂系统标定和C型臂X射线投影失真图像校正作为一个整体，基于摄像机透视成像模型对C型臂系统进行标定和对C型臂X射线投影失真图像进行校正，改变了在传统方法中先离线进行图像校正然后进行线性标定，使系统步骤变得简洁方便，容易在线使用。

Description

一种基于透视成像模型标定的C型臂图像校正方法

(一)技术领域

本发明涉及一种基于透视成像模型标定的C型臂图像校正方法，它主要是把C型臂系统成像参数标定和C型臂X射线投影失真图像校正作为一个整体，基于透视成像模型来标定C型臂系统获得C型臂系统的内外参数，再利用标定C型臂系统获得的畸变参数，实现对C型臂X射线投影失真图像校正的一种图像校正方法，属于图像处理领域。

(二)背景技术

图像校正是一种常见的图像处理方式，如遥感图像校正后进行图像镶嵌，航拍图像运用前进行校正，摄像机图像校正进行定位或三维重建。

目前C型臂X射线投影图像存在三种类型的失真变形：枕形失真，S形失真，图像局部失真。影像增强器的输入屏为弧形曲面，人体组织器官透视后残余的X线投影到该曲面上，无法到达预定的位置从而导致最终生成的图像存在枕形失真，影像增强器中加速后的电子束受地球及周围电子设备磁场影响发生偏转和平移，从而导致输出图像产生“S”形失真，X线图像成像系统器件的加工安装误差，尤其是MOS在面阵CCD半导体衬底上的排列定位误差，导致输出图像产生局部变形失真。此外，C型臂影像增强器的重量会造成C型臂电枢中心产生一定程度的偏移，从而导致图像出现偏移。C型臂x射线投影图像的失真变形，使图像中病人的解剖结构信息发生扭曲变形，在直接用于手术导航和器械定位时，会影响定位的精度，影响手术的结果，严重时会出现医疗事故。因此，基于C型臂X射线投影图像导航的手术器械定位和三维重建辅助手术系统必须对C型臂X射线投影失真图像进行校正。

卡斯帕尔(Casperson)定性分析与定量测量了影像增强器所成的图像，指出枕形失真是C型臂x射线投影图像变形的主要形式。纳赛尔(Nassir)研究了C型臂的成像模型，提出C型臂的成像模型与摄像机的针孔模型相似。因此C型臂X射线投影图像的失真模型可以用摄像机的失真模型模拟。

传统方法在处理基于C型臂X射线投影图像导航时，分成了两个阶段：C型臂X射线投影失真图像离线校正和C型臂系统标定。此方法有优点也有缺点。

C型臂X射线投影失真图像校正有多种方法，通常可分为全局法和局部法。全局法是对整幅图像使用同一个最优的单一变换模型，以确定理想图像与失真图像中对应像素之间的空间映射关系和灰度值映射关系。具体的全局校正算法包括刚体变换、仿射变换、透视变换以及多项式变换等；局部法则是在图像中局部控制点的基础上，构造映射函数，以确定理想图像中局部区域像素与失真图像中对应的局部区域像素之间的空间映射关系和灰度值映射关系。在传统方法中，首先离线对C型臂X射线投影失真图像进行校正，而且大部分采用的是高阶多项式拟合法，精度有时能达到很高。但是这种方法需要建立理想图像，也就是说高阶多项式拟合图像校正实际是拟合失真图像像素点和对应理想图像像素点的函数关系，但是实际运用中绝对理想的理想图像是很难得到的，一般都是用近似理想图像替代，这就造成图像校正的结果与近似理想图像的好坏有很大关系。传统方法中，图像校正后再使用简单的线性模型对C型臂系统进行标定来获取C型臂系统的内外参数。优点是计算简单，缺点是由于C型臂系统中影像增强器和标定模板的自身重量影响，造成标定模板中心，影像增强器中心和放射源不在同一直线上。如果继续使用线性模型标定C型臂系统，那么获得的C型臂系统参数将不准确，严重影响导航的精度。

(三)发明内容

本发明的目的是提供一种基于透视成像模型标定的C型臂图像校正方法，该方法不会因标定模板中心，影像增强器中心和放射源不在同一直线上而出现很大误差，也不需要先进行C型臂X射线投影失真图像校正再进行C型臂系统标定，而且理论上可以在线实现C型臂系统标定的同时实现C型臂X射线投影失真图像校正，从而可以节省时间，实现实时性。

本发明是一种基于透视成像模型标定的C型臂图像校正方法，它是在以下硬件载体上实现的：PC机、双目摄像机、支架、C型臂系统、标记板、标定模板以及固定卡爪。以上硬件的组合方式如图7所示：双目摄像机是建立参考坐标系，确定标定模板上标记点的三维空间坐标；C型臂系统是采集标定模板X射线投影图像；支架是支撑和固定双目摄像机；固定卡爪是固定标定模板在影像增强器外壳上；标记板是建立标定模板中标记点的坐标系；PC机是处理输入的C型臂X射线模板投影失真图像和图像中标记点对应的摄像机坐标系下标记点的三维空间坐标，从而实现C型臂系统标定和C型臂X射线投影失真图像校正。

其中，PC机可选用个人计算机或工控机。个人计算机作为悉知的普及性计算机性能越来越高，已经完全满足这个方法的要求。而工控机性能介于一般，但稳定性比个人计算机要高

其中，双目摄像机可选由单摄像机组合或双目摄像机。两个单摄像机组合需标定相对位置，双目摄像机相对关系是已知的，直接可以使用。如北京微视(MV)工业双目摄像机、深圳博安盈科技有限公司iBoin双目摄像机等等。

其中，支架可以定做。

其中，C型臂系统有大C型臂和小C型臂。如上海冠东/上海杰瑞公司的介入型小C臂X光机，GE公司LVC型大C臂X光机。

其中，标记板可以是圆形也可以是方形格等。

其中，立体标定模板可以使用不透过X射线的材料制作。比如铝块和电路板。材料质地越轻，对影像增强器的成像影响就越小。

其中，固定卡爪与C型臂罩配套，目的是为了固定标定模板在影像增强器的外壳上且可防损坏影像增强器。

本发明是一种基于透视成像模型标定的C型臂图像校正方法，其实现方法流程框图(即C型臂系统标定和C型臂X射线投影失真图像校正的整体方案示意图)如图8所示，基于上述流程的硬件组装(如图7)：双目摄像机固定在支架上，支架放置在病床的一侧，标记板粘贴于标定模板的一侧，标定模板通过固定卡爪固定在C型臂影像增强的外壳上，且有标记板的一侧面向双目摄像机，C型臂系统采集的X射线投影图像和以双目摄像机坐标系为基准的模板上标记点的三维坐标输送到另一个隔离间的PC机进行处理。本发明一种基于C型臂透视成像模型标定和C型臂图像校正方法，它由以下具体步骤实现：

步骤一：设计C型臂系统标定的标定模板(如图9)。具体为：测量C型臂影像增强器外壳直径大小，以C型臂影像增强器外壳直径451mm为直径制作一个C型臂罩(如图1)。C型臂罩上安装有两层校准靶，其间距为100mm；下层校准靶直径为380mm，且在直径为350mm的圆形区域内，有以间距20mm均匀规则分布的225个圆形小孔(如图2)。其中圆心和两个相互垂直的直径上有四个小孔的直径为6mm，其余直径为2mm；上层校准靶直径为380mm，且在直径为350mm的圆形区域内分布有9个直径为2mm的圆形小孔。其中圆心小孔垂直对应上层校准靶圆心孔，其余8个小孔分成两组，每组相距圆心距离分别是和

(如图3)。

步骤二：获取标记点的失真坐标。在失真图像上通过图像处理获得，具体为：

①.对标定模板X射线投影图像使用中值滤波去噪平滑处理。

②.使用形态学中的开运算提取标定模板背景图像。

③.提取标定模板X射线投影图像中的标记点。用标定模板X射线投影图像减去标定模板背景图像得到只含圆形标记点的图像。随后经过二值化，开运算和闭运算去除噪声得到只含圆形标记点的二值图像。

④.计算出这些圆形标记点的质心作为标记模板上标记点在X射线投影图像中的图像坐标(如图6)。

步骤三：C型臂系统标定。包括标定模板上标记点的三维空间坐标、透视成像模型、畸变模型和失真坐标通过系统标定获得透视成像参数，具体为：

i)标定模板上标记点的三维空间坐标。此步骤由双目摄像机定位标记板和不同坐标系之间的平移完成，具体为：

②「双目摄像机固定在支架上，支架放置在病床的一侧，标记板粘贴于标定模板的一侧，标定模板通过固定卡爪安装在C型臂影像增强器的外壳上，且有标记板的一侧面向双目摄像机(如图7)。

②标记板的X轴垂直上层校准靶的圆形平面，Z轴指向上层校准靶的圆心，Y轴与上层校准靶的圆边相切，标记板的原点O₁为X轴，Y轴和Z轴的交点。设定新坐标系原点在标记板原点与Y方向单元小格的中间，且以这点为原点O建立坐标系，X轴垂直上层校准靶的圆形平面，方向与标记板X轴单位矢量平行，Z轴指向上层校准靶的圆心，Y轴方向与标记板Y轴单位矢量平行。原点O₁在摄像机的坐标系中为(X，Y，Z)，Y方向中离原点一个单元格的坐标在摄像机的坐标系中为(X₂，Y₂，Z₂)，X方向中离原点一个单元格的坐标在摄像机的坐标系中为(X₁，Y₁，Z₁)，新坐标系原点O在摄像机的坐标系中为((X₂-X)/2，(Y₂-Y)/2，(Z₂-Z)/2)。把新坐标系移到上层校准靶的圆心，且以这个坐标系为参考，逐步确定上层校准靶中小孔的圆心在摄像机坐标系下的三维空间坐标。然后再把新坐标系移到下层校准靶的圆心，且以这个坐标系为参考，逐步确定下层校准靶中小孔的圆心在摄像机坐标系下的三维空间坐标。关系式如下：

(X₀₁，Y₀₁，Z₀₁)＝(X_标，原+X_摄，Y_标，原+Y_摄，Z_标，原+Z_摄)    (式8)

其中：(X₀₁，Y₀₁，Z₀₁)是标定板原点在摄像机坐标系下的坐标；

(X_摄，Y_摄，Z_摄)是标定板原点在摄像机坐标系下的偏移量；

假设的原点O：

(X₀，Y₀，Z₀)＝(X₀₁+X₂，Y₀₁+Y₂，Z₀₁+Z₂)/2    (式9)

X，Y方向的单位矢量：

$\stackrel{\→}{x}=\frac{(X_1-X_{01},Y_1-Y_{01},Z_1-Z_{01})}{\sqrt{{(X_1-X_{01})}^2+{(Y_1-Y_{01})}^2+{(Z_1-Z_{01})}^2}}$ (式10)

$\stackrel{\→}{y}=\frac{(X_2-X_{01},Y_2-Y_{01},Z_2-Z_{01})}{\sqrt{{(X_2-X_{01})}^2+{(Y_2-Y_{01})}^2+{(Z_2-Z_{01})}^2}}$ (式11)

坐标系平移到标定模板上层平面的圆心：

$\stackrel{\→}{z}=\stackrel{\→}{y}\×\stackrel{\→}{x}$ (式12)

$({X_0}^{\′},{Y_0}^{\′},{Z_0}^{\′})=(X_0,Y_0,Z_0)-\stackrel{\→}{z}\×225$ 式13)

其中：(X₀′，Y₀′，Z₀′)是坐标系在标定模板上层平面圆心处的原点

$({X_0}^{\′\′},{Y_0}^{\′\′},{Z_0}^{\′\′})=({X_0}^{\′},{Y_0}^{\′},{Z_0}^{\′})-\stackrel{\→}{x}\×112$ (式14)

其中：(X₀″，Y₀″，Z₀″)是坐标系在标定模板下层平面圆心处的原点

ii)透视成像模型标定，畸变模型和失真坐标通过系统标定获得透视成像参数，具体步骤为：

①.提取C型臂标定模板X射线投影图像在图像中心区域的标记点，并用Tsai两步法进行标定。该方法有两步：第一步利用最小二乘法求解超定线性方程，给出外部参数；第二步使用优化搜索求解内部参数。摄像机坐标系与世界坐标系关系式表示矢量形式：

[x_wY_d y_wY_d z_wY_d Y_d -x_wX_d -y_wX_d -z_wX_d]·[r₁/T_y r₂/T_y r₃/T_y T_x/T_y r₄/T_y r₅/T_y r₆/T_y]’＝X_d

其中：x_w，y_w，z_w世界坐标；

x_d，y_d，z_d摄像机坐标‘

r_i(i＝1，2，3，4，5，6)旋转矩阵分量；

T_y为Y方向的平移分量；(式1)

实际图像坐标到摄像机坐标的变换：

u_d＝s_xd_x′^-1X_d+u₀，v_d＝d_y ^-1Y_d+v₀    (式2)

其中：d_x为摄像机在X方向的像素间距；

d_y为摄像机在Y方向的像素间距；

s_x为图像尺度因子；

(u₀，v₀)为光学中心；

图像成像的理想位置与图像成像失真偏移的实际位置关系：

x′＝x+δ_x，y′＝y+δ_y，δ_x(x，y)＝k₁x(x²+y²)，δ_y(x，y)＝k₂y(x²+y²)    (式3)

其中：(x，y)为图像成像的理想位置；

(x′，y′)为图像成像失真偏移的实际位置；

δ_x，δ_y为非线性畸变值；

k₁，k₂为径向失真畸变参数；

②.在整个图像中使用Levenberg-Marquadt迭代优化算法，精确求解C型臂系统的内外参数R，T和k₁，k₂，p₁，p₂，s₁，s₂。这里迭代的初始值设置为Tsai法得到的外参数值R，T和k₁，k₂，p₁＝0，p₂＝0，s₁＝0，s₂＝0，标记点是从图像中提取的所有的标记点。图像成像失真的偏移量：

δ_x(x，y)＝k₁x(x²+y²)+(p₁(3x²+y²)+2p₂xy)+s₁(x²+y²)，δ_y(x，y)＝k₂y(x²+y²)+(p₂(3y²+x²)+2p₁xy)+s₂(x²+y²)

其中：p₁，p₂为切向失真畸变参数；

s₁，s₂为薄透镜失真畸变参数；(式4)

实际真实坐标(x′，y′)和理想真实坐标(x，y)的关系改成向量相乘的形式：

AP＝X′-X    (式5)

其中： $A=\left[\begin{array}{cccccc}x(x^2+y^2)& 0& {3x}^2+y^2& 2\mathrm{xy}& x^2+y^2& 0\\ 0& y(x^2+y^2)& 2\mathrm{xy}& x^2+{3y}^2& 0& x^2+y^2\end{array}\right]$

P＝[k₁ k₂ p₁ p₂ s₁ s₂]′，X′＝[x′y′]′，X＝[x y]′

用Levenberg-Marquadt算法优化精确求解，就是求目标函数：

$\min {\left|\right|\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(m_i-m(H,P,R,T,M_i))\left|\right|}^2$ (式6)

假设模板图像上有n个标记点，其中m_i是第i个像素点图像坐标，M_i是第i个像素点对应的三维空间坐标，H内部参数矩阵。

③再次在整个图像中使用Levenberg-Marquadt迭代优化算法，精确求解C型臂系统的内参数f，u₀，v₀。迭代的初始值设置为Tsai法得到的内参数值f，u₀，v₀。

④循环迭代这个步骤直至迭代次数完成或达到误差限而终止程序运行且输出结果。

步骤四：图像校正。由标定的透视成像模型，对失真图像进行校正；具体步骤是从透视成像模型标定C型臂系统中获得的失真参数对失真图像进行图像校正，也就是根据图像成像的理想位置与图像成像失真偏移的实际位置关系，把标定获得的失真畸变参数代入其中，获得校正图像的理想位置，最后进行灰度映射得到校正图像。如果所得到的校正图像的理想位置不是整数，那需要采用插值算法对图像进行灰度赋值。这里采用的是最近邻插值算法，关系如下：

f(x)＝f(x_k)，如果1/2(x_k-1+x_k)＜x＜1/2(x_k+x_k+1)    (式7)

本发明与现有技术相比的优点在于：传统的C型臂X射线投影失真图像的校正方法是先离线建立一个理想图像，再用多阶(4阶或5阶)多项式拟合失真图像像素点和理想图像像素点之间的关系进行图像校正。这个方法的难点是需要建立一个理想图像，理想图像越好，校正效果也就越好。但实际上绝对理想的理想图像是不可能得到的，因此在实际运用中就常使用近似理想图像来代替。近似理想图像越接近理想图像，图像校正就越好，反之就越差。在传统方法中，一般是先进行C型臂X射线投影失真图像校正后再用线性模型对C型臂系统标定获取系统的内外参数。但是由于C型臂影像增强器和校正模板的重量，造成C型臂电抠偏移，放射源，C型臂影像增强器和标定模板的中心不在同一直线，图像将发生偏移，使用线性模型标定C型臂系统时将获得不准确的系统参数。换句话说即使出厂时能够建立比较准确的近似理想图像但时间一长这个近似的理想图像也将变得不准确，再用线性模型标定结果将变得更不准确。本发明就解决了这个问题。基于透视成像模型标定的C型臂图像校正方法不同于传统的先离线校正图像再进行C型臂系统标定的方法，这个方法是在线把C型臂系统标定和图像校正作为一个整体，先使用透视成像模型对C型臂系统进行标定再用标定获得到的C型臂系统内外参数对图像进行校正。其中标定时先使用Tsai法非线性透视成像模型和图像中心区域的标记点对C型臂系统进行标定来获取C型臂系统的内外参数值，然后以这作为初始值，用图像中所有的标记点和优化算法优化来精确地获得C型臂系统的内外参数值。

(四)附图说明：

图1C型臂罩示意图(a-主视图b-俯视图)

图2下层校准靶示意图(a-俯视图b-剖面图)

图3上层校准靶示意图(a-俯视图b-剖面图)

图4固定卡爪示意图(a-主视图b-俯视图)

图5标记板示意图

图6C型臂X射线标定模板投影图中提取的标记点及质心示意图

图7C型臂标定和图像校正系统示意图

图8本发明所述方法的流程框图

图9标定模板立体示意图

图10放射源到探测器124cm，旋转0度时，标定模板的原始图和校正图

图11放射源到探测器124cm，旋转30度时，标定模板的原始图和校正图

图12放射源到探测器120cm，旋转0度时，直尺的原始图和校正图

图中符号说明如下：

O₁-标记板的坐标原点；O-新坐标系的原点，位于标记板原点与Y方向单元小格的中间；X-X轴，Y-Y轴，Z-Z轴，垂直XO₁Y平面并指向里；

R——C型臂系统的旋转矩阵；T——C型臂系统的平移矩阵；f——焦距；u₀，v₀——C型臂X射线投影图像的中心；k₁，k₂——径向失真畸变参数；p₁，p₂——离心失真畸变参数；s₁，s₂——薄透镜失真畸变参数；

1-计算机；2-影像增强器；3-双目摄像机；4-支架；5-标记板；6-标定模板；7-病床；8-放射源；9-C型臂；10-下层校准靶；11-上层校准靶；12-标记点。

(五)具体实施方式

下面通过具体实施例结合所选用的具体硬件载体对本发明一种基于透视成像模型标定的C型臂图像校正方法，进行详细叙述：

第一种情况：在C型臂旋转角度为0度和放射源到探测器距离为124cm时。

(1)设计C型臂系统标定的标定模板。选择一种X射线不能穿透的材料制作一个双层的立体标定模板。本实施例采用的是密度相对轻的铝，以减小标定模板的自身重量对系统的影响。如图9，标定模板有上下两层，其中上层校准靶分布有直径为2mm的9个小孔，分布规律如图3所示。下层层校准靶分布有221个小孔，小孔相互间距为20mm。在两条相互垂直的直径上，离圆心分别为40mm和80mm的小孔直径为6mm，圆心的直径为6mm，其余的小孔直径都为2mm(如图2)。两层校准靶之间的间距为100mm。

(2)组装C型臂标定和图像校正系统。具体为：双目摄像机固定在支架上，支架放置在病床的一侧，标记板粘贴于标定模板的一侧，标定模板通过固定卡爪固定在C型臂影像增强的外壳上，且有标记板的一侧面向双目摄像机，C型臂采集的X射线投影图像和以双目摄像机坐标系为基准的模板上标记点的三维坐标输送到另一个隔离间的PC机进行处理。

(3)建立标定模板中圆形标记点在摄像机坐标系下的空间三维坐标；此步骤由双目摄像机定位标记板和不同坐标系之间的平移完成，具体为：

标记板的X轴垂直上层校准靶的圆形平面，Z轴指向上层校准靶的圆心，Y轴与上层校准靶的圆边相切，标记板的原点O₁为X轴，Y轴和Z轴的交点。设定新坐标系原点在标记板原点与Y方向单元小格的中间，且以这点为原点O建立坐标系，X轴垂直上层校准靶的圆形平面，方向与标记板X轴单位矢量平行，Z轴指向上层校准靶的圆心，Y轴方向与标记板Y轴单位矢量平行。原点O₁在摄像机的坐标系中为{31.74642，28.73861，675.62539}，Y方向中离原点一个单元格的坐标在摄像机的坐标系中为{39.57450，28.12457，666.93914}，X方向中离原点一个单元格的坐标在摄像机的坐标系中为{31.14742，5.00666，677.09865}，新坐标系原点O在摄像机的坐标系中为{35.6605，28.4316，671.2823}。原点O在摄像机的坐标系中为，X方向，Y方向和Z方向的单位矢量分别为 $\stackrel{\→}{x}=\{-0.0252,-\mathrm{0.9978,0.0619}\},$ $\stackrel{\→}{y}=\{0.6685,-0.0524,-0.7418\},$ $\stackrel{\→}{z}=\{-0.7434,-0.0227,-0.6684\}.$ 把原点O建立的坐标系沿Z轴平移到上层圆形平面的圆心，且在摄像机中坐标系的坐标为{202.9289，33.5453，821.6637}。以这个坐标系为基准确立上层标记模板上标记点圆心在摄像机坐标系下的空间位置。最后把圆心的坐标沿X轴平移到下层圆形模板的圆心，且在摄像机坐标系的坐标为{205.7495，145.2946，814.7264}。以这个坐标系为基准确定下层模板的标记点在摄像机坐标系下的空间坐标。

建立C型臂X射线标定模板投影图像中标记点的图像坐标和标定模板标中标记点的三维空间坐标之间的几何映射关系；此步骤由透视成像模型标定C型臂系统完成，具体为：

(4)基于透视成像模型的方法对C型臂系统进行标定。

①.提取C型臂标定模板X射线投影图像在图像中心区域的29个标记点，并用Tsai两步法对C型臂系统进行标定：第一步利用最小二乘法求解超定线性方程，给出外部参数：

$R=\left[\begin{array}{ccc}0.4615& -0.0448& -0.8860\\ 0.8904& 0.0165& 0.4548\\ -0.0057& -0.9988& 0.0476\end{array}\right]$

T_x＝634.8773，T_y＝-553.6835

第二步使用优化搜索求解内部参数：在使用优化搜索算法时，k₁，k₂的初值为0.(u₀，v₀)的初值为图像的中心(512.5，512.5)，T_z，f的初值为下式计算所得：

$[y_i-d_y(v_{\mathrm{di}}-v_0)]\·\left[\begin{array}{c}f\\ T_z\end{array}\right]=w_i{d}_y(v_{\mathrm{di}}-v_0)$

其中，y_i＝r₄x_wi+r₅y_wi+r₆z_wi+T_y，f＝1090.29mm，T_z＝1126.99。优化搜索之后，

      w_i＝r₇x_wi+r₈y_wi+r₉z_wi

f＝1019.31mm，T_z＝1080.08，u₀＝512.54，v₀＝512.20，k₁＝-9.6217×10^-6，k₂＝-8.0710×10^-6。

②.在整个图像中使用Levenberg-Marquadt迭代优化算法，精确求解C型臂系统的内外参数R，T和k₁，k₂，p₁，p₂，s₁，s₂。这里迭代的初始值设置为Tsai法得到的外参数值R，T和k₁，k₂，p₁＝0，p₂＝0，s₁＝0，s₂＝0，标记点是从图像中提取的所有的标记点。

$\min {\left|\right|\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(m_i-m(H,P,R,T,M_i))\left|\right|}^2$

这里初始值为R，T，K(k₁，k₂，p₁，p₂，s₁，s₂)计算的

$R=\left[\begin{array}{ccc}0.5176& -0.0485& -0.9825\\ 0.9608& 0.0191& 0.5125\\ -0.3130& -0.9793& 0.3255\end{array}\right],T=\left[\begin{array}{c}702.4840\\ -615.3309\\ 989.7990\end{array}\right]$ k₁＝-5.1760×10^-6，k₂＝-5.2563×10^-6

p₁＝8.8384×10^-5，p₂＝7.6464×10^-5，s₁＝-1.4860×10^-4，s₂＝-1.2559×10^-4。

③.在整个图像中(u₀，v₀)，f作为初始值，用Levenberg-Marquadt算法优化精确求解C型臂系统的(u₀，v₀)，f。

$\min {\left|\right|\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(m_i-m(H,P,R,T,M_i))\left|\right|}^2$ 计算的u₀＝512.5376，v₀＝512.2044，f＝1019.31mm.

4.循环上述步骤直至迭代次数完成或失真参数不再出现明显变化，这里最后得到k₁＝-5.1759×10^-6，k₂＝-5.2564×10^-6，p₁＝8.8386×10^-5，p₂＝7.6459×10^-5，s₁＝-1.4861×10^-4，s₂＝-1.2559×10^-4。f＝1019.31mm，u₀＝512.5439，v₀＝512.2055。

$R=\left[\begin{array}{ccc}0.5202& -0.0488& -0.9875\\ 0.9657& 0.0192& 0.5151\\ -0.3146& -0.9843& 0.3272\end{array}\right],T=\left[\begin{array}{c}706.0644\\ -618.4685\\ 994.8177\end{array}\right].$

(5)建立C型臂X射线投影失真图像的图像坐标和理想图像的图像坐标之间的映射关系；此步骤由从透视成像模型标定C型臂系统中获得的失真参数对失真图像进行图像校正完成；也就是根据图像成像的理想位置与图像成像失真偏移的实际位置关系，把标定获得的失真畸变参数代入其中，获得校正图像的理想位置，最后进行灰度映射得到校正图像。图像成像的理想位置与图像成像失真偏移的实际位置关系，

x′＝x+δ_x，y′＝y+δ_y，δ_x(x，y)＝k₁x(x²+y²)，δ_y(x，y)＝k₂y(x²+y²)

其中：(x，y)为图像成像的理想位置；

(x′，y′)为图像成像失真偏移的实际位置；

δ_x，δ_y为非线性畸变值；

k₁，k₂为径向失真畸变参数；

δ_x(x，y)＝k₁x(x²+y²)+(p₁(3x²+y²)+2p₂xy)+s₁(x²+y²)，δ_y(x，y)＝k₂y(x²+y²)+(p₂(3y²+x²)+2p₁xy)+s₂(x²+y²)

如果所得到的校正图像的理想位置不是整数，那需要采用插值算法对图像进行灰度赋值。这里采用的是最近邻插值算法，关系如下：

f(x)＝f(x_k)，如果1/2(x_k-1+x_k)＜x＜1/2(x_k+x_k+1)

校正的标定模板图像如图10.

第二种情况：在C型臂旋转角度为30度和放射源到探测器距离为124cm时。

(1)设计C型臂系统标定的标定模板。选择一种X射线不能穿透的材料制作一个双层的立体标定模板。本实施例采用的是密度相对轻的铝，以减小标定模板的自身重量对系统的影响。如图8，标定模板有上下两层，其中上层校准靶分布有直径为2mm的9个小孔，分布规律如图3所示。下层层校准靶分布有221个小孔，小孔相互间距为20mm。在两条相互垂直的直径上，离圆心分别为40mm和80mm的小孔直径为6mm，圆心的直径为6mm，其余的小孔直径都为2mm(如图2)。两层校准靶之间的间距为100mm。

(2)组装C型臂标定和图像校正系统。具体为：双目摄像机固定在支架上，支架放置在病床的一侧，标记板粘贴于标定模板的一侧，标定模板通过固定卡爪固定在C型臂影像增强的外壳上，且有标记板的一侧面向双目摄像机，C型臂采集的X射线投影图像和以双目摄像机坐标系为基准的模板上标记点的三维坐标输送到另一个隔离间的PC机进行处理。

(3)建立标定模板中圆形标记点在摄像机坐标系下的空间三维坐标；此步骤由双目摄像机定位标记板和不同坐标系之间的平移完成，具体为：

标记板的X轴垂直上层校准靶的圆形平面，Z轴指向上层校准靶的圆心，Y轴与上层校准靶的圆边相切，标记板的原点O₁为X轴，Y轴和Z轴的交点。设定新坐标系原点在标记板原点与Y方向单元小格的中间，且以这点为原点O建立坐标系，X轴垂直上层校准靶的圆形平面，方向与标记板X轴单位矢量平行，Z轴指向上层校准靶的圆心，Y轴方向与标记板Y轴单位矢量平行。原点O₁在摄像机的坐标系中为{206.57662，41.27837，560.99809}，Y方向中离原点一个单元格的坐标在摄像机的坐标系中为{213.35762，46.37188，552.59965}，X方向中离原点一个单元格的坐标在摄像机的坐标系中为{215.92559，19.99725，555.73404}，新坐标系原点O在摄像机的坐标系中为{209.9671，43.8251，556.7989}。原点O在摄像机的坐标系中为，X方向，Y方向和Z方向的单位矢量分别为 $\stackrel{\→}{x}=\{0.3923,-0.8929,-0.2209\},$ $\stackrel{\→}{y}=\{\mathrm{0.5681,0.4267},-0.7036\},$ $\stackrel{\→}{z}=\{-0.7226,-0.1505,-0.6747\}.$ 把原点O建立的坐标系沿Z轴平移到上层圆形平面的圆心，且在摄像机中坐标系的坐标为{372.5446，77.6956，708.6078}。以这个坐标系为基准确立上层标记模板上标记点圆心在摄像机坐标系下的空间位置。最后把圆心的坐标沿X轴平移到下层圆形模板的圆心，且在摄像机坐标系的坐标为{328.6099，177.7045，733.3458}。以这个坐标系为基准确定下层模板的标记点在摄像机坐标系下的空间坐标。

建立C型臂X射线标定模板投影图像中标记点的图像坐标和标定模板标中标记点的三维空间坐标之间的几何映射关系；此步骤由透视成像模型标定C型臂系统完成，具体为：

(4)基于透视成像模型的方法对C型臂系统进行标定。

①.提取C型臂标定模板X射线投影图像在图像中心区域的29个标记点，并用Tsai两步法对C型臂系统进行标定：第一步利用最小二乘法求解超定线性方程，给出外部参数：

$R=\left[\begin{array}{ccc}0.3733& 0.3980& -0.8380\\ 0.8352& 0.2589& 0.4852\\ 0.4101& -0.8810& -0.2358\end{array}\right]$

T_x＝422.4281，T_y＝-674.4778

第二步使用优化搜索求解内部参数：在使用优化搜索算法时，k₁，k₂的初值为0.(u₀，v₀)的初值为图像的中心(512.5，512.5)，T_z，f的初值为下式计算所得：

$[y_i-d_y(v_{\mathrm{di}}-v_0)]\·\left[\begin{array}{c}f\\ T_z\end{array}\right]=w_i{d}_y(v_{\mathrm{di}}-v_0)$

其中，y_i＝r₄x_wi+r₅y_wi+r₆z_wi+T_y，f＝1076.1798mm，T_z＝1200.7733。优化搜索之后，

      w_i＝r₇x_wi+r₈y_wi+r₉z_wi

f＝1011.82mm，T_z＝1159.9592，u₀＝512.58，v₀＝512.04，k₁＝-9.7856×10^-6，k₂＝-8.1368×10^-6。

②.在整个图像中使用Levenberg-Marquadt迭代优化算法，精确求解C型臂系统的内外参数R，T和k₁，k₂，p₁，p₂，s₁，s₂。这里迭代的初始值设置为Tsai法得到的外参数值R，T和k₁，k₂，p₁＝0，p₂＝0，s₁＝0，s₂＝0，标记点是从图像中提取的所有的标记点。

$\min {\left|\right|\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(m_i-m(H,P,R,T,M_i))\left|\right|}^2$

这里初始值为R，T，K(k₁，k₂，p₁，p₂，s₁，s₂)计算的

$R=\left[\begin{array}{ccc}0.4012& -0.4292& -0.9054\\ 0.8773& 0.2691& 0.5243\\ 0.0933& -1.0592& 0.0647\end{array}\right],T=\left[\begin{array}{c}456.9692\\ -718.6952\\ 1120.4216\end{array}\right]$ k₁＝-5.1328×10^-6，k₂＝-5.2965×10^-6

p₁＝1.1022×10^-4，p₂＝9.3400×10^-5，s₁＝-1.7920×10^-4，s₂＝-1.7742×10^-4。

③.在整个图像中(u₀，v₀)，f作为初始值，用Levenberg-Marquadt算法优化精确求解C型臂系统的(u₀，v₀)，f。

$\min {\left|\right|\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(m_i-m(H,P,R,T,M_i))\left|\right|}^2$ 计算的u₀＝512.5806，v₀＝512.0445，f＝1011.82mm

④.循环上述步骤直至迭代次数完成或失真参数不再出现明显变化，这里最后得到k₁＝-5.1327×10^-6，k₂＝-5.2965×10^-6，p₁＝1.1022×10^-4，p₂＝9.3390×10^-5，s₁＝-1.7922×10^-4，s₂＝-1.7741×10^-4。f＝1011.82mm，u₀＝512.5892，v₀＝512.0748。

$R=\left[\begin{array}{ccc}0.4011& -0.4291& -0.9051\\ 0.8770& 0.2690& 0.5242\\ 0.0932& -1.0589& 0.0647\end{array}\right],T=\left[\begin{array}{c}456.8393\\ -718.4960\\ 1120.4216\end{array}\right].$

(5)建立C型臂X射线投影失真图像的图像坐标和理想图像的图像坐标之间的映射关系；此步骤由从透视成像模型标定C型臂系统中获得的失真参数对失真图像进行图像校正完成；也就是根据图像成像的理想位置与图像成像失真偏移的实际位置关系，把标定获得的失真畸变参数代入其中，获得校正图像的理想位置，最后进行灰度映射得到校正图像。图像成像的理想位置与图像成像失真偏移的实际位置关系，

x′＝x+δ_x，y′＝y+δ_y，δ_x(x，y)＝k₁x(x²+y²)，δ_y(x，y)＝k₂y(x²+y²)

其中：(x，y)为图像成像的理想位置；

(x′，y′)为图像成像失真偏移的实际位置；

δ_x，δ_y为非线性畸变值；

k₁，k₂为径向失真畸变参数；

δ_x(x，y)＝k₁x(x²+y²)+(p₁(3x²+y²)+2p₂xy)+s₁(x²+y²)，δ_y(x，y)＝k₂y(x²+y²)+(p₂(3y²+x²)+2p₁xy)+s₂(x²+y²)

如果所得到的校正图像的理想位置不是整数，那需要采用插值算法对图像进行灰度赋值。这里采用的是最近邻插值算法，关系如下：

f(x)＝f(x_k)，如果1/2(x_k-1+x_k)＜x＜1/2(x_k+x_k+1)

校正的标定模板图像如图11.

在C型臂旋转角度为0度和放射源到探测器距离为120cm时，使用标定C型臂系统获得的参数校正X射线投影直尺失真图，如图12。

Claims (1)

1、一种基于透视成像模型标定的C型臂图像校正方法，它是在以下硬件载体上实现的：PC机、双目摄像机、支架、C型臂系统、标记板、标定模板以及固定卡爪；双目摄像机固定在支架上，支架放置在病床的一侧，标记板粘贴于标定模板的一侧，标定模板通过固定卡爪固定在C型臂系统影像增强器的外壳上，且有标记板的一侧面向双目摄像机，C型臂系统采集的X射线投影图像和以双目摄像机坐标系为基准的模板上标记点的三维坐标输送到另一个隔离间的PC机进行处理；其特征在于：该方法具体步骤如下：

步骤一：设计C型臂系统标定的标定模板具体为：测量C型臂影像增强器外壳直径大小，以C型臂影像增强器外壳直径451mm为直径制作一个C型臂罩；C型臂罩上安装有两层校准靶，其间距为100mm；下层校准靶直径为380mm，且在直径为350mm的圆形区域内，有以间距20mm均匀规则分布的225个圆形小孔；其中圆心和两个相互垂直的直径上有四个小孔的直径为6mm，其余直径为2mm；上层校准靶直径为380mm，且在直径为350mm的圆形区域内分布有9个直径为2mm的圆形小孔，其中圆心小孔垂直对应上层校准靶圆心孔，其余8个小孔分成两组，每组相距圆心距离分别是和

步骤二：获取标记点的失真坐标在失真图像上通过图像处理获得，具体为：

①「.对标定模板X射线投影图像使用中值滤波去噪平滑处理；

②.使用形态学中的开运算提取标定模板背景图像；

③.提取标定模板X射线投影图像中的标记点；用标定模板X射线投影图像减去标定模板背景图像得到只含圆形标记点的图像，随后经过二值化，开运算和闭运算去除噪声得到只含圆形标记点的二值图像；

④.计算出这些圆形标记点的质心作为标记模板上标记点在X射线投影图像中的图像坐标；

步骤三：C型臂系统标定包括标定模板上标记点的三维空间坐标、透视成像模型、畸变模型和失真坐标通过系统标定获得透视成像参数，具体为：

i)标定模板上标记点的三维空间坐标；此步骤由双目摄像机定位标记板和不同坐标系之间的平移完成，具体为：

①「双目摄像机固定在支架上，支架放置在病床的一侧，标记板粘贴于标定模板的一侧，标定模板通过固定卡爪安装在C型臂影像增强器的外壳上，且有标记板的一侧面向双目摄像机；

②标记板的X轴垂直上层校准靶的圆形平面，Z轴指向上层校准靶的圆心，Y轴与上层校准靶的圆边相切，标记板的原点O₁为X轴，Y轴和Z轴的交点；设定新坐标系原点在标记板原点与Y方向单元小格的中间，且以这点为原点O建立坐标系，X轴垂直上层校准靶的圆形平面，方向与标记板X轴单位矢量平行，Z轴指向上层校准靶的圆心，Y轴方向与标记板Y轴单位矢量平行，原点O₁在摄像机的坐标系中为(X，Y，Z)，Y方向中离原点一个单元格的坐标在摄像机的坐标系中为(X₂，Y₂，Z₂)，X方向中离原点一个单元格的坐标在摄像机的坐标系中为(X₁，Y₁，Z₁)，新坐标系原点O在摄像机的坐标系中为((X₂-X)/2，(Y₂-Y)/2，(Z₂-Z)/2)；把新坐标系移到上层校准靶的圆心，且以这个坐标系为参考，逐步确定上层校准靶中小孔的圆心在摄像机坐标系下的三维空间坐标，然后再把新坐标系移到下层校准靶的圆心，且以这个坐标系为参考，逐步确定下层校准靶中小孔的圆心在摄像机坐标系下的三维空间坐标；关系式如下：

(X₀₁，Y₀₁，Z₀₁)＝(X_标，原+X_摄，Y_标，原+Y_摄，Z_标，原+Z_摄)    (式8)

其中：(X₀₁，Y₀₁，Z₀₁)是标定板原点在摄像机坐标系下的坐标；

(X_摄，Y_摄，Z_摄)是标定板原点在摄像机坐标系下的偏移量；

假设的原点O：

(X₀，Y₀，Z₀)＝(X₀₁+X₂，Y₀₁+Y₂，Z₀₁+Z₂)/2    (式9)

X，Y方向的单位矢量：

$\stackrel{\→}{x}=\frac{(X_1-X_{01},Y_1-Y_{01},Z_1-Z_{01})}{\sqrt{{(X_1-X_{01})}^2+{(Y_1-Y_{01})}^2+{(Z_1-Z_{01})}^2}}$ (式10)

$\stackrel{\→}{y}=\frac{(X_2-X_{01},Y_2-Y_{01},Z_2-Z_{01})}{\sqrt{{(X_2-X_{01})}^2+{(Y_2-Y_{01})}^2+{(Z_2-Z_{01})}^2}}$

(式11)

坐标系平移到标定模板上层平面的圆心：

$\stackrel{\→}{z}=\stackrel{\→}{y}\×\stackrel{\→}{x}$ (式12)

$({X_0}^{\′},{Y_0}^{\′},{Z_0}^{\′})=(X_0,Y_0,Z_0)-\stackrel{\→}{z}\×225$ (式13)

其中：(X₀′，Y₀′，Z₀′)是坐标系在标定模板上层平面圆心处的原点

$({X_0}^{\′\′},{Y_0}^{\′\′},{Z_0}^{\′\′})=({X_0}^{\′},{Y_0}^{\′},{Z_0}^{\′})-\stackrel{\→}{x}\×112$ (式14)

其中：(X₀″，Y₀″，Z₀″)是坐标系在标定模板下层平面圆心处的原点

ii)透视成像模型标定，畸变模型和失真坐标通过系统标定获得透视成像参数，具体步骤为：

①.提取C型臂标定模板X射线投影图像在图像中心区域的标记点，并用Tsai两步法进行标定；该方法有两步：第一步利用最小二乘法求解超定线性方程，给出外部参数；第二步使用优化搜索求解内部参数；摄像机坐标系与世界坐标系关系式表示矢量形式：

[x_wY_d y_wY_d z_wY_d Y_d -x_wX_d -y_wX_d -z_wX_d]·[r₁/T_y r₂/T_y r₃/T_y T_x/T_y r₄/T_y r₅/T_y r₆/T_y]′＝X_d

其中：x_w，y_w，z_w世界坐标；

x_d，y_d，z_d摄像机坐标‘

r_i(i＝1，2，3，4，5，6)旋转矩阵分量；

T_y为Y方向的平移分量；(式1)

实际图像坐标到摄像机坐标的变换：

u_d＝s_xd_x′^-1X_d+u₀，v_d＝d_y ^-1Y_d+v₀    (式2)

其中：d_x为摄像机在X方向的像素间距；

d_y为摄像机在Y方向的像素间距；

s_x为图像尺度因子；

(u₀，v₀)为光学中心；

图像成像的理想位置与图像成像失真偏移的实际位置关系：

x′＝x+δ_x，y′＝y+δ_y，δ_x(x，y)＝k₁x(x²+y²)，δ_y(x，y)＝k₂y(x²+y²)    (式3)

其中：(x，y)为图像成像的理想位置；

(x′，y′)为图像成像失真偏移的实际位置；

δ_x，δ_y为非线性畸变值；

k₁，k₂为径向失真畸变参数；

②.在整个图像中使用Levenberg-Marquadt迭代优化算法，精确求解C型臂系统的内外参数R，T和k₁，k₂，p₁，p₂，s₁，s₂；这里迭代的初始值设置为Tsai法得到的外参数值R，T和k₁，k₂，p₁＝0，p₂＝0，s₁＝0，s₂＝0，标记点是从图像中提取的所有的标记点；图像成像失真的偏移量：

δ_x(x，y)＝k₁x(x²+y²)+(p₁(3x²+y²)+2p₂xy)+s₁(x²+y²)，δ_y(x，y)＝k₂y(x²+y²)+(p₂(3y²+x²)+2p₁xy)+s₂(x²+y²)

其中：p₁，p₂为切向失真畸变参数；

s₁，s₂为薄透镜失真畸变参数；    (式4)

实际真实坐标(x′，y′)和理想真实坐标(x，y)的关系改成向量相乘的形式：

AP＝X′-X    (式5)

其中： $A=\left[\begin{array}{cccccc}x(x^2+y^2)& 0& 3x^2+y^2& 2\mathrm{xy}& x^2+y^2& 0\\ 0& y(x^2+y^2)& 2\mathrm{xy}& x^2+3y^2& 0& x^2+y^2\end{array}\right]$

P＝[k₁ k₂ p₁ p₂ s₁ s₂]′，X′＝[x′y′]′，X＝[x  y]′

用Levenberg-Marquadt算法优化精确求解，就是求目标函数：

$\min {\left|\right|\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(m_i-m(H,P,R,T,M_i))\left|\right|}^2$ (式6)

假设模板图像上有n个标记点，其中m_i是第i个像素点图像坐标，M_i是第i个像素点对应的三维空间坐标，H内部参数矩阵；

③再次在整个图像中使用Levenberg-Marquadt迭代优化算法，精确求解C型臂系统的内参数f，u₀，v₀；迭代的初始值设置为Tsai法得到的内参数值f，u₀，v₀；

④循环迭代这个步骤直至迭代次数完成或达到误差限而终止程序运行且输出结果；

步骤四：图像校正由标定的透视成像模型，对失真图像进行校正；具体步骤是用从透视成像模型标定C型臂系统中获得的失真参数对失真图像进行图像校正，也就是根据图像成像的理想位置与图像成像失真偏移的实际位置关系，把标定获得的失真畸变参数代入其中，获得校正图像的理想位置，最后进行灰度映射得到校正图像；如果所得到的校正图像的理想位置不是整数，那需要采用插值算法对图像进行灰度赋值；这里采用的是最近邻插值算法，关系如下：

f(x)＝f(x_k)，如果1/2(x_k-1+x_k)＜x＜1/2(x_k+x_k+1)    (式7)。

Images