立体影像调度的实时核线重采样方法研究

桂坡坡，叶勤，邓非

(1.同济大学 测绘与地理信息学院，上海 200092；2.武汉大学 测绘学院，武汉 430079)

1 引 言

随着摄影测量学、计算机技术的发展，数字摄影测量系统无论是理论上还是实践应用中都取得了重大突破与应用。在数字摄影测量双像测图获取立体信息中，核线影像占有着重要地位。核线影像是原始立体影像对沿核线进行重采样而形成，在核线影像上同名像点被约束在同一行，因此是计算机立体显示的基础。核线影像大大增强计算机的立体显示效果，提高立体量测与匹配的精度及速度，是数字摄影测量数据处理中自动获取数字地面高程的基础，因此在数字摄影测量中具有重要意义。

现有以共线方程为基础的摄影测量理论已相对完善，国内、外主流的商用数字摄影测量系统，如JX4、Virtuozo、Inpho等[1]。在进行立体测图获取信息时，核线影像的生成流程都采用在相对定向完成后、测图前对整个像对进行核线重采样的方式；即在定向模块中，完成相对定向后即对立体像对进行核线重采样，生成核线影像。测图时，测图模块直接调用核线影像进行显示、计算分析。这种策略需要留出专门的核线影像生成时间、以及专门的存储空间；若要快速实时地实现立体显示以及立体信息获取，则这种处理策略无法满足要求。

但是随着摄影测量的快速发展，目前用于摄影测量或辅助摄影测量的数据源越来越多，数字摄影测量工作站一方面面临着多源数据联合进行摄影测量的问题，另一方面面临实时快速的生成产品缩短内业成图周期的挑战。而影像对数据网上立体发布以及某些场合如实时街景立体测量中实时化的要求，这些都迫切需要在保证量测精度的同时提高摄影测量系统的工作效率、实时性(这包括定向完成后立体显示的实时性)，为此要求实时快速地显示与生成产品，缩短内业处理的时间，此时对涉及的立体影像全部生成核线影像从时间上已不允许。

因此，本文针对当前数字摄影测量系统常用工作流程存在的实时性差、耗时长并且要有大量空间来存储生成的核线影像的缺点，运用实时核线重采样的思想，在定向后对定向模块下的影像对不再生成核线影像，而是在测图过程中，采用立体显示加速策略，在内存中实时快速生成所需要影像区域的核线影像，同时将影像分层分块，采用线程调度的方法提高影像的刷新效率，在实施立体测图过程中实现核线影像的实时重采样和满足视觉要求的调用。

2 核线重采样矩阵

核线重采样算法中，确定核线重采样矩阵既是生成核线影像的关键也是实时核线重采样的基础。文献[2-4]中给出了几种核线重采样算法，但是这些方法的数据针对性较强，普适性不高。文献[1]及文献[4]中提到的基于相对方位元素的核线生成算法在数字摄影测量系统中有很强的应用性。计算机视觉中，常采用匹配的策略通过计算本质矩阵进行对极几何的恢复，实时性较强，如文献[5-7]，但是限于航空影像的图幅及拍摄高度，匹配的效率与正确率都将下降，并且恢复的核线影像对地形图测制精度较低。本文从普适性及精度的角度出发采用了文献[8]中的两步旋转法。

根据摄影测量基本原理[8]，在已知立体影像对的两张像片间的相对定向元素后，两张像片相对于所选的摄影测量参考坐标系的姿态角(φ,ω,κ)已知。从原始影像得到核线影像只需要两步旋转变换，在重采样分辨率保持不变的情况下，必然存在重采样矩阵Rn使得原始影像下像点坐标(x，y)到核线影像下像点坐标(x′，y′)存在唯一映射，即

X′=RnX

(1)

图1 原始影像转换核线影像示意图

第一步：将原始影像旋转到绝对垂直位置。这里绝对垂直位置是指影像平面与所选的相对定向摄影测量参考坐标系下XOY平面平行的位置。此步旋转的旋转矩阵通过像片相对于摄影测量参考坐标系的3个角元素来确定，即可完成原始影像旋转到绝对垂直位置的转换。

X⊥=RX

(2)

其中，

(3)

φ,ω,κ为影像拍摄时刻的相机相对于摄影测量参考坐标系的姿态角，c，s分别为三角函数cos，sin的缩写。

第二步：将影像由绝对垂直位置旋转到核线影像位置，此步涉及绕基线旋转影像，此时旋转矩阵用Rb表示。

X′=RbX⊥

(4)

(5)

