CN103827917A - 用于使用棋盘图案的一幅或多幅图像的自动相机校准的方法和装置 - Google Patents

Abstract

本发明涉及一种为了最佳图像质量而校准并确定镜头围绕其对称轴的旋转范围、由此更新投影模型并使用并行处理实时校正图像的径向失真的高精度方法、模型和装置。在此提供的解决方案依赖于光学设备的完全几何校准，并在随后实时呈现透视校正图像，该光学设备诸如通常用于医疗和一般工业的相机。校准包括确定将每个参数分配至对应入射光的3D方向的合适映射函数的参数。这类解决方案的实际实现是非常直接的，仅要求相机捕捉轻易可用的校准目标的单个视图，并且计算机实现的处理流水线使用计算平台的并行执行能力实时运行，该校准目标可装配在专用校准装置内部。

Description

用于使用棋盘图案的一幅或多幅图像的自动相机校准的方法和装置

技术领域

本发明涉及一种为了最佳图像质量而校准并确定镜头围绕其对称轴的旋转范围、由此更新投影模型并使用并行处理实时校正图像的径向失真的高精度方法、模型和装置。

背景技术

通常需要光学仪器产生给定对象的几何组成图像，其中3D中的每个点被投影到该图像中的点内。该图像通常根据某些预定义成像模型（在此情况下假设为针孔模型）形成。实际光学***与理想性态的偏差导致在所得图像中引入像差。在此像差是会出现在彩色和灰度级成像设备两者中的镜头失真。其特性在考虑定桶形或枕形失真效应的情况下主要是像素的径向几何移位。镜头失真是许多相对直接的图像分析任务中的棘手问题。它完全劣化了实际图像中所得测量值的精度，在该实际图像中正常透视相机中的像素会沿径向遭受较大位移。此外，它还会影响视觉感知，即通过改变深度暗示。

解决镜头失真的一种已知方法是结合相机***的完全校准。完全相机校准实质上意味着获取相机的内部和外部参数[8]。许多应用不获取完整的内部相机参数以及相对于某一校准目标的相对定向。非度量（Non-metric）应用的不同之处在于在无需估计整个相机的内部和外部参数的情况下确定失真。

在1960年代，摄影测量领域内工作的Duane C.Brown提出了一种基于直线必须被成像为直线的常识来确定镜头失真的方法。由Brown在1971年发布并由Fryer和Brown在1986年扩展的这一技术变为已知的“铅垂线”法。完整的历史性回顾由Clarke和Fryer在1998年给出。该技术由提供铅垂线法简化版本的机器视觉社区所采用，诸如Prescott和McLean（1997b）（同样参见专利Prescott和McLean（1997a）XP000683415），Haneishi等人（1995a）（同样参见专利Haneishi等人（1995b）XP000527216），Poulo和Gorman（1999）US6002525以及Asari等人（1999），这些文献都描述了使用同线点图像校正失真的类似常识。由于这些方法仅估计失真，因此它们有时被称为非度量校准。

基于“铅垂线”校准的一个固有问题在于优化/搜索必须相对于直线参数（未知的）和失真参数（也是未知的）进行。也可以利用例如Devernay和Faugeras（2001）、Tzu-Hung（2003）以及Bing（1999）EP0895189的替换方法，其中迭代调整失真参数以最小化拟合至失真线分段的线。由于目标误差和失真参数之间不存在足够有效的数学关系，因此没有分析推导可用。这导致收敛缓慢并且除非如Swaminathan和Nayar（2000）所述那样采取特定步骤，否则对升高的失真水平会变得不稳定。在此非度量方法中，Swaminathan和Nayar（2000）在失真空间而非平常的非失真空间内重新阐述了目标函数。

另一种方法已由Ahmed和Farag（2001）提出，其中使用检测线的曲率来估计导出的失真式的参数。然而，模拟结果示出了存在噪声的糟糕性能，并且实际结果缺乏量化估计。

校准失真的一种更为标准的方法是同时估计相机的外源和本征参数。Tsai的方法（Tsai，1987）涉及在考虑到一组控制点的3D位置的情况下通过迭代数字优化方案同时估计单失真参数和某些内部参数，诸如焦距。这一方法的缺点在于它需要知晓3D控制点，但反过来对所有但简单的失真轮廓线仅提供相对低的精度。这一原理的算术变化已由若干作者提出，诸如针对镜头失真使用更为合适模型的Weng等人（1992）以及Wei和Ma（1994）。这些方法也需要知晓3D控制点。

失真校正图像的渲染在描述了需要3D控制点或多图像组的2D控制点的Heikkila和Silven（1997），wHeikkila（2000）中得到了深入研究。基于多组2D控制点的替换方法也由Zhang（1998，2000）以及Sturm和Maybank（1999）推进。该技术通过替换的线性最小二乘解（随后在包括所有估计参数的数字最小化中迭代调整）来解决失真。这些技术的相对复杂性由于包括了镜头失真而增加。

另一方面，存在其中仅需要移除失真而非本征和外源参数整体的许多情况。一个例子在于实际图像中的多视图几何的估计，其中已经具体提出了顺应镜头失真的技术。Zhang（1996）深入研究了同时估计失真参数和基本矩阵的可能性。结果是如果噪声较低且失真较高则是可能的。Fitzgibbon（2001）（具有专利Fitzgibbon（2003）GB2380887）、Micusik和Pajdla（2003）以及Barreto和Daniilidis（2004）使用会导致多项式本征值问题的替换失真模型以及更为可靠的失真和几何估计。Stein（1997）采取了相反的方法，并且使用基本矩阵估计中的误差作为目标误差来估计失真参数。

