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

摘要

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

说明书

技术领域

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

背景技术

通常需要光学仪器产生给定对象的几何组成图像，其中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\times 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(\phi \right)& r\sin \left(\phi \right)& -r(w_{i-1,x}\cos \left(\phi \right)+w_{i-1,y}\sin \left(\phi \right))\\ -\sin \left(\phi \right)& \cos \left(\phi \right)& w_{i-1,x}\sin \left(\phi \right)-w_{i-1,y}\cos \left(\phi \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(\chi \right)~S^{-1}\left(\begin{array}{c}\rho \cos \left(\theta \right)\\ \rho \sin \left(\theta \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_{\alpha ,q^{\prime }}=\left(\begin{array}{ccc}\cos \left(\alpha \right)& \sin \left(\alpha \right)& (1-\cos \left(\alpha \right))q_x^{\prime }-\sin \left(\alpha \right)q_y^{\prime }\\ -\sin \left(\alpha \right)& \cos \left(\alpha \right)& \sin \left(\alpha \right)q_x+(1-\cos \left(\alpha \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_{{\alpha }_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}{\alpha }_{i+1}\\ {\stackrel{·}{\alpha }}_{i+1}\\ q_{i+1}\end{array}\right)=\left(\begin{array}{cc}T& 0_{2\times 3}\\ 0_{3\times 2}& I_3\end{array}\right)\left(\begin{array}{c}{\alpha }_i\\ {\stackrel{·}{\alpha }}_i\\ q_i\end{array}\right)$

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

$T=\left(\begin{array}{cc}1& \mathrm{\delta 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}_{{\alpha }_i,q_i}& 0_{3\times 3}\\ 0_{3\times 3}& R_{{\alpha }_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{\Gamma }_{\xi }\left(R_{-{\alpha }_i,q_i^{\prime \prime }}{K}_y^{-1}y\right).$

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

$q_i^{\prime \prime }~{\left(\begin{array}{ccc}{q}_{i,x}^{\prime \prime }& q_{i,y}^{\prime \prime }& 1\end{array}\right)}^T~{\Gamma }_{\xi }^{-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}^{\prime \prime }\\ 0& f& -f{q}_{i,y}^{\prime \prime }\\ 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以及GeForce8800GTX GPU的工作站。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 Parallel  Processors: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 Calibration  Procedure for Oblique-Viewing Endoscope,”in MICCAI,2003,pp. 373–381.

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

[4]N.Fukuda,Y.Chen,M.Nakamoto,and T,“A scope cylinder  rotation tracking method for oblique-viewing endoscopes without  attached 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(Scientific  Computation and Engineering Series).The MIT Press,2008.

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

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

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

[10]K.VijayanAsari,S.Kumar,and D.Radhakrishnan,“A new  approach for nonlinear distortion correction in endoscopic images  based on least squares estimation,”IEEE Transactions on Medical  Imaging,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 virtual  guidance of endoscopy,”IEEE Transactions on Medical Imaging,vol. 20,no.7,pp.605–617,2001.

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

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

[14]A.Fitzgibbon,“Simultaneous linear estimation of multiple  view 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 with  Exchangeable 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/calib  doc/index.html#ref

[17]Z.Zhang,“Flexible camera calibration by viewing a plane  from 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 an  image-enhanced endoscopy system,”IEEE Transactions on Medical  Imaging,vol.21,no.12,pp.1524–1535,Jan2002.

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

[20]J.Mallon and P.F.Whelan,“Which pattern？Biasing  aspects 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,”in  BMVC,2009.

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

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

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