首先绕Y轴旋转绝对垂直影像Φ，其次绕Z轴旋转K。为使核线影像的大小最小，采用式(6)所得Ω角度继续绕X轴旋转。式(6)中ϖ′和ϖ″分别代表立体影像左右片与相对定向所选参考坐标系X轴的旋转角(如果参考坐标系是地面已知坐标系，则为左右片外方位角元素中绕X轴的旋转角)。

(6)

Rb=RΚ*RF*RΩ

(7)

因此有：

(8)

令

(9)

(10)

(11)

式(10)给出了由原始影像点坐标计算核线影像点坐标的正算式，式(11)给出的是由核线影像点坐标计算原始影像上点坐标的反算式。

由式(9)知核线采样矩阵只与影像对左右片之间的相互位置关系元素有关。对于像对中左右影像，绕基线的旋转矩阵Rb是共有的，Rb与各自像片旋转矩阵相乘即可以得到对应影像的核线重采样矩阵。根据式(10)、式(11)，进行核线影像坐标与原始影像坐标转换时，需要知道原始影像与核线影像的主距，实际解算中，一般假设两者主距相等，同为原始影像的主距f。

3 基于立体显示加速策略的核线影像实时性采样

在实时生成核线影像测图前，需对影像进行统一预处理，生成金字塔影像。影像定向完成后，确定影像对4个角点在核线影像上的相应位置。在实时测图前，不对已定向的影像对进行核线重采样生成核线影像，而是直接在测图过程的显示环节，在内存中完成视场区域内影像的核线重采样。由于只对视场区域内的影像进行核线生成，所以涉及的计算数据比整个影像的核线生成明显减少，计算速度快；但由于视场会随着测图区域的变化而移动，所以在实时性采样过程中要研究立体显示加速策略。

整个基于实时性采样的立体测图显示过程为：影像静态预处理、实时预处理、基于立体显示加速策略的影像块实时采样。最终实现立体显示、为测图做准备。

3.1 影像预处理

预处理是影像块实时采样前的数据准备工作，以提高实时立体显示时影像的加载速度。预处理包括静态预处理与实时预处理。

3.1.1 静态预处理

静态预处理是对原始影像建立影像金字塔，以原始影像为单位进行，与像对信息无关。影像金字塔是以影像宽与高减半的方式向上创建，每一层金字塔影像的宽和高都是其下一层金字塔影像的一半且是其上一层金字塔的两倍。在实时显示的调用过程中，通过视口显示的影像缩放比确定需要调用影像的金字塔层级，如当缩放倍数为8时，根据缩放倍数调用金字塔第3层影像。

影像金字塔的建立除可以提高影像的调度速度，在大范围粗略显示时还可以减少了内存的冗余处理。当显示范围很大、影像缩小不需要精细的显示影像时，可不必对原始影像对进行核线重采样，而是直接调用金字塔上层影像。这减少了内存中处理的像素数据、缩短了处理时间，保证了缩放的流畅性。

3.1.2 实时预处理

实时预处理是为实时显示影像对，进行立体观察做准备。需对影像对的加载或切换进行操作，为了保证视觉的连续性与稳定性，这些处理有很强的时间要求被称为实时预处理。它主要包括计算重采样矩阵及核线影像分块等。核线重采样矩阵只与影像对中左右片的相互位置关系元素有关，当影像相对定向完成后，左右片的核线重采样矩阵唯一确定。实时显示影像对不变时，核线重采样矩阵保持不变；因此实时预处理过程是以像对为单位进行的，在显示影像对改变时，重采样矩阵也相应发生变化，所以该预处理被称为实时预处理。

在相对定向完成、已知立体影像对的两张像片间的相对定向元素后，核线重采样矩阵也已知，此时根据正算式(10)计算左右像片影像4个角点(0，0)(0，H)(W，H)(W，0)对应在核线影像上的坐标(X1，Y1)(X2，Y2)(X3，Y3)(X4，Y4)，确定角点中X、Y的最大最小值，分别以X、Y最大最小值之差确定影像的宽与高，以X、Y最小值为影像的左下点。在影像对的立体显示中，通用的三维渲染引擎(如OpenGL等)对影像大小为2的整数次方的支持较好，同时综合考虑立体刷新的效率与速度，研究中采用256*256的分块大小，将确定的核线影像范围分成若干影像块，右侧及底部不足整数块时以实际大小处理。

经过影像预处理，任一像点坐标对应在任意金字塔层上的坐标、以及在核线影像上对应的分块号和块内坐标关系都被确定。因此已知当前块上某点坐标，通过块号以及金字塔层级即可以算出在核线影像上的坐标，根据反算公式(11)可计算出该点在原始影像上对应坐标，并进行前方交会等运算。