存在失真校准的替换方法，其中由于失真自由场景的取代而丢弃控制点对应。这些场景随后由相机***成像，随后进行图像对准处理以校正失真。Lucchese和Mitra（2003）描述了一种其中失真图像被扭曲直到其（关于强度）与基准图像配准为止。使用粗略滤波以找出配准的类似技术在Tamaki（2002）（具有专利Tamaki等人（2002）US2002057345）中描述，而Sawhney和Kumar（1999）（同样具有专利Kumar等人（1998）WO9821690）则描述了一种不需要非失真基准图像的配准技术。代替地，多个图像与马赛克图像（例如由若干不同视图组合得到的全景图像）的生成配准，并同时估计失真。这些技术基于极高的计算开销，在Tamaki（2002）的引用中为20分钟。

最后一类非度量校准方法基于的是频域中失真引起的高阶相关。Farid和Popescu（2001）描述了一种技术，但是相比于常规相机校准技术其性能不佳，并且看上去还多少依赖于图像内容。Yu（2004）使用替换的失真模型进一步发展了这种方法，并且报告了在源图像是常规校准目标的情况下使用常规相机校准完成的精度方法。最后，由Gary（2003）US2003141433提出了一种建立弯曲CCD阵列的途径。镜头失真轮廓线由一系列线段复制到阵列内，由此得到的图像看上去没有失真。

发明内容

在此提供的解决方案依赖于光学设备（诸如通常用于医疗和一般工业的相机）的完全几何校准，并在随后实时呈现透视校正图像。校准包括确定将每个参数分配至对应入射光的3D方向的合适映射函数的参数。这类解决方案的实际实现是非常直接的，仅要求相机捕捉容易获得的校准目标（可装配到专用校准装置内部）的单个视图，并且计算机实现的处理流水线使用计算平台的并行执行能力实时运行。

本发明还涉及在计算机介质上具体化的计算机程序，其在由处理器执行时实现所有的方法步骤。

对本发明更为全面的理解将从随后的描述、附图和权利要求中得到。

附图说明

被包括用于提供对本发明进一步理解的附图并入并组成本应用一部分，所述附图例示了本发明的实施例并且与描述一起用于阐明本发明的原理。在附图中：

图1（a）示出了组成管道镜的具有CCD相机的刚性内窥镜。该***能够对难以接近或受限的人体体腔进行可视化，并被广泛用于手术和诊断。图1（b）是通过30°倾斜视角关节内窥镜获取的膝关节内部图像。

图2示出了根据平面棋盘图案的单幅图像校准内窥镜相机。图2（a）的校准图像使用图2（b）所示的设置以非限制性的方式（没有关于相机和平面的相对姿态的假设）获取。呈现的设置由控制光照条件的丙烯酸盒组成。平面栅格是在无需用户介入的情况下使得能够鲁棒且准确检测图像边角的的背光。

图3是示出了被提出用于校准相机并校正图像/视频的径向失真的***的不同模块的方案。左手侧涉及仅在将内窥镜镜头连接至CCD相机之后执行一次并要求使用图2中的设置获取单幅校准图像的初始化程序。右手侧则示出了在每个帧时刻i执行以校正内窥镜视频流中的径向失真的处理流水线。

图4示出了边界轮廓的跟踪。从左到右的方案涉及在时刻i获取的原始图像、使用仿射变换S的扭曲图像、以及通过从笛卡尔坐标改变至球形坐标而获取的极性图。红色虚线覆盖图涉及在先边界估计结果，而十字则是当前轮廓S映射Ω_i-1至单位元中的检测点。紫色阴影区是遭受极性扭曲并且在其中执行搜索的图像区。

图5示出了（a）在镜头经历完全旋转时主点c、三角标记p和边界中心w的绘图。（b）和（c）则例示了被考虑用于更新内窥镜相机的本征矩阵的模型和假设。镜头将虚拟图像投影至由CCD相机成像到I内的平面I′上。镜头围绕平行于光轴且与平面I和I′正交的轴l进行相对旋转。c表示主点而q则是导出图像旋转的中心。

图6示出了（a）在三角标记被正确检测时以及（b）没有标记信息时的图像旋转中心q_i的计算。

图7示出了作为三角位移α的函数的RMS再投影误差。蓝线表示没有旋转补偿（其值是图像边角的RMS位移）情况下的误差。红线则示出了使用本发明模型（式11）时的再投影误差。对于α＝0，描绘的值是SIC校准的再投影误差。

图8示出了具有镜头探针旋转的内窥镜视频序列的径向失真校正。原始的和扭曲的帧分别在顶行和底行中呈现。（a）示出了用于执行初始校准的基准位置（α＝0）。（b）比较了在没有（左）和具有（右）镜头旋转运动补偿的情况下的失真校正结果。（c）和（d）则是示出了就感知场景而言的RD校正的益处。由于选择维持中央图像区中的尺寸，扭曲的帧更大。

图9示出了用于不同输出图像分辨率的每帧执行时间。输入图像具有1600×1200像素的分辨率，而输出图像具有可变方形大小。右侧的表格示出了相对于单纯CPU解决方案使用优化的混合CPU+GPU方法实现的加速。

图10示出了用于混合CPU+GPU计算单元的特定情况的径向失真校正算法步骤。红框标识CUDA内核调用。绿点框表示分配的设备存储器。除了在CPU上计算的紫框之外，所有处理都在GPU上执行。

图11示出了在SIC和bouguet校准工具盒之间的校准参数比较。

图12示出了针对不同图像分辨率使用全***运行优化的混合CPU+GPU解决方案所实现的执行时间（毫秒单位）和每秒帧数（fps）。

图13示出了根据单幅图像（软件屏幕快照）的校准。

图14示出了校准盒方案。（1）–使用不透明的防水材料建造以避免光从外部进入的盒上隔室；（2）–用于镜头探针进入的开口；（3）–最小化从外部进入的光的不透明材料膜；（4）–图像获取设备；（5）–其中安装了漫射光源的盒的下隔室；（6）–上下隔室之间的隔离物；（7）–由透明或不透明材料制成且被置于隔离物顶部上的已知图案的平面网格。

具体实施方式

本发明提出了用于镜头旋转情况下投影参数的校准和在线更新以及实时径向失真校正的完全解决方案。该解决方案包括图3的方案和图14的装置示出的模块和框，其中操作者从使用图2的设置获取棋盘图案的单幅图像开始。

在此呈现的***最初旨在用于典型地利用安装在标准CCD相机上的管道镜镜头探针的医疗内窥镜。然而，本发明的解决方案可以扩展至具有可以从精确几何校准和实时图像校正中获益的显著径向失真的相机的其他应用领域。这些领域的例子包括通常依赖于由安装了相机的鱼眼镜头、广角镜头、迷你镜头或低质量光学器件获取的图像/视频的监视、工业检查、汽车、机器人导航等。

在本专利申请中公开的装置包括用于获取已知图案的单幅图像的校准盒（图14）。获取的图像被用作新颖校准算法（图13）的输入，该算法将提供执行实时和独立于镜头旋转的图像失真校正所需的相机参数。

校准图像使用由不透明材料建造以避免光从外部进入的盒来获取（图13）。该盒具有用于镜头探针进入的一个或多个开口（图13-3），并且每个开口由用于最小化从外部进入的光的橡胶膜（图13-2）所围绕。该盒由透明或半透明材料的隔离物（图13-6）分隔成两个隔室。由透明或半透明材料制成的具有已知图案的平面网格（图13-7）被置于隔离物上的上隔室（图13-1）内。用于平面网格的散射背光照明的光源则被置于下隔室（图13-5）内。上隔室可由用于更好地重复将执行视觉检查的部件的流体填充。该装置包括具有顺序和并行处理能力且与用于图像获取的相机（图13-4）相连的计算机单元，诸如中央处理单元和图形处理单元。还需要连接至该计算机单元的显示器和接口设备用于可视化和输入用户命令。

用于在镜头旋转的情况下投影参数的校准和在线更新以及使用并行处理的实时图像修正/校正的方法总结如下：检测校准帧中的边角，并且在无需进一步用户介入的情况下估计相机本征K₀和径向失真ξ。在该简要的初始化步骤之后，对每幅获取的图像执行图3右侧所示的处理流水线。在每个帧时刻i，检测边界轮廓Ω_i以及三角标记p_i的位置。检测结果用作考虑到校准帧中的轮廓Ω_i以及三角标记p_i时估计滤波器中的输入，并在随后估计由于镜头探针相对于相机头的可能运动引起的相对图像旋转。这一2D旋转由基于新的自适应投影模型用作更新相机校准的输入的角度αi和固定点qi参数化。最后，当前几何校准K_i，ξ被用来扭曲输入帧并校正径向失真。类似于图像扭曲和边界检测，该处理流水线以计算密集任务实时运行，其使用GPU的并行执行能力而被有效执行。

A.一幅或多幅棋盘图案图像的获取

用户介入被限制为使用图13中详示的图2（b）的校准盒获取一幅或多幅校准图像。相机相对于棋盘平面的姿势是任意的，并且图13的盒旨在控制用于图像边角自动检测的照明条件。获取的图像与图2（a）类似。

B.校准图像中网格边角的检测以及图像平面对应的建立

重要的是意识到在具有严重径向失真的图像中定位边角会由于直线的弯曲而很成问题。在帧获取之后的处理步骤如下：

1）局域化其中失真更不明显的中央图像区域内12个边角的最小值。这一操作通过使用标准图像处理技术的探试算法进行。图像边角x和网格点g之间的一致对应通过计数棋盘图案的方形来实现。

2）图像平面对应被用于估计将平面网格坐标中的点映射至图像坐标中的点的函数，即使用类DLT方法估计

（B部分中的式7）。

3）使在图像***中生成边角假设（hypotheses）的单应性（homography）来将棋盘图案投影至图像平面。这些假设通过应用标准图像边角查找器来确认和细化。

4）迭代执行步骤2）和3）直到认为建立的对应的数目令人满意为止。

C.使用至少一幅校准图像内的图像平面对应确定相机本征参数

图3所示初始化程序的目的在于确定镜头探针处于基准位置Ω₀，p₀时的本征校准矩阵K₀和径向失真ξ。相机校准是文献中有几百篇引用文献的充分研究的主题。使用最广泛的软件或许是利用Zhang的方法[17]以根据3幅平面网格图像的最小值校准一般性相机的Bouguet的工具盒[16]。不幸的是，Bouguet工具盒不符合本申请的可用性要求。在实践中，用于实现精确结果的输入图像数远大于3，并且使用RD检测图像中的网格边角需要实质性的用户介入。本发明的方法可以使用最少一幅图像来工作。

若干作者已经解决了医疗内窥镜本征校准的特定问题[9]–[11]，[18]，[19]。然而，这些方法要么无法在OR中实际使用，要么它们利用圆点图案来使得能够进行校准点的自动检测，这降低可结果精度[20]。

装备有管道镜的相机是具有复杂光学布置的复合***。投影在中央[5]，并且图像失真由所谓的分隔模型良好描述[21]，[14]，[13]。

Barreto等人示出了按照在上述投影模型的相机能够根据从任意位置获取的平面棋盘图案的单幅图像来校准。设g是平面同质坐标中表达的棋盘中的点，x是图像平面中的对应点，且H是平面和相机之间相对姿态的单应性编码[8]。棋盘平面内的点g通过如下函数映射到图像内：

x～K₀Γ_ξ(Hg).

其中Γ_ξ是计及图像径向失真的非线性投影函数[21]，而K₀是具有如下结构的本征参数矩阵：

$K_0~\left(\begin{array}{ccc}\mathrm{af}& \mathrm{sf}& c_x\\ 0& a^{-1}f& c_y\\ 0& 0& 1\end{array}\right)---\left(3\right)$

其中f、a和s各自表示焦距、纵横比和偏斜，而c＝(c_x，c_y)^T则是图像主点的非同质坐标。

设

和是根据二阶Veronese映射[13]而具有对应于点g和x的提升表示(lifted representation)的6维同质矢量。其可被证明为

$\hat{x}~{\hat{H}}_{6\×6}\hat{g}---\left(7\right)$