3.2 立体显示加速策略

为提高海量遥感影像的读取与显示速度，影像的分层与分块是常用策略，文献[9][10][11]都采用了这种处理方式。本文中针对航空摄影立体影像对进行分层分块，同时完成两张影像的调度，通过当前测标位置的Z值及影像对的相对位置关系，确定对应的影像块。

在显示算法实现过程中，为提高影像的显示刷新速率，将每一个影像块作为一个线程在后台进行调度。每一帧中都要判断是否有影像块移入或移出当前视口，对移入移出当前视口的影像块分别进行加载和卸载。为保障影像刷新的流畅性，在实现中添加缓存器，其实质是一个存储列表，存储着下一帧可能出现在当前视口内的所有影像块，即沿当前显示影像块四周的影像块。如图2所示，图中每小格代表一个影像块，其中白色为当前显示影像块，灰色区域为缓存器中保存的影像块。

影像块的加载与卸载过程是与缓存器进行交互，从当前缓存器中进行影像块的加载，卸载的影像块并不立即清除，而是重新保存在缓存器中。这种设置避免了程序直接与硬盘数据进行通信，可加快数据的处理，提高显示速度。同时，缓存器通过线程在后台根据当前位置信息不断计算更新可能在下一帧出现的影像块。

图2 当前视口影像块与缓存器中影像块关系图

3.3 影像块实时采样

影像块的采样作为一个线程在后台进行运行，因此核线影像生成的实时性一方面体现在核线影像存储在当前内存中，另一方面体现在核线影像的获取是根据当前视口而决定的。

影像块的采样算法流程如表1所示。

表1 影像块的采样算法

根据当前块号计算其在核线影像上的角点坐标(x1′，y1′)(x2′，y2′)(x3′，y3′)(x4′，y4′)反算至原图像上(x1，y1)(x2，y2)(x3，y3)(x4，y4)，确定4点中x、y的最大最小值，分别以x、y最大最小值之差确定影像的宽与高，以x、y最小值为影像的左上点，在对应影像金字塔层上取值。

采样生成显示影像块的核线影像时，遍历该核线影像块、对影像块上每点利用式(11)反算其在原图像上坐标，采用双线性插值法进行点的灰度插值，得到此时需显示影像块的核线影像。

4 实验分析

图3 立体像对

实验用计算机配置为：8GB内存、Intel Core i7(2.9GHz)CPU。实验数据如图3所示，是由胶片扫描得到的数字原始影像，主距f为210.681mm，影像大小皆为8328*8375，像素大小为0.02799mm/pixel。该立体像对的基线向量为(322.95，-9.32，1.15)，各分量根据摄影比例尺调整到与实地大小相接近，单位为m。左右影像的内标定系数及姿态角如表2所示，其中，a10、a11、a12、a20、a21、a22为内定向参数。

表2 像对左右片的各项参数

由于式(2)中的算式都是建立在以像主点为原点的像平面坐标系上的，因此若基于扫描影像进行核线重采样，则应该对获得的影像像素坐标(i，j)进行内定向，将像素坐标转换为像平面坐标系中的坐标。

内定向公式为：x=a10+a11i+a12jy=a20+a21i+ka22j，此时可将像点的像空间坐标表示为如下形式：

(12)

根据表1结合式(9)和式(12)计算左片的核线重采样矩阵为：

表3是由立体像对重叠区域选取6个地面模型点(由于左片摄站点为模型的原点，其在地面坐标系中的坐标X、Y、Z大约为：500934、543180、910，所以地面模型点的X、Y会出现负值，而地面点的Z都为负值、在-900m左右)，来验证算法正确性。因为相对定向后，像片相对于参考坐标系的方位元素已知，由地面模型点的坐标利用共线条件方程式计算该点对的左右像点坐标，进而利用式(2)中的方法求得对应的核线影像上点的坐标。根据核线影像上下视差已消除这一特点验证算法的正确性。

表3 模型点坐标及其对应核线影像上下视差

计算这6点上下视差的中误差为Δ=0.00124135pixel。由于像片拍摄的高度约为900m，地面分辨率约为0.12m/pixel，0.12Δ非常小，因此算法误差可以忽略不计，即算法精度完全满足测图需求。

根据核线重采样矩阵生成全幅核线影像的全局用时166737.36ms。

表4是通过原始影像4个角点正算到核线影像上坐标，分别取最大值与最小值，计算在影像快速显示时，将影像按256*256分块的行列分块数。

表4 核线信息

表4中左下点坐标是以像主点为原点的像素坐标。

以(322.54，-10.12)对应的模型点出现在视口内为例来说明立体显示的策略：