其中

是6×6矩阵。由于每个图像平面对应对提升的单应性施加3个线性限制，因此可以使用类似直接线性变换(DLT)的方法而从12个点对应(x，g)中估计矩阵

考虑到则本征参数矩阵K₀、失真参数ξ和原始单应性H能够被直接因式分解。该初始相机校准可以使用标准迭代非线性优化而被进一步细化。

D.圆形和帧区域之间的边界轮廓检测

找出分割圆形和帧区域的图像轮廓(参见图1(b))不仅对于有意义可见内容的划界很重要，而且对于相对于相机头推断镜头探针的旋转也很重要。提出的用于跨连续帧跟踪边界轮廓和三角标记的方法与Fukuda等人[4]及Stehle等人[15]的工作有关。首先通过使用常规图像处理技术提取三角标记来推导倾斜视角内窥镜中的镜头旋转。该方法假设圆形区域的位置不在操作中改变，这一假设通常不成立，并且也不清楚它是否能够实时运行。Stehle等人提出了通过将二次曲线拟合至径向中检测到的边缘点来跨帧跟踪边界轮廓。其与本发明算法的主要差别在于本发明使用混合的串行+并行实现来渲染极性图像并且沿水平线执行搜索。这一策略实现了使用低计算成本的鲁棒和精度。

由于镜头探针相对于相机头移动，因此轮廓位置跨帧改变，而这可以使用初始离线估计来阻止。边界检测必须在每个帧时刻执行，其对所选算法的计算复杂性加以限制。若干问题阻止了分段圆形区域的单纯方法的使用：光通常传播至帧区域；取决于成像场景和光照条件，圆形图像会具有暗区；以及通常存在会影响分段性能的高亮、镜面和饱和。

合理的是假设待确定的曲线通常是具有5个自由度（DOF）的椭圆Ω[8]。于是，本发明提出在实现鲁棒性和准决定论运行时（quasi-deterministic runtime）之前使用这一形状跟踪跨帧的边界。设Ω_i-1如图4的（a）所示是帧时刻i-1处的曲线估计。这当前帧i的边界轮廓被如下更新：

1）在部分或完全包含边界的最后估计结果的原始图像（图4的（b））中定义环形区。

2）渲染其中在先边界估计被映射至中央垂线（图4的（c））的新条带图像。周知总是存在将任意椭圆映射至其中心是原点的单位圆的仿射变换[8]。这一变换S如下给出：

$S~\left(\begin{array}{ccc}r\cos \left(\mathrm{\φ}\right)& r\sin \left(\mathrm{\φ}\right)& -r(w_{i-1,x}\cos \left(\mathrm{\φ}\right)+w_{i-1,y}\sin \left(\mathrm{\φ}\right))\\ -\sin \left(\mathrm{\φ}\right)& \cos \left(\mathrm{\φ}\right)& w_{i-1,x}\sin \left(\mathrm{\φ}\right)-w_{i-1,y}\cos \left(\mathrm{\φ}\right)\\ 0& 0& 1\end{array}\right),$

其中r是Ωi-1的短轴和长轴之比，φ是长轴和水平方向之间的角度，而（wi-1，x，wi-1，y）则是圆锥中心wi-1的非同质坐标。变换S被用于生成图4的（b）的中间值，而极性图像（图4的（c））则通过施加从笛卡尔坐标到球形坐标的变化来获得。

3）使用图像处理技术来检测条带图像中的边界点。图4的（c）所示条带图像中的边界点通过应用1维边缘滤波器（例如，高斯的拉普拉斯算子）并从左到右扫描水平线以找出变换来检测。

4）由球形坐标χ=（ρ，θ）表达的边缘点通过式1的函数FS被映射回原始图像点。

$x~F_S\left(\mathrm{\χ}\right)~S^{-1}\left(\begin{array}{c}\mathrm{\ρ}\cos \left(\mathrm{\θ}\right)\\ \mathrm{\ρ}\sin \left(\mathrm{\θ}\right)\\ 1\end{array}\right),---\left(1\right)$

当前二次边界Ωi是最终使用鲁棒二次拟合[23]并使用RANSAC[24]来估计的，由此避免可能异常值的不利影响。

5）重复执行从1）至4）的步骤直到边界轮廓估计收敛为止。

被呈现用以估计圆形图像的边界轮廓的步骤使用计算单元的并行处理能力来加速。

E.在镜头经历选择运动时评估镜头探针相对于相机的角位移

在D部分描述的方法的正确收敛之后，边界轮廓被映射到极性图像（图4的（c））的中央垂线内，由此使得能够通过扫描略偏右的辅助垂线并选择具有最大强度的像素位置来进行三角标记的鲁棒检测，。三角标记的位置被馈入随机滤波器，后者相对于每帧处的校准位置推导当前角位移。

F.基于镜头探针的校准和角位移更新当前帧的投影模型

镜头和相机头之间的相对运动引起校准参数的改变，由此阻止校正基本的恒定模型的使用[5]。存在少数提出了这一问题解决方案的工作[2]–[4]，[15]，[22]，但是它们中的大多数都具有需要附加仪器来检测镜头旋转的缺点[2]，[3]，[22]。仅依赖于图像信息来推定相对运动的少数方法要么缺乏鲁棒性[4]，要么无法更新全组的相机参数[15]。该部分提出了由试验以及由光学器件布置的观念性理解同时驱动的用于倾斜视角内窥镜的新本征相机模型。假设探针围绕与图像平面正交但无需与光轴重合的轴旋转。更新相机模型所需的镜头旋转参数由接收关于边界轮廓和三角标记的图像信息作为输入的鲁棒EKF估计。本发明的动态校准方案相比于[15]具有两个重要优点：（i）整个投影模型作为镜头旋转的函数（而非仅仅是RD轮廓线曲线）被更新；以及（ii）镜头旋转能够在没有三角标记的情况下估计（参见图1（b））。令人信服的试验得到的在场景中重新投影的3D点确证了本方法。

为了评定相机模型中相机头和镜头之间的相对旋转的效果，本发明在完全旋转镜头探针的同时获取10幅校准图像。使用C部分中的方法为每个角位置估计相机校准，并且边界Ω和三角标记两者都如D和E部分所述。图8的（a）描绘了用于主点c、边界中心w以及镜头标记p的结果。由于这三个参数描述了几近完美的同心轨道，因此通过围绕与图像平面正交的轴旋转来建模镜头旋转对相机本征值看上去是合理的。这一想法已由Wu等人[3]推进，但是他们考虑的是轴总是通过主点，而这如我们所做的试验所示是通常无法成立的假设。

图5的（b）和（c）的方案旨在给出有关描述镜头旋转对相机本征值的影响的提议模型的想法。让我们假设管道镜镜头将虚拟图像投影至位于远端的平面I′。我们可以将I′看作是通过相机目镜可直接看见的图像。K_c是该虚拟投影的本征矩阵，而c′则是光轴与平面相交的点位置。现在假设相机头连接至目镜，以使得CCD平面I完美平行于I′并与光轴相交。对I的投影具有本征值K_h，其具有作为c′的图像的主点c。由此，如果K_h是本征矩阵，则估计

$K_0~\left(\begin{array}{ccc}f& 0& c_x\\ 0& f& c_y\\ 0& 0& 1\end{array}\right)$

于是其可被分解为

$K_0~K_h{K}_c~\left(\begin{array}{ccc}{f}_h& 0& c_x\\ 0& f_h& c_y\\ 0& 0& 1\end{array}\right)\left(\begin{array}{ccc}{f}_c& 0& 0\\ 0& f_c& 0\\ 0& 0& 1\end{array}\right)$

其中f_c是管道镜镜头的焦距，而f_h是将度量单元转换成像素的相机头的焦距。

现在考虑镜头探针以角度α围绕轴l旋转(图8的(c))。l被假设为与虚拟平面I′正交，但无需与镜头轴重合。在此情况下，点c′描述具有幅度α的圆弧，并且由于I和I′平行，因此其图像也平行c。由相机头和旋转管道镜形成的复合光学***的本征矩阵变为

K～K_hR_α，q″K_c  (8)

其中R_α，q′是围绕点q′旋转α的平面旋转，其中轴l与I′相交。

$R_{\mathrm{\α},q^{\′}}=\left(\begin{array}{ccc}\cos \left(\mathrm{\α}\right)& \sin \left(\mathrm{\α}\right)& (1-\cos \left(\mathrm{\α}\right))q_x^{\′}-\sin \left(\mathrm{\α}\right)q_y^{\′}\\ -\sin \left(\mathrm{\α}\right)& \cos \left(\mathrm{\α}\right)& \sin \left(\mathrm{\α}\right)q_x+(1-\cos \left(\mathrm{\α}\right))q_y\\ 0& 0& 1\end{array}\right),---\left(9\right)$

q′的位置显然未通过旋转改变，并且对其图像q～K_hq′也是一样。考虑K_h的具体结构，可以按如下方式重写式8

K～R_α，qK_hK_c

～R_α，qK₀

我们刚刚已经导出用于适应镜头探针旋转并与图8的(a)的查看一致的内窥镜相机的投影模型。初始化程序估计任意基准位置(α＝0)处的相机校准K₀，ξ。在特定帧时刻i，本征参数矩阵变为

$K_i~R_{{\mathrm{\α}}_i,q_i}{K}_0---\left(10\right)$

其中α_i是镜头的相对角位移，而qi是在旋转期间保持固定的图像点。由于径向失真是镜头的特征，因此参数ξ不受相对于相机头的相对运动影响。于是，从式2，随后是在相机坐标帧内表示的一般3D点x成像到：

x～K_iΓ_ξ(I₃0)X  (11)

本征参数矩阵在每个帧时刻的更新需要知晓相对角位移α_i和图像旋转中心q_i。我们现在描述如何从边界轮廓Ω和三角标记p的位置推导出这些参数(部分III)。

设w_i和w₀分别是当前和基准帧内边界轮廓Ω_i和Ω₀的中心。类似地，p_i和p₀是两图像中三角标记的位置。我们假设w_i，w₀和p_i，p₀都由其我们旨在估计其参数的平面旋转

相关。这一情形如图6的(a)所示，图中可以轻易看出旋转中心q_i必须是由w_i，w₀和p_i，p₀定义的线段的二等分线的交叉。一旦q_i已知，旋转角α_i的估计就微不足道。只要三角标记未知（如果其不存在或无法被检测），q_i的估计就需要最小三个不同的边界轮廓（图6的（b））。

为了避免限定不足的情况并增加对w和p测量中误差的鲁棒性，使用随机EKF[12]估计旋转参数。状态转换假设用于运动和固定旋转中心的恒定速度模型。该式对于状态变量是线性的

$\left(\begin{array}{c}{\mathrm{\α}}_{i+1}\\ {\stackrel{\·}{\mathrm{\α}}}_{i+1}\\ q_{i+1}\end{array}\right)=\left(\begin{array}{cc}T& 0_{2\×3}\\ 0_{3\×2}& I_3\end{array}\right)\left(\begin{array}{c}{\mathrm{\α}}_i\\ {\stackrel{\·}{\mathrm{\α}}}_i\\ q_i\end{array}\right)$

其中T依赖于帧获取间隔δt

$T=\left(\begin{array}{cc}1& \mathrm{\δt}\\ 0& 1\end{array}\right).$

测量式对于α_i和q_i是非线性的

$\left(\begin{array}{c}{w}_i\\ p_i\end{array}\right)=\left(\begin{array}{cc}{R}_{{\mathrm{\α}}_i,q_i}& 0_{3\×3}\\ 0_{3\×3}& R_{{\mathrm{\α}}_i,q_i}\end{array}\right)\left(\begin{array}{c}{w}_0\\ p_0\end{array}\right),$