假设当前测标Z值为-902.41m，计算点(322.54，-10.12，-902.41)对应在核线影像上的像点坐标及对应块信息如表5所示。

根据当前影像缩放比计算当前影像的金字塔层级。现假设处于第0级，即在原始影像上取值。若当前影像不在第0级时需要将坐标归算到对应层级的影像上。分别将块的4个角点反算至原始影像上，取其外包围矩形内数据，内存中对块逐像素采样。图4中给出了表4信息对应的左右核线影像区块图。区块调度过程中，受限于分块方式及当前测标Z值，同一次调入的左右核线区块范围不一定相同，甚至不会出现相同区域，但当前视口内所有左核线与右核线区块覆盖范围的总和具有重叠区域。

表5 影像块信息

图4 区块采样图像

每个区域的重采样时间约为4.35ms。假设当前窗口的刷新率为30，则完成一帧所用的时间约为33ms；同时刷新所有当前显示在视口内的左右节点影像会造成卡屏，切换不流畅，因此利用线程根据当前位置不断刷新可能出现或可能移出的节点、能很好的保证切换的流畅性。

上面的实验分析是在完成相对定向后进行的处理，模型依据的坐标系是以左片摄站点为原点的参考坐标系。当已知立体影像对的外方位元素时，上述方法与策略也可以进行立体影像对的实时核线影像生成与显示，此时依据的是地面参考坐标系，但是计算的方法、公式以及显示的策略与上面完全一致。

5 结束语

本文针对数字摄影测量工作站中核线影像的生成流程，提出了基于实时核线重采样思想的立体影像调度方法，该方法摒除了传统数字摄影测量系统定向模块中预先生成核线重采样过程，而是在测图显示过程中对影像实时进行核线重采样。因此可简化操作流程，同时减少硬盘存储空间、降低文件的组织复杂度。通过本文的算法实现与验证表明，该方法在精度上满足一般测图目的的要求，在显示时间与速度上在现有的机器配置条件下，满足实时显示、观测的要求，因此方法具有很强的可行性。

本文的研究方法可以为大测区内多像对之间的切换提供策略支持，即可以完成像对间的无缝漫游最终达到无缝测图的目的。文献[12]指出只有网络与计算机集群处理技术充分地应用于新一代的数字摄影测量系统中，才能使数字摄影测量发展到一个新的台阶。本文的研究对基于网络的立体影像显示同样适用，在立体像对的相对位置关系(或者外方位元素已知时)，不需经过核线重采样过程，即可在客户端或者浏览器中查看基于核线影像的立体模型。因此也是影像对网络立体显示的基础之一。

参考文献：

[1] 张祖勋，张剑清.数字摄影测量学[M].武汉：武汉大学出版社，1997：233-255.

[2] 袁修孝，吴珍丽.基于POS数据的核线影像生成方法[J].武汉大学学报(信息科学版)，2008，33(6)：560-564.

[3] 李华光，陈鹰.基于数字纠正的核线影像生成方法[J].山东建筑工程学院学报，2004，19(4)：45-48.

[4] 王宁，刘源璋，吴平，等.基于相对方位元素的核线影像生成方法[J].测绘与空间地理信息，2011，34(3)：254-257.

[5] BJORKMANN M，EKLUNDH J O.Real-time epipolar geometry estimation and disparity[C].The Proceedings of the Seventh IEEE International Conference on.Computer Vision，1999(1)：234-241.

[6] BJORKMAN M，EKLUNDH J O.A real-time system for epipolar geometry and ego-motion estimation[C].Proceedings of IEEE Conference on.Computer Vision and Pattern Recognition，2000(2)：506-513.

[7] CHUM O，WERNER T，MATAS J.Epipolar geometry estimation via ransac benefits from the oriented epipolar constraint[C].Proceedings of the 17th International Conference on Pattern Recognition，2004(1)：112-115.

[8] Toni Schenk.数字摄影测量学——背景、基础、自动定向过程[M].武汉：武汉大学出版社，2009：194-201.

[9] 余粉香，王光霞，万刚，等.大数据量遥感影像的快速调度与显示[J].海洋测绘，2006，26(2)：27-30，40.

[10] 程伟，王波，石剑龙，等.海量图像数据快速显示方法及实现[J].地理与地理信息科学，2011，27(1)：109-110.

[11] 张涵斐，黄忠红，孟永军，等.海量遥感影像的存储与快速调度显示方法[J].测绘与空间地理信息，2011，34(3)：36-39，43.

[12] 张祖勋，张剑清，张力.数字摄影测量发展的机遇与挑战[J].武汉测绘科技大学学报，2000，25(1)：23-27.