其中只要三角标记检测失败就丢弃最后两个式子。

通过将网格边角重新投影至为不同角度α获取的棋盘图案的图像上来证实该提出的模型。针对基准位置执行SIC校准[21]，以使得能够确定矩阵K0、RD参数ξ和网格点的3D参数。随后，在不移动镜头探针的情况下小心旋转相机头以保持与校准图案的相对姿态。使用图6中的几何构造为每个帧估计图像旋转中心q_i和角位移α_i。最终，使用式11将3D网格点投影到帧上，并且测量到实际图像边角位置的误差距离。

图7描绘了针对不同角位移α的重新投影误差的RMS。其值在2到5个像素之间变化，并且没有查看到***性行为。本发明仅关注本征参数，而[2]和[3]考虑的则是本征和外源校准两者，并且利用附加仪器。虽然无法做出直接比较，但是需要说明的是在仅考虑接近图像中心的点的情况下本发明的重新投影误差小于[3]且等于[2]。根据上文且针对所有推测，试验结果明确证实了提出的模型。

G.通过图像扭曲实时校正径向失真

该部分讨论作为可视化***最终输出的校正透视图像的渲染。如[7]中指出的，应该使用逆映射法执行具体变换下的图像的有效扭曲。于是，我们必须导出将期望的无失真图像中的点y映射至原始失真帧中的点x的函数F。从式11可得

$F\left(y\right)~K_i{\mathrm{\Γ}}_{\mathrm{\ξ}}\left(R_{-{\mathrm{\α}}_i,q_i^{\′\′}}{K}_y^{-1}y\right).$

K_y规定无失真图像的特定特征（例如，中心、分辨率），R_-αi，q″将扭曲结果旋转回原始定向，而q″则是旋转中心q_i的背投影。

$q_i^{\′\′}~{\left(\begin{array}{ccc}{q}_{i,x}^{\′\′}& q_{i,y}^{\′\′}& 1\end{array}\right)}^T~{\mathrm{\Γ}}_{\mathrm{\ξ}}^{-1}\left(K_i^{-1}{q}_i\right).$

为了在中心区域中维持对象的标度，我们扩展图像的外周并保持无失真中心区域的尺寸。这通过从D部分的边界轮廓的半径计算扭曲图像的尺寸u来实现。设rd是原点和点K_i-1p₀（扭曲半径）之间的距离。期望的图像尺寸u由

给出，其中f是相机焦距，而r_u是无失真的半径。因此，矩阵K_y必须是

$K_y~\left(\begin{array}{ccc}f& 0& -f{q}_{i,x}^{\′\′}\\ 0& f& -f{q}_{i,y}^{\′\′}\\ 0& 0& 1\end{array}\right)$

其中扭曲图像的中心是图像旋转中心q_i被映射的点。图8示出了用于内窥镜视频序列的某些帧的RD校正结果。例子清楚地显示了场景感知的改善以及在图像几何结构校正期间计及镜头探针旋转的重要性（图8的（b））。

经校正图像的渲染需要高性能计算资源以实时处理图像。本发明提出的RD校正算法则适于在并行处理单元（在此情况下为GPU）上并行计算。相比于如下4种方法，本发明提升了自身的串行+并行解决方案的影响：

（i）单纯基于CPI的解决方案；

（ii）使用OpenMP2directives[6]的假想CPU版本，其能被用于共享存储器架构，诸如x86系列的常规COTS多核CPU。

（iii）未经优化的混合CPU+GPU版本；以及

（iv）经优化的混合CPU+GPU版本。

完全运行***的试验设置包括装备有2.40GHz的

CoreTM2Quad CPU、4Gbyte的RAM以及GeForce8800GTXGPU的工作站。GPU具有128个流处理器（SP），各自以1.35GHz运行并且总视频RAM为768Mbyte。源代码以GCC版本4.3.4编写，并且使用CUDA版本3.0。

图9比较了使用上述解决方案实现的执行时间（执行单帧所需的时间）。作为示出GPU要比CPU快得多的度量的加速表示出了单纯CPU和经优化的混合CPU+GPU执行之间之比。图9给出的比较示出了CPU无法像GPU那样高效地处理大型扭曲操作（算法中最为耗时的任务）。该表示出了对于输出尺寸在3000x3000像素之上的图像，不仅GPU加速执行，CPU性能会由于高速缓存未命中损失而恶化。

图10示出了径向失真校正算法的混合串行+并行实现步骤以及在CPU和GPU上运行的混合程序的开发中适用的实现策略的细节。GPU上算法的并行化被分成三个主要步骤：

(i)图像转换：将图像分割成其RGB通道并执行灰度级转换。

(ii)边界估计：灰度级图像被限至并行处理单元存储器并且沿着轮廓映射至球形坐标（式1），启动增强内核。

(iii)径向失真校正：来自图像转换阶段的三个图像通道被结合到并行处理单元存储器，并且启动RD校正内核。

这是CPU中计算的中间步骤（在图10中隔离）。针对椭圆拟合至的轮廓点扫描极性图像的每条线。椭圆参数被馈送至在F部分讨论的EKF，并且校准参数被更新并作为自变量传递至失真校正内核。该步骤在CPU而非GPU上实现是由于其涉及处理的串行特性。

优化的混合CPU+GPU处理方法依赖于数据预对准程序步骤，后者允许执行结合[1]已知的每线程组单存储器访问。如果相反地要由CPU执行的数据未对准，则不执行合并访问，并且针对每组线程进行若干存储器交易，由此大幅降低了内核下执行的性能。

虽然没有使用图像的阿尔法通道，但是需要履行用于完全执行结合访问的存储器布置请求。传输数据量的增加在传输时间中引入了10.6%的害处，而实现的结合则将内核执行时间缩减了66.1%。总之，该结合实现相比于未经优化的混合CPU+GPU实现节省了38.7%的计算时间。

图12的表格则示出了在处理150帧视频序列时图10中每个步骤所需时间（时间表示用于处理每种分辨率的所有帧的平均时值）。图像转换（I.C.）时间还包括来自主机的输入图像到设备的全局存储器的传输时间。边界估计（B.E.）时间则考虑边界点的提取，椭圆拟合以及在CPU上执行的EKF估计处理（图14中的紫框）。注意到都在GPU上计算的输入图像到球形空间的映射以及轮廓增强滤波器消耗表II中表示的边界估计计算时间的1至3毫秒。径向失真校正（R.D.C.）还考虑从设备的全局存储器到主机的传输时间。

参考文献

[1]D.Kirk and W.-m.Hwu,Programming Massively ParallelProcessors:A Hands-on Approach.Morgan Kaufmann,2010.

[2]T.Yamaguchi,M.Nakamoto,Y.Sato,Y.Nakajima,K.Konishi,M.Hashizume,T.Nishii,N.Sugano,H.Yoshikawa,K.Yonenobu,and S.Tamura,“Camera Model and CalibrationProcedure for Oblique-Viewing Endoscope,”in MICCAI,2003,pp.373–381.

[3]C.Wu,B.Jaramaz,and S.G.Narasimhan,“A FullGeometric and Photometric Calibration Method for Oblique-viewingEndoscope,”International Journal of Computer Aided Surgery,vol.15,pp.19–31,2010.

[4]N.Fukuda,Y.Chen,M.Nakamoto,and T,“A scope cylinderrotation tracking method for oblique-viewing endoscopes withoutattached sensing device,”Software Engineering and Data Mining,no.1,pp.684–687,2010.

[5]J.Barreto,J.Santos,P.Menezes,and F.Fonseca,“Ray-based Calibration of Rigid Medical Endoscopes,”in OMNIVIS,2008.

[6]B.Chapman,G.Jost,and R.Van Der Pass,Using OpenMP:Portable Shared Memory Parallel Programming(ScientificComputation and Engineering Series).The MIT Press,2008.

[7]L.Velho,A.C.Frery,and J.Gomes,Image Processing forComputer Graphics and Vision.Springer Publishing Company,Incorporated,2008.

[8]R.I.Hartley and A.Zisserman,Multiple View Geometry inComputer Vision,2nd ed.Cambridge University Press,ISBN:0521540518,2004.

[9]W.Smith,N.Vakil,and S.Maislin,“Correction of Distortionin Endoscope Images,”IEEE Transactions on Medical Imaging,vol.11,no.1,pp.117–122,Jan1992.

[10]K.VijayanAsari,S.Kumar,and D.Radhakrishnan,“A newapproach for nonlinear distortion correction in endoscopic imagesbased on least squares estimation,”IEEE Transactions on MedicalImaging,vol.18,no.4,pp.345–354,1999.

[11]J.Helferty,C.Zhang,G.McLennan,and W.Higgins,“Videoendoscopic distortion correction and its application to virtualguidance of endoscopy,”IEEE Transactions on Medical Imaging,vol.20,no.7,pp.605–617,2001.

[12]S.M.Bozic,Digital and Kalman Filtering.London:EdwardArnold,1979.

[13]J.P.Barreto,“A unifying geometric representation forcentral projection systems,”Comput.Vis.Image Underst.,vol.103,no.3,pp.208–217,2006.

[14]A.Fitzgibbon,“Simultaneous linear estimation of multipleview geometry and lens distortion,”in CVPR,vol.1,2001,pp.I–125–I–132vol.1.

[15]T.Stehle,M.Hennes,S.Gross,A.Behrens,J.Wulff,and T.Aach,“Dynamic Distortion Correction for Endoscopy Systems withExchangeable Optics,”in Bildverarbeitungfür die Medizin2009.Berlin:Springer,2009,142–146.

[16]J.-Y.Bouguet.Camera Calibration Toolbox for Matlab.[Online].Available:http://www.vision.caltech.edu/bouguetj/calibdoc/index.html#ref

[17]Z.Zhang,“Flexible camera calibration by viewing a planefrom unknown orientations,”in ICCV,1999,pp.666–673.

[18]R.Shahidi,M.Bax,C.Maurer,J.Johnson,E.Wilkinson,B.Wang,J.West,M.Citardi,K.Manwaring,and R.Khadem,“Implementation,calibration and accuracy testing of animage-enhanced endoscopy system,”IEEE Transactions on MedicalImaging,vol.21,no.12,pp.1524–1535,Jan2002.

[19]C.Wengert,M.Reeff,P.Cattin,and G.Sz′ekely,“Fullyautomatic endoscope calibration for intraoperative use,”inBildverarbeitungfür die Medizin2006,2006,pp.419–423.

[20]J.Mallon and P.F.Whelan,“Which pattern？Biasingaspects of planar calibration patterns and detection methods,”Pattern Recognition Letters,vol.28,no.8,pp.921–930,Jan2007.

[21]J.Barreto,J.Roquette,P.Sturm,and F.Fonseca,“Automatic Camera Calibration Applied to Medical Endoscopy,”inBMVC,2009.

[22]S.D.Buck,F.Maes,A.D’Hoore,and P.Suetens,“Evaluation of a novel calibration technique for optically trackedoblique laparoscopes,”Proceedings of the10th internationalconference on Medical image computing and computer-assistedintervention-Volume Part I,pp.467–474,2007.

[23]A.Fitzgibbon,M.Pilu,and R.Fisher,“Direct least squarefitting of ellipses,”IEEE Transactions on Pattern Analysis andMachine Intelligence,vol.21,no.5,pp.476–480,May1999.

[24]M.A.Fischler and R.C.Bolles,“Random SampleConsensus:A Paradigm for Model Fitting with Applications to ImageAnalysis and Automated Cartography,”Commun.ACM,vol.24,no.6,pp.381–395,1981.

Claims (13)

1.一种设备，包括如下部分：

(i)由不透明材料制成以避免光从外部进入的盒(1)，其具有允许镜头探针进入的一个或多个开口(2)，每个开口由用于将从外部进入的光最小化的橡胶膜(3)围绕；该盒由透明或半透明材料的隔离物(6)分隔成两个隔室，一个开口隔室能够由用于更好地重复将执行视觉检查的部件的流体填充；

(ii)由透明或半透明材料制成且被置于隔离物顶部上的具有已知图案的平面网格(7)；

(iii)安装在下盒隔室(5)中且用于平面网格的漫射背光照明的光源；

(iv)具有顺序和并行处理能力且与用于图像获取的相机相连的计算机单元，诸如中央处理单元和图形处理单元；

(v)连接至所述计算机单元的显示设备；以及

(vi)连接至所述计算机单元的用于输入用户命令的接口设备。

2.如权利要求1所述的设备，包括平面棋盘图案，该平面棋盘图案通过所述相机查看并且能够作为结果包括位于两个连续对角线白色方形的中心处并标识位于这两个连续对角线白色方形之间的网格边角的可区分特征。

3.如权利要求1所述的设备，还能够包括用于查看所述平面图案的多个可能进入点，这些进入点根据镜头探针的类型具有多种形状和尺寸并且用于从不同视角查看所述图案。

4.一种用于自动相机校准、根据可能的镜头探针旋转更新投影模型以及实时图像失真校正的方法，所述方法执行如下操作：

(i)使用校准盒棋盘图案的一幅或多幅图像的用户命令获取；

(ii)检测校准图像中的网格边角并建立图像平面对应；

(iii)使用至少一幅校准图像中的图像平面对应确定相机本征参数，诸如焦距、主点、偏斜和纵横比以及镜头失真；

(iv)在诸如内窥镜相机的医疗或工业设备中连续获取图像/视频；

(v)在适用时检测圆形和帧区域之间的边界轮廓、三角标记的图像坐标，该操作在每个帧时刻执行；

(vi)在适用时使用(v)的输出鲁棒确定镜头探针经受的相对于相机头的角位移；

(vii)考虑到(iii)的初始相机校准来更新本征校准信息，并且在适用时使用用于描述倾斜视角内窥镜或其镜头相对于相机头旋转的任何其他相机设备的新投影模型更新(vi)的相对旋转信息；

(viii)使用(vii)的当前相机校准执行用于校正径向失真的实时图像修正。

5.如权利要求4所述的方法，其中操作(ii)包括在具有棋盘图案的平面的边角点和校准图像的边角点之间自动建立对应，并且还包括如下步骤：

a)使用图像处理技术检测中央图像区域中的边角或种子点；

b)推断该边角或种子点在网格平面内的坐标；

c)使用相机参数、镜头失真以及网格平面和相机之间的相对姿态的隐式编码确定将平面网格坐标中的点映射至图像坐标中的点的函数；

d)使用在(c)中确定的函数将平面网格边角投影到图像内；

e)通过围绕投影的平面网格边角的局部搜索定位图像中的新边角；

f)迭代执行步骤(c)至(e)，直到认为建立的边角对应的数目令人满意为止。

6.如权利要求4所述的方法，其中操作(iii)包括使用权利要求4中建立的边角对应从单幅校准图像估计相机本征参数和镜头失真，并且还包括如下步骤：

a)计算将平面网格内的提升点坐标映射至图像平面中的提升点坐标的函数；

b)分解(a)的映射函数以确定焦距、主点、纵横比、偏斜、失真参数、以及网格平面和相机之间的相对姿态、旋转和平移。

7.如权利要求6所述的方法，其能够由用于实现如下一个或多个目的的迭代非线性优化的附加步骤所补充：

a)细化估计的参数，

b)融合由不同校准图像提供的估计结果，

c)使用与为如权利要求5和6所述的计算方法假设的模型函数不同的模型函数估计径向失真。

8.如权利要求4所述的方法，其中操作(v)包括无监督地实时跟踪圆形图像(CI)的边界轮廓并在三角标记可见时检测三角标记，并且包括如下步骤：

a)在部分或完全包含CI边界的最后估计结果的原始图像中定义环形区；

b)渲染其中在先边界估计结果被映射至中央垂线的新条带图像；

c)使用图像处理技术检测条带图像中的边界点；

d)将检测到的点映射回原始图像并使用鲁棒的估计结果拟合二次曲线；

e)重复步骤(a)至(d)直到边界轮廓估计结果收敛为止；

f)通过扫描略偏右的条带图像中的辅助垂线，选择具有最大强度的像素位置并将该像素位置映射回原始图像，来鲁棒检测三角标记。

9.如权利要求4所述的方法，其中操作(vi)包括在镜头探针经历旋转运动时估计相对于镜头位置的角位移，用以获取用作权利要求5,6和7的方法中的输入的校准帧，并利用在适用时由权利要求8针对基准位置获取的边界轮廓和三角标记以及利用视频的连续获取，并且对于当前帧还包括如下步骤：

a)如权利要求8所述在适用时检测边界轮廓和三角标记；

b)将(a)的信息馈入随机滤波器，该随机滤波器推断当前角位移，由此推断使用专门图像信息确定的相对旋转。

10.如权利要求4所述的方法，其中操作(vii)包括在镜头探针经历旋转运动时，在考虑到在如权利要求5,6和7所述获取的初始基准位置处的相机本征校准以及如权利要求9所述获取的镜头探针相对于该基准位置的当前角位移的情况下，推断用于当前帧的本征校准和投影模型，该本征校准和投影模型包括焦距、主点、纵横比、偏斜和失真参数。

11.如权利要求4所述的方法，其中操作(viii)包括基于并行处理并利用常规计算单元***通过图像扭曲实时校正径向失真，所述方法包括：

a)将数据从主机计算单元传输到并行处理单元；

b)在***并行处理单元上处理图像；

c)将最终处理的图像传输回主机计算单元；

d)将结果可视化。

12.一种在计算机介质上具体化的计算机程序，在由处理器执行时实现权利要求4至11中任一项的方法步骤。

13.使用并行处理单元的并行计算能力加速如权利要求8至11所述的方法。

Images