CN102663957B - 一种交互式三维城市全景地图的自动生成方法 - Google Patents

Abstract

本发明涉及一种交互式三维城市全景地图的自动生成方法。该方法基于城市形态及格式塔理论，通过建立估计城市空间认知的能量函数偏离项和遮挡项，把道路扩宽、视点提高、建筑位移和建筑高度降低等手段与优化方法结合，利用线性化约束条件的能量函数，实时求解最小化能量函数，从而实现了交互式三维数字城市全景地图的自动生成。该方法提供了一种可视化速度快、城市场景变形小的三维数字城市全景地图的制作方法，在城市导航、空间位置服务等方面具有广阔的应用前景。

Description

一种交互式三维城市全景地图的自动生成方法

一、技术领域

本发明涉及一种交互式三维城市全景地图的自动生成方法，属于空间信息技术领域。

二、背景技术

三维虚拟地理环境在空间位置服务和科学研究等领域具有巨大的应用价值。三维电子导航地图与二维电子导航地图相比，前者对用户来说更加直观、判读轻松、可识别性强，也提供给了第一人称视角。但目前三维电子地图多为透视投影，这种投影是基于对人眼和相机对客观世界的观察模式而设计的，它能真实的在三维空间环境中反映客观世界的形态，缺陷在于，该投影会导致城市中的景观相互遮挡，如用户感兴趣特征(如道路、建筑、湖泊与标注等)常被更靠近视点的物体遮挡，这种遮挡在建筑密集的城市中尤为普遍。

尽管目前有很多方法在一定程度上解决了感兴趣特征被遮挡带来的判读困难，但是仍无法很好的达到三维城市导航地图可读性等问题。例如，在三维电子地图中嵌入二维电子导航图，能够在一定程度上弥补被遮挡的感兴趣区域的信息，但会使用户在三维场景漫游时，同时判读二维电子地图，增加了用户负担并占用了有限的屏幕空间。再如，基于导航语音提示可以令用户预知被遮挡特征的信息，但目前导航信息难以提供复杂情况下的空间信息。

为了避免三维城市电子导航地图中感兴趣特征被遮挡，本发明实现了交互式三维城市全景地图的自动生成方法，在保证交互速率的前提下保证感兴趣特征的可见性，同时尽可能减少了对场景的形变，以保持与原透视投影下城市外观的相似度，使用户更加方便、准确的了解其当前位置、周边环境以及感兴趣区域位置等地理信息，将更多的精力投入到路径查找、驾驶或旅行决策等任务中去。本发明在城市导航、空间位置服务等方面具有广阔的应用前景。

三、发明内容

1、目的：三维城市电子地图中感兴趣特征常被遮挡，无法满足车辆导航的需要，而传统遮挡剔除的方法会对场景造成较大形变。本发明的目的是基于城市形态和格式塔理论，实现了多种遮挡消除策略的能量最小化方法，提供一种交互式三维城市全景地图的自动生成方法，保证感兴趣特征的可见性，避免视点或场景的改变，满足城市导航以及基于空间位置服务等方面的需求。

2、技术方案：

一种交互式三维城市全景地图的自动生成方法，其特征在于，包括如下步骤：

步骤一：基于ε距离邻域和β-骨架规则识别街区内建筑空间分布模式

所要构建DT图的结点为建筑足迹(footprint)多边形的重心，结点之间的距离为建筑足迹多边形间的最小距离；先构造结点的Delaunay三角网，从Delaunay三角网内筛选连接边，建立结点的最小生成树MST；下一步，基于ε距离邻域规则，从Delaunay三角网内选择距离小于阈值ε，且满足β＜1时，把β-骨架规则的非MST边加入子图中，得到满足ε距离邻域条件的DT子图；对于DT子图中度为1的结点，在Delaunay三角网中寻找其满足β＜1时，β-骨架规则最短的非MST边，并将该边加入子图中。

步骤二：提取城市建筑格式塔聚类

设定合理的阈值以找出潜在的格式塔聚类，进而用经训练的支持向量机分类器加以判别；整个建筑的格式塔聚类提取分为2步：(1)确定潜在的建筑格式塔聚类，(2)提取具有格式塔特征的建筑。

步骤三：基于光线跟踪的城市感兴趣特征可见性计算

利用光线跟踪实现城区道路可见性判断。将一条光线和与其相交的遮挡物作为一个遮挡对。每个遮挡对记录有光线、遮挡物及光线击中点参照于遮挡物重心的相对坐标。对于光线与建筑的相交测试，先判断光线在二维上是否相交，快速排除不相交的建筑，再在三维上判断光线是否与建筑的各个面相交。如有交点，则建筑为遮挡物，取距离参考点最近的交点为击中点，否则建筑不造成遮挡；采用动态的KD树结构加快二维平面上光线跟踪。

步骤四：建立城市全景图的约束条件和能量项

(1)建立城市全景图的能量函数偏离项

为了保证格式塔聚类的共同命运特征，建立街区内格式塔聚类位移约束。对道路形状、道路与建筑冲突进行约束，以保证街区与道路的相对位置，此外，约束了建筑高度与视点高度，最后针对漫游时出现的视图抖动与闪烁现象，建立时空一致性约束。

(2)建立城市全景图能量函数的遮挡项

定义x为包含建筑重心与道路条带顶点的二维位置的向量，h为建筑高度向量，对于一个遮挡对{p_r，B_i}，参考点p_r的可见性表示为x的函数。设造成遮挡的建筑B_i高度为h_i，足迹重心为v_b′，则视线在建筑上的击中点为且p_r由其相邻道路顶点表示p_r＝(1-s)v_i0′+sv_(i+1)0′，有可见性函数的二次形式：

$D_{\{p_r,B_i\}}(x,H)\≈{l_{\mathrm{LOS}}}^T\left(x\right)[H(p_{B_i}\left(x\right)-p_r\left(x\right))-h_i(p_v-p_r\left(x\right))]---\left(1\right)$

式中p_v为视点二维位置，H为其高度。

考虑所有遮挡对集合O，将(1)式离散化表示为x，h和H的函数(2)，使感兴趣道路可见。

由此能量函数的遮挡项表示为：

$E_{\mathrm{occ}}(x,h,H)=\left\{\begin{array}{cc}-D(x,h,H),& D(x,h,H)<0\\ 0,& D(x,h,H)\&GreaterEqual;0\end{array}\right.---\left(3\right)$

步骤五：城市全景图能量函数的最优化求解

(1)计算能量函数的最小二乘形式

综合各约束，得到最终的能量函数偏离项

E_dev(x，h，H)＝(E_shp(x)+E_rela(x)+E_ges(x))+η₁E_road(x)+η₂E_dist(x)

(4)

+η₃(E_h1(h)+E_h2(h))+η₄E_H(H)+η₅E_temp(x，h，H)

上述各项都具有正定二次形式，将各项最小二乘形式有：

$E_{\mathrm{dev}}(x,h,H)={\left|\right|\left(\begin{array}{ccc}{A}_{\mathrm{shp}}& & \\ A_{\mathrm{rela}}& & \\ A_{\mathrm{ges}}& & \\ {\widetilde{\mathrm{\&eta;}}}_1{A}_{\mathrm{road}}& & \\ {\widetilde{\mathrm{\&eta;}}}_2{A}_{\mathrm{dist}}& & \\ {\widetilde{\mathrm{\&eta;}}}_5{A}_{\mathrm{temp}1}& & \\ {\widetilde{\mathrm{\&eta;}}}_5{A}_{\mathrm{temp}2}& & \\ & {\widetilde{\mathrm{\&eta;}}}_3{B}_{h1}& \\ & {\widetilde{\mathrm{\&eta;}}}_3{B}_{h2}& \\ & {\widetilde{\mathrm{\&eta;}}}_{5}I& \\ & & {\widetilde{\mathrm{\&eta;}}}_4\\ & & {\widetilde{\mathrm{\&eta;}}}_5\end{array}\right)\left(\begin{array}{c}x\\ h\\ H\end{array}\right)\left(\begin{array}{c}{b}_{\mathrm{shp}}\\ b_{\mathrm{rela}}\\ b_{\mathrm{ges}}\\ {\widetilde{\mathrm{\&eta;}}}_1{b}_{\mathrm{road}}\\ {\widetilde{\mathrm{\&eta;}}}_2{b}_{\mathrm{dist}}\\ {\widetilde{\mathrm{\&eta;}}}_5{b}_{\mathrm{temp}1}\\ {\widetilde{\mathrm{\&eta;}}}_5{b}_{\mathrm{temp}2}\\ {\widetilde{\mathrm{\&eta;}}}_3{q}_{h1}\\ {\widetilde{\mathrm{\&eta;}}}_3{q}_{h1}\\ {\widetilde{\mathrm{\&eta;}}}_5\\ {\widetilde{\mathrm{\&eta;}}}_4\\ {\widetilde{\mathrm{\&eta;}}}_5\end{array}\right)\left|\right|}^2---\left(5\right)$

其中A_(·)和B_(·)分别对应各能量项与二维位置x和建筑高度h高度相关的矩阵，而b_(·)和q_(·)则是与对应各能量项与二维位置x和建筑高度h高度相关的向量，而表示原权重系数η_(·)的平方根。

(2)数值求解

1)能量函数的线性化

对于E_occ(x，h，H)，求解能最小化E_occ(x，h，H)消除遮挡。约束为令各遮挡对的可见性函数 $D_{\{p_r,B_i\}}(x,h,H)\&GreaterEqual;0,$ 有：

$\stackrel{\&OverBar;}{x}=\arg \underset{\stackrel{\&OverBar;}{x}\&Element;P}{\min }{\stackrel{\&OverBar;}{x}}^{T}H\stackrel{\&OverBar;}{x}-2{\stackrel{\&OverBar;}{x}}^T{\stackrel{\&OverBar;}{A}}^{T}d+d^{T}d$

(6)

式(6)线性化后

$s.t.G\left(\begin{array}{c}x\\ h\end{array}\right)\&GreaterEqual;c---\left(7\right)$

$A=\left(\begin{array}{c}{A}_{\mathrm{shp}}\\ A_{\mathrm{rela}}\\ A_{\mathrm{ges}}\\ {\widetilde{\mathrm{\&eta;}}}_1{A}_{\mathrm{road}}\\ {\widetilde{\mathrm{\&eta;}}}_2{A}_{\mathrm{dist}}\\ {\widetilde{\mathrm{\&eta;}}}_5{A}_{\mathrm{temp}1}\\ {\widetilde{\mathrm{\&eta;}}}_5{A}_{\mathrm{temp}2}\end{array}\right),b=\left(\begin{array}{c}{b}_{\mathrm{shp}}\\ b_{\mathrm{rela}}\\ b_{\mathrm{ges}}\\ {\widetilde{\mathrm{\&eta;}}}_1{b}_{\mathrm{road}}\\ {\widetilde{\mathrm{\&eta;}}}_2{b}_{\mathrm{dist}}\\ {\widetilde{\mathrm{\&eta;}}}_5{b}_{\mathrm{temp}1}\\ {\widetilde{\mathrm{\&eta;}}}_5{b}_{\mathrm{temp}2}\end{array}\right),B=\left(\begin{array}{c}{\widetilde{\mathrm{\&eta;}}}_3{B}_{h1}\\ {\widetilde{\mathrm{\&eta;}}}_3{B}_{h2}\\ {\widetilde{\mathrm{\&eta;}}}_{5}I\end{array}\right),q=\left(\begin{array}{c}{\widetilde{\mathrm{\&eta;}}}_3{q}_{h1}\\ {\widetilde{\mathrm{\&eta;}}}_3{q}_{h1}\\ {\widetilde{\mathrm{\&eta;}}}_4\end{array}\right).$

2)用Kalman滤波方法实现视点高度的平滑变化

为了满足视点的时间一致性，使用Kalman滤波方法令视点高度平滑变化。

3、优点及功效：针对三维城市场景复杂遮挡的处理，以前的方法只是简单提高视点或扩宽道路，引入了较大变形保证可见性。本发明提出了结合多种遮挡消除策略的能量最小化方法，不仅保证了感兴趣特征的可见性，而且避免了对视点或场景特征大幅度的改变。所设计的DT子图能有效发掘街区内建筑群的排列模式；设计了一系列能量函数保证了建筑、道路的相似性并维持了时间一致性；在求解时采用了带约束能量函数，约束求解时可用预计算矩阵分解及矩阵系数沿对角线块状独立分布的特性，达到快速求解的目的；能够自动生成交互式三维城市全景地图，满足城市导航的需要。

四、附图说明

图1三维城市全景地图自动生成方法的流程示意图

图2道路的条带表示示意图

图3建筑与道路间约束关系示意图

图4(a)原视图

图4(b)本发明得到的全景地图

图5(a)原视图

图5(b)本发明得到的全景地图

图6(a)原视图

图6(b)本发明得到的全景地图

图7(a)原视图

图7(b)本发明得到的全景地图

五、具体实施方式

本发明涉及一种交互式三维城市全景地图的自动生成方法，该方法的具体步骤如下：

步骤一：基于ε距离邻域和β-骨架规则识别街区内建筑空间分布模式

本发明利用介于Delaunay三角网(DelaunayTriangulation)与最小生成树之间的数据类型DT子图识别街区内建筑空间分布模式。为此引入ε距离邻域规则：在Delaunay三角网各条边中，若其距离小于某指定阈值ε，则该边予以保留，作为其两端结点i，j之间存在连结边。通过此规则构造的DT子图将保留局部边连接，保持局部的拓扑特性。

但是引入ε距离邻域规则得到的图会存在两个缺点：一是可能产生奇异三角形，即存在两个较小锐角的三角形，而奇异三角形容易在优化时导致数值问题；另外是对于部分孤立的结点有可能出现一条连结边的情形，即结点的度为1，导致其移动自由度过大，进而在位移时对与其邻近的结点造成较大的扰动。对于前者，引入另一种图类型β-骨架在β＜1时的构图规则，用于消除这种奇异三角形：对于平面上任意两结点i，j，若在邻域B内存在其它结点，则结点i，j之间具有连结边，此处邻域B定义为两个半径为d(i，j)/β且过结点i，j的圆的交集，避免了存在其他结点k使得夹角∠ikj大于某特定值θ时，i，j相连结的情形，当β＜1时，此处

θ＝π-arcsinβ，β∈(0，1).(1)

因此，引入β＜1时，β-骨架规则可有效避免奇异三角形的产生。

对于孤立的结点可能会出现一条连结边的情形，本发明引入流形学习中构建图的k最近邻规则，即，对于每个结点至少与其最近的k个结点相连接。此处k＝2，显然，结点k最近邻在必定在结点在Delaunay三角网的相邻结点中，得到的图为DT子图。

下面构造DT子图。首先定义所要构建图的结点为建筑足迹(footprint)多边形的重心，结点之间的距离为建筑足迹多边形间的最小距离。先构造结点的Delaunay三角网，从Delaunay三角网内筛选连接边，基于Delaunay三角网建立结点的最小生成树MST。下一步，基于ε距离邻域规则，从Delaunay三角网内选择距离小于阈值ε，且满足β＜1时，β-骨架规则的非MST边加入子图中，得到满足ε距离邻域条件的DT子图。对于DT子图中度为1的结点，在Delaunay三角网中寻找其满足β＜1时，β-骨架规则最短的非MST边，并将该边加入子图中。至此，DT子图便构造完成。一般的，ε＝75米，

步骤二：提取城市建筑格式塔聚类

目前，对格式塔聚类的分类都是采用人工设定阈值，将满足格式塔聚类与不满足格式塔聚类区分开的方法。本发明采用机器学习算法，将视觉上直观的判别，转化为数值上人为难以确定的特征判别函数。首先设定合理的阈值以找出潜在的格式塔聚类，进而用经训练的支持向量机分类器加以判别。整个建筑的格式塔聚类提取分为如下2步。

(1)确定潜在的建筑格式塔聚类

将一个潜在格式塔聚类G定义为建筑B_i的集合：G＝{B₀，B₁，...B_n}。利用种子填充算法找出潜在聚类，以任意建筑B₀为原点开始，初始时G＝{B₀}，沿着某可行方向依次将建筑B₁，...B_n加入G中。每次要加入新建筑时，只考虑与G中最后一幢建筑B_n相邻近的建筑，判断其是否符合当前聚类G的格式塔特征，如符合则作为第n+1个元素加入G中，否则终止对聚类G的扩大。聚类中建筑B_i只与B_i-1和B_i+1存在空间邻近关系，在聚类时，有可能出现与建筑B_i-1邻近的k个建筑都符合当前聚类G的格式塔特征，此时，得到k个新的潜在聚类G：G∪B_i，再分别继续对k潜在聚类进行扩张。这样从B₀扩张将有可能得到m个潜在聚类，这m个潜在聚类之间也可能两两符合同一类格式塔特征，于是从B₀处将这两个聚类合并。

将新建筑B_i加入潜在格式塔聚类时，只考虑邻近性、相似性与共向性。对于邻近性特征，用建筑足迹多边形顶点构建Delaunay三角网，限制足迹多边形上边必须是Delaunay三角网的边。如果三角网中存在顶点属于两幢不同建筑的三角形，则认为建筑可通视，即建筑有空间邻近关系。

对于相似性特征，则用面积一致性做粗略地近似。建筑面积相似性可用B_i面积Area(B_i)与聚类中建筑平均面积的比值R_area表示：

$R_{\mathrm{area}}=\frac{\mathrm{Area}\left(B_i\right)}{\stackrel{\&OverBar;}{A}}---\left(2\right)$

只要R_area满足R_area∈(1/ε_area，ε_area)，则认为B_i和B_i+1可能相似，ε_area＝2.0。

对于共向性特征，可分别考虑建筑主方向一致性以及建筑线状延伸一致性。前者主要是指单幢建筑的主要朝向，其计算方法将下面给出。主要朝向一致性同样用当前建筑B_i的主朝向与聚类中建筑朝向的均值的夹角θ_Orient表示：

${\mathrm{\&theta;}}_{\mathrm{Orient}}=\arccos \left(\frac{O_{B_i}\&CenterDot;\stackrel{\&OverBar;}{O}}{\left|\right|O_{B_i}\left|\right|\left|\right|\stackrel{\&OverBar;}{O}\left|\right|}\right).---\left(3\right)$

只要θ_orient＜ε_orient，则B_i满足朝向一致性，这里ε_orient设为30°。

对于线状方向一致性，则计算聚类中初始建筑B₀的重心到当前建筑B_i的重心的方向d_i，与首幢建筑B_i的重心到聚类中上一幢建筑B_i-1的重心方向d_i-1，对于d_i与d_i-1间夹角θ_dir，只要有θ_dir＜ε_dir则认为B_i满足线状方向一致性，这里ε_dir设为15°。

在考虑两个潜在聚类G_i与G_j的合并时，同样采用与加入新建筑相似的判断方法。对于面积一致性与朝向一致性，采用与(2)和(3)式相同的形式与阈值比较两类的平均面积与朝向。关键是检查线状方向一致性，设G_i与G_j的第二幢建筑分别为B_1i和B_1j，当B_1i和B_1j并非同幢建筑时，检查B₀重心到B_1i重心的方向d_i与B_1j重心到B₀重心的方向d_j之间夹角θ_dir，同样若θ_dir＜ε_dir则认为满足线状方向一致性。

(2)提取具有格式塔特征的建筑

本发明基于以下观察：建筑墙面在与其本身朝向重合的两个建筑朝向方向上具有最小投影，即建筑的主要墙面一般决定了建筑朝向。而当墙面与两个建筑朝向方向都成45°夹角时，累加投影长度达到最大，这正是要计算的方向。因此，目标函数可以通过对墙面在两个朝向方向上的投影做线积分，即求和得到。设主朝向方向为d，可最小化如下目标函数

$d=\arg \underset{d}{\min }\underset{p_i\&Element;F}{\mathrm{\&Sigma;}}\left(\right|p_i\&CenterDot;d|+|p_i\&CenterDot;d^{\&perp;}\left|\right).---\left(4\right)$

其中p_i表示建筑足迹线段矢量，F为建筑足迹，d^⊥表示与d相垂直的次朝向。对于(4)式，可采用Newton法求解。

下面给出计算具有n幢建筑的潜在聚类G各项特征的方法：

(1)面积差异：设潜在聚类G中建筑的面积均值为标准差为σ_S，面积差异δ_area可量化表达为 ${\mathrm{\&delta;}}_s={\mathrm{\&sigma;}}_s/\stackrel{\&OverBar;}{S}.$

(2)高度差异：设潜在聚类G中建筑的高度均值为标准差为σ_h，高度差异δ_height可量化表达为 ${\mathrm{\&delta;}}_{\mathrm{height}}={\mathrm{\&sigma;}}_h/\stackrel{\&OverBar;}{h}.$

(3)相似差异：相似性是一种难以量化的指标。这里采取ShapeContext的方法计算两幢建筑B_i与B_j足迹之间的相似度s_ij。其中每幢建筑B_i都与聚类中其余建筑B_j得到一个相似度，将这些相似度的均值定义为建筑B_i的相似度：

$s_i=\frac{1}{n-1}\underset{i\&NotEqual;j}{\mathrm{\&Sigma;}}{s}_{\mathrm{ij}}.---\left(5\right)$

此处取最小建筑相似度作为聚类的相似差异δ_simil：δ_simil＝1-min_i∈Gs_i，这是为了避免G中存在与其它建筑足迹极其不相似的建筑。

(4)朝向差异：直接采用各建筑主朝向方位角θ的标准差σ_θ作为朝向差异δ_orient，方位角θ被规则化在[0，π)区间内。

(5)线状排列不一致性：线状不一致性采用建筑重心之间连线的折角来确定。设有建筑B_i，1＜i＜n，B_i重心处的折角可表示为B_i重心到B_i-1重心处的方向-d_i-1与B_i重心到B_i+1重心处的方向d_i-1之间的夹角α_i，则聚类G中共有n-1个这样的折角。对于B_i的排列不一致性v_i可表示为v_i＝π-α_i。则聚类G的线状排列不一致性δ_align用v_i的均值定义。

(6)间隔差异：设L_ij表示相邻两幢建筑B_i与B_j之间的最小距离，显然聚类G中可以得到n-1个这样的距离。计算距离L_ij的均值与方差σ_L，则间隔差异δ_interv可量化为

(7)排列倾斜度：排列倾斜度定义为G内建筑排列方向与G内建筑主要朝向的夹角大小。建筑的排列方向可用G内建筑重心的拟合直线方向表示，设p为归一化的拟合直线方向，(d，d^⊥)为一对归一化的建筑朝向，排列倾斜度可表示为

δ_incline＝max(p·d，p·d^⊥)(6)

组合这上述7个特征得到了聚类G的特征向量。对于训练样本的特征向量x，利用支持向量机学习了判别函数D(x)。这里判别函数的符号D(x)表示了类别，其大小则表示了距离分类面的远近程度，对于正类，即格式塔聚类而言，距离分类面越远代表它符合格式塔聚类的程度越高。

对于潜在格式塔聚类G，尽管其本身可能未能被判别为格式塔聚类，但是有可能存在其子集0≤c＜c+m-1≤n，m≥3为格式塔聚类的情形。因此，需要按规模从大到小判定其各子集，直到判别出格式塔聚类为止。对于同规模子集同时被判别为格式塔聚类时，取判别函数最大值D(x)的子集为格式塔聚类，对于剩下的补集，继续递归地寻找潜在的格式塔聚类。

步骤三：基于光线跟踪的城市感兴趣特征可见性计算

利用光线跟踪实现城区道路可见性判断。将一条光线和与相交的遮挡物作为一个遮挡对。每个遮挡对记录有光线、遮挡物及光线击中点参照于遮挡物重心的相对坐标。

对于光线与建筑的相交测试，可先判断光线在二维上是否相交，快速排除不相交的建筑，再在三维上判断光线是否与建筑的各个面相交。如有交点，则建筑为遮挡物，取距离参考点最近的交点为击中点，否则建筑不造成遮挡。采用动态的KD树结构加快二维平面上光线跟踪。因为建筑是按街区分割的，而且建筑位移不会离开基本静态的街区，所以以街区为单位进行相交测试。首先判断光线是否与街区相交，若建筑与街区包围盒相交，则再判断建筑是否与街区多边形相交。而在街区内，虽然建筑个数较少，为了避免对街区所有建筑足迹做相交测试，可以建立另一层加速结构，本发明采用了排序表加速结构。对二维平面的两个轴建立链表，链表的每个结点记录建筑足迹多边形包围盒AABB(AxisAlignedBoundingBox)在当前轴上的最小或最大值。在光线跟踪时，根据光线在当前街区的进入点，从该坐标轴对应链表中找到相应结点，从该结点开始，判断光线是否与AABB包围盒相交，直到坐标值大于光线在街区上出口坐标值的结点为止。

步骤四：建立城市全景图的约束条件和能量项

(1)建立城市全景图的能量函数偏离项

在城市空间认知基础上，保持视图相似性的线性约束。从单个街区出发，对街区的形状进行约束，以保证建筑的相对位置与相对高度。为了保证格式塔聚类的共同命运特征，建立街区内格式塔聚类位移约束。对道路形状、道路与建筑冲突进行约束，以保证街区与道路的相对位置，此外，约束了建筑高度与视点高度，最后针对漫游时出现的视图抖动与闪烁现象，建立时空一致性约束。

1)约束城市街区内建筑的分布模式

为了位移时尽量保持反映街区内建筑分布模式的DT子图的形状，要求局部的位移变化是平滑、局部一致的。如果DT子图上任何位置都满足此特性，则将保持DT子图的全局形状。

设建筑位移函数为d(x)：□²→□²，是二维平面位置x的函数，令d(x)的变化，即Jacobian为0：

$J_{d\left(x\right)}=\frac{\&PartialD;d}{\&PartialD;x}=0.---\left(7\right)$

设d(x)＝(d_x(x)，d_y(x))^T，将(7)展开：

$J_{d\left(x\right)}=\left(\begin{array}{cc}\frac{\&PartialD;d_x\left(x\right)}{\&PartialD;x}& \frac{\&PartialD;d_x\left(x\right)}{\&PartialD;y}\\ \frac{\&PartialD;d_y\left(x\right)}{\&PartialD;x}& \frac{\&PartialD;d_y\left(x\right)}{\&PartialD;y}\end{array}\right)=0.---\left(8\right)$

为了令整个DT子图有此特性，即在整个二维域Ω上的位移d(x)尽可能满足(7)式，求解如下最小泛函d(x)：

$E\left(d\right)=\underset{\mathrm{\&Omega;}}{\&Integral;}{\left|\right|J_d\left|\right|}_F^2\mathrm{dx}.---\left(9\right)$

此处，||·||_F为矩阵Frobenius范数，展开有，

$E\left(d\right)=\underset{\mathrm{\&Omega;}}{\&Integral;}{\left|\right|\frac{\&PartialD;d}{\&PartialD;x}\left|\right|}^2+{\left|\right|\frac{\&PartialD;d}{\&PartialD;y}\left|\right|}^2\mathrm{dxdy}.---\left(10\right)$

求解得到，

Ad＝0.(11)

此处，Δ为Laplace算子。对于图域的离散化，采用流形学习中权重图上的离散化方法：

$\underset{j\&Element;N_i}{\mathrm{\&Sigma;}}{w}_{\mathrm{ij}}(d_j-d_i)=0.---\left(12\right)$

N_i为i的邻域结点集合，w_ij为连接结点i和j的边的权重，这里采用了热核函数，近似为，

$w_{\mathrm{ij}}=\exp (-\frac{{\left|\right|x_j-x_i\left|\right|}^2}{t}).---\left(13\right)$

对于DT子图，也近似采用(13)式作为离散化时的权重。

因为要得到关于结点位移位置的能量函数，设对于结点i，位移前后的位置分别为v_i和v_i′

v_i′＝v_i+d_i.(14)

代入(12)式后有

$\underset{j\&Element;N_i}{\mathrm{\&Sigma;}}{w}_{\mathrm{ij}}({v_i}^{\&prime;}-{v_j}^{\&prime;})=\underset{j\&Element;N_i}{\mathrm{\&Sigma;}}{w}_{\mathrm{ij}}(v_i-v_j).---\left(15\right)$

为了令(15)式成立，对于v_i′的能量项e(v_i′)，设l_ij＝v_i-v_j、l_ij′v_i′-v_j′有

$e_{\mathrm{shp}\left\{v_i^{\&prime;}\right\}}=\underset{j\&Element;N_i}{\mathrm{\&Sigma;}}{w}_{\mathrm{ij}}{\left|\right|{l_{\mathrm{ij}}}^{\&prime;}-l_{\mathrm{ij}}\left|\right|}^2.---\left(16\right)$

令e(v_i′)＝0以最小化e(v_i′)后，得到(15)式。从(16)式看出，能量项e(v_i′)具有保持连接边长度特性，从而在某种程度上避免位移时建筑间的碰撞。对建筑集合设x为包含各结点位置v_i′的向量，得到建筑集合的能量为，

2)约束建筑的相对位置

为了避免相对位移，要保证局部单个建筑与邻近建筑的相对位置，所以只需考虑DT子图中满足ε距离邻近的连接边。对于单个建筑，即DT子图中的结点i满足ε距离邻近的结点集合为N_ε(i)，要求它们之间的相对位移为0，得到

$\underset{j\&Element;N_{\mathrm{\&epsiv;}}\left(i\right)}{\mathrm{\&Sigma;}}(d_j-d_i)=0.---\left(18\right)$

建立关于v_i′的能量函数，将(14)式代入(18)式后有，

$\underset{j\&Element;N_{\mathrm{\&epsiv;}}\left(i\right)}{\mathrm{\&Sigma;}}{v_i}^{\&prime;}-{v_j}^{\&prime;}=\underset{j\&Element;N_{\mathrm{\&epsiv;}}\left(i\right)}{\mathrm{\&Sigma;}}{v}_i-v_j.---\left(19\right)$

(19)式等号右边为常数，可看作结点i的关于其邻近结点相对位置，

$\underset{j\&Element;N_{\mathrm{\&epsiv;}}\left(i\right)}{\mathrm{\&Sigma;}}{v}_i-v_j={\mathrm{\&delta;}}_i.---\left(20\right)$

(20)式试图保持局部的几何特征，希望结点i位移后的相对位置δ_i′与原始相对位置一致，得到能量项：

$e_{\mathrm{rela}\left\{v_i\right\}}=\underset{j\&Element;N_{\mathrm{\&epsiv;i}}}{\mathrm{\&Sigma;}}{\left|\right|{{\mathrm{\&delta;}}_i}^{\&prime;}-{\mathrm{\&delta;}}_i\left|\right|}^2.---\left(21\right)$

考虑到所有建筑有

ω_i为结点i与相邻结点相对位置的权重。一般的，周边建筑距离较小的建筑具有较大的权重。权重ω_i与周边建筑平均距离相关，因此，

${\mathrm{\&omega;}}_i\&Proportional;\exp (-\frac{{\stackrel{\&OverBar;}{d}}_i^2}{t}),{\stackrel{\&OverBar;}{d}}_i=\frac{1}{\left|N_{\mathrm{\&epsiv;}}\left(i\right)\right|}\underset{j\&Element;N_{\mathrm{\&epsiv;}}\left(i\right)}{\mathrm{\&Sigma;}}d(B_i,B_j).---\left(23\right)$

d(B_i，B_j)表示建筑B_i与B_j之间最近距离。

3)约束建筑格式塔共同命运

对于同一格式塔聚类g_i内建筑，使其按共同方向t∈□²位移，可最小化如下能量：

$e_{\mathrm{ges}\left\{g_i\right\}}=\underset{j\&Element;g_i}{\mathrm{\&Sigma;}}{\left|\right|v_i+t-{v_i}^{\&prime;}\left|\right|}^2.---\left(24\right)$

对格式塔聚类集合内所有元素，则要求其满足最小化(23)式，有

4)道路约束

道路的约束包括道路位置和形状两方面。首先将道路几何形状定义为条带形，即由沿道路中轴线{v₁，v₂，...，v_n}两边的顶点{v₁₀，v₁₁，v₂₀，v₂₁，...，v_n0，v_n1}构成，如图2所示。定义顶点下缀的首位表示道路顶点的序号，第二位表示顶点在道路的左手边或右手边。显然有

$\frac{v_{i0}+v_{i1}}{2}\&equiv;v_i.---\left(26\right)$

这里只改变条带上的顶点v_i0和v_i1，而视v₁为常量。

为了不移动道路的中轴线，建立如下能量函数约束道路位置。v_i为道路中轴线位置，有：

$e_{\mathrm{cen}\{v_{i0},v_{i1}\}}={\left|\right|\frac{{v_{i0}}^{\&prime;}+{v_{i1}}^{\&prime;}}{2}-v_i\left|\right|}^2.---\left(27\right)$

尽管要扩宽感兴趣的道路，或要在建筑位移时压缩次要道路，但仍要求道路有一定“弹性”以保持其原始宽度。设道路原始宽度为w，有：

$e_{\mathrm{wid}\{v_{i0},v_{i1}\}}={\left|\right|({v_{i0}}^{\&prime;}-{v_{i1}}^{\&prime;})-w\frac{v_{i0}-v_{i1}}{\left|\right|v_{i0}-v_{i1}\left|\right|}\left|\right|}^2.---\left(28\right)$

为了使道路形状不发生改变，需要保持道路条带上的线段与中轴线平行。设道路的一节中轴线段为v_i，v_i+1，与其平行的直线满足直线方程

n_i ^Tp+c＝0，p＝(x，y)^T.(29)

n_i为与垂直的单位矢量，c为与中轴线距离相关的常数。以中轴线左手边的条带线段v_i0v_(i+1)0为例，当其满足(25)式时，有

n_i ^T(v_(i+1)0′-v_i0′)＝0.(30)

最小化如下能量，

$e_{\mathrm{para}\{v_{i0},v_{(i+1)0}\}}={\left|\right|{n_i}^T({v_{(i+1)0}}^{\&prime;}-{v_{i0}}^{\&prime;})\left|\right|}^2.---\left(31\right)$

设所有道路中轴线顶点个数为n，所有道路条带顶点的集合为S，综合(27)、(28)和(31)得到道路约束能量有，

$E_{\mathrm{road}}\left(x\right)=\underset{i}{\overset{n}{\mathrm{\&Sigma;}}}(e_{\mathrm{cen}\{v_{i0},v_{i1}\}}\left(x\right)+e_{\mathrm{wid}\{v_{i0},v_{i1}\}}\left(x\right)+\mathrm{\&gamma;}\underset{\{v_{i\&CenterDot;},v_{(i+1)\&CenterDot;}\}\&Element;S}{\mathrm{\&Sigma;}}{e}_{\mathrm{para}\{v_i\&CenterDot;,v_{(i+1)\&CenterDot;}\}}\left(x\right)).---\left(32\right)$

其中γ＝10.0。

5)约束建筑与道路之间的距离

为了令街区内建筑与街区周边所有道路都产生相互作用，对于街区周边任意一条道路，必须至少建立一个距离约束。若该道路没有处于最近距离ε范围的建筑，则取离该道路最近的建筑建立距离约束。如图3所示，设与道路相邻的建筑为B_i，其足迹的重心位置为v_b，与其相邻的道路线段为v_i0v_(i+1)0，两者间最近距离为d。建筑足迹上与道路最近点为p_b＝v_b+t，道路上最近点为p_r＝(1-s)v_i0+sv_(i+1)0，要求位移后能保持最近距离d，有，

$e_{\mathrm{dist}\{v_b,v_{i0},v_{(i+1)0}\}}={\left(\right||{v_b}^{\&CenterDot;}+t-(1-s){v_{i0}}^{\&prime;}-{{\mathrm{sv}}_{(i+1)0}}^{\&prime;}||-d)}^2.---\left(33\right)$

设所有满足如上邻近关系集合为则能量函数为，

6)约束建筑高度

设建筑高度为h_i，h是描述所有建筑高度的向量，对于建筑集合有，

为表示建筑重要性的权重，在[0.5，2]区间范围内。

为了保持邻近建筑相对高度，定义建筑i与邻近建筑相对高度为：

$\underset{j\&Element;N_{\mathrm{\&epsiv;i}}}{\mathrm{\&Sigma;}}{h}_i-h_j={\mathrm{\&delta;}}_{\mathrm{hi}}.---\left(36\right)$

对令相对高度一致有，

其中ω_i为结点i与相邻结点相对位置的权重。

7)时间一致性约束

为了保持建筑位移在时间上平滑变化的，建立建筑位移的时间一致性约束。以上标(·)^t表示t时刻的状态，为此令DT子图G上的结点在当前帧的位移变化与上一帧一致，有：

$E\left(d^t\right)=\underset{\mathrm{\&Omega;}}{\&Integral;}{\left|\right|J_d^t-J_d^{t-1}\left|\right|}_F^2\mathrm{dx}.---\left(38\right)$

其中所表示的上一帧位移变化是已知的，求解泛函(37)式得到一个Poisson方程：

Δd^t＝Δd^t-1.(38)

根据(12)至(14)式对等号两边离散化后有

$\underset{j\&Element;N_i}{\mathrm{\&Sigma;}}{w}_{\mathrm{ij}}({v_i}^{\&prime;t}-{v_j}^{\&prime;t})=\underset{j\&Element;N_i}{\mathrm{\&Sigma;}}{w}_{\mathrm{ij}}({v_i}^{\&prime;t-1}-{v_j}^{\&prime;t-1}).---\left(39\right)$

设l_ij′^t＝v_i′^t-v_j′^t，

$e_{\mathrm{temp}1\left\{v_i^{\&prime;}\right\}}=\underset{j\&Element;N_i}{\mathrm{\&Sigma;}}{w}_{\mathrm{ij}}{\left|\right|{l_{\mathrm{ij}}}^{\&prime;t}-{l_{\mathrm{ij}}}^{\&prime;t-1}\left|\right|}^2.---\left(40\right)$

为了保持相邻帧之间道路的宽度，有：

$e_{\mathrm{temp}2\left\{v_{i\&CenterDot;}^{\&prime;t}\right\}}={\left|\right|{v_{i\&CenterDot;}}^{\&prime;t}-{v_{i\&CenterDot;}}^{\&prime;t-1}\left|\right|}^2.---\left(41\right)$

为了保持建筑高度，令当前所有建筑的高度与上一帧一致，有：

$e_{\mathrm{temp}3\left\{h_i^{\&prime;t}\right\}}={({h_i}^{\&prime;t}-{h_i}^{\&prime;t-1})}^2.---\left(42\right)$

最后，令视点高度与上一帧一致有，

e_temp4＝(H′^t-H′^t-1)².(43)

综合(40)至(43)式，对所有建筑的集合有即道路条带上的顶点S，时间一致性约束表示为：

(2)建立城市全景图能量函数的遮挡项

可见性约束要涉及建筑与道路顶点的位置，建筑高度以及视点位置。定义x为包含建筑重心与道路条带顶点的二维位置的向量，h为建筑高度向量。对于一个遮挡对{p_r，B_i}，参考点p_r的可见性可表示为x的函数。设造成遮挡的建筑B_i高度为h_i，足迹重心为v_b′，为p_r、之间的向量，则视线在建筑上的击中点为且p_r可由其相邻道路顶点表示p_r＝(1-s)v_i0′+sv_(i+1)0′，有

$D_{\{p_r,B_i\}}(x,H)=H\left|\right|p_{B_i}\left(x\right)-p_r\left(x\right)\left|\right|-h_i\left|\right|p_v-p_r\left(x\right)\left|\right|---\left(45\right)$

式中p_v为视点二维位置，H为其高度。不妨假设平移后的p_h(x)仍处于视线p_r(x)p_v附近，为此对(45)乘以

$\frac{{(p_{B_i}\left(x\right)-p_r\left(x\right))}^T(p_v-p_r\left(x\right))}{\left|\right|p_h\left(x\right)-p_r\left(x\right)\left|\right|\left|\right|p_v-p_r\left(x\right)\left|\right|}$

得到

$D_{\{p_r,B_i\}}(x,H)=H\frac{{(p_v-p_r\left(x\right))}^T(p_{B_i}\left(x\right)-p_r\left(x\right))}{\left|\right|p_v-p_r\left(x\right)\left|\right|}-h_i\frac{{(p_{B_i}\left(x\right)-p_r\left(x\right))}^T(p_v-p_r\left(x\right))}{\left|\right|p_{B_i}\left(x\right)-p_r\left(x\right)\left|\right|}---\left(46\right)$

让p_h(x)-p_r(x)近似与(p_v-p_r(x))平行，有

$\frac{p_{B_i}\left(x\right)-p_r\left(x\right)}{\left|\right|p_{B_i}\left(x\right)-p_r\left(x\right)\left|\right|}\&ap;\frac{p_v-p_r\left(x\right)}{\left|\right|p_v-p_r\left(x\right)\left|\right|}=l_{\mathrm{LOS}}\left(x\right).$

最终近似得到可见性函数(45)式的二次形式：

$D_{\{p_r,B_i\}}(x,H)\&ap;{l_{\mathrm{LOS}}}^T\left(x\right)[H(p_{B_i}\left(x\right)-p_r\left(x\right))-h_i(p_v-p_r\left(x\right))]---\left(48\right)$

考虑所有遮挡对集合将(48)式离散化表示为x，h和H的函数(49)，使感兴趣道路可见。

由此能量函数的遮挡项可表示为：

$E_{\mathrm{occ}}(x,h,H)=\left\{\begin{array}{cc}-D(x,h,H),& D(x,h,H)<0\\ 0,& D(x,h,H)\&GreaterEqual;0\end{array}\right..---\left(50\right)$

步骤五：城市全景图能量函数的最优化求解

(1)计算能量函数的最小二乘形式

先求解位移方向t，令(24)式为t的函数，有，

$t=\arg \mathrm{mi}{n}_t\underset{i\&Element;g_i}{\mathrm{\&Sigma;}}{\left|\right|v_i+t-{v_i}^{\&prime;}\left|\right|}^2---\left(51\right)$

得到

$t=\frac{1}{\left|g_i\right|}\underset{j\&Element;g_i}{\mathrm{\&Sigma;}}{v_j}^{\&prime;}-v_j.---\left(52\right)$

将(52)式代入(24)式得到有正定二乘形式的能量项

$\begin{array}{c}{e}_{\mathrm{ges}\left\{g_i\right\}}=\underset{i\&Element;g_i}{\mathrm{\&Sigma;}}{\left|\right|(v_i-\frac{1}{\left|g_i\right|}\underset{j\&Element;g_i}{\mathrm{\&Sigma;}}{v}_j)-({v_i}^{\&prime;}-\frac{1}{\left|g_i\right|}\underset{j\&Element;g_i}{\mathrm{\&Sigma;}}{v_j}^{\&prime;})\left|\right|}^2\\ =\frac{1}{\left|g_i\right|}\underset{i\&Element;g_i}{\mathrm{\&Sigma;}}{\left|\right|\underset{j\&Element;g_i,j\&NotEqual;i}{\mathrm{\&Sigma;}}(v_i-v_j)-\underset{j\&Element;g_i,j\&NotEqual;i}{\mathrm{\&Sigma;}}\left|\right|}^2\end{array}.---\left(53\right)$

(53)式可以保持结点的相对位置。

对于(33)式，设过线段v_i0′v_(i+1)0′的直线方程为

n_i0′^Tp+c_i0′＝0，p＝(x，y)^T.(54)

此处n_i0′为单位法向量。则点p′到线段v_i0′v_(i+1)0′的距离为n_i0′^Tp+c′。因为(31)式约束了线段v_i0′v_(i+1)0′的方向，可近似认为n_i0′为常数且等于(31)中的n_i，则c_i0′为c_i0′＝-n_i ^Tv_i0，最终得到(33)式的二次函数近似：

$e_{\mathrm{dist}\{v_b,v_{i0},v_{(i+1)0}\}}\&ap;{({n_i}^T({v_b}^{\&prime;}+t-{v_{i0}}^{\&prime;})-d)}^2.---\left(55\right)$

综合各约束，得到最终的能量函数偏离项

E_dev(x，h，H)＝(E_shp(x)+E_rela(x)+E_ges(x))+η₁E_road(x)+η₂E_dist(x)

+η₃(E_h1(h)+E_h2(h))+η₄E_H(H)+η₅E_temp(x，h，H)

上述各项都具有正定二次形式。实验中，η₁、η₂和η₃被设置为0.1，而η₄和则η₅被设置为1.0。将各项最小二乘形式有：

$E_{\mathrm{dev}}(x,h,H)={\left|\right|\left(\begin{array}{ccc}{A}_{\mathrm{shp}}& & \\ A_{\mathrm{rela}}& & \\ A_{\mathrm{ges}}& & \\ {\widetilde{\mathrm{\&eta;}}}_1{A}_{\mathrm{road}}& & \\ {\widetilde{\mathrm{\&eta;}}}_2{A}_{\mathrm{dist}}& & \\ {\widetilde{\mathrm{\&eta;}}}_5{A}_{\mathrm{temp}1}& & \\ {\widetilde{\mathrm{\&eta;}}}_5{A}_{\mathrm{temp}2}& & \\ & {\widetilde{\mathrm{\&eta;}}}_3{B}_{h1}& \\ & {\widetilde{\mathrm{\&eta;}}}_3{B}_{h2}& \\ & {\widetilde{\mathrm{\&eta;}}}_{5}I& \\ & & {\widetilde{\mathrm{\&eta;}}}_4\\ & & {\widetilde{\mathrm{\&eta;}}}_5\end{array}\right)\left(\begin{array}{c}x\\ h\\ H\end{array}\right)\left(\begin{array}{c}{b}_{\mathrm{shp}}\\ b_{\mathrm{rela}}\\ b_{\mathrm{ges}}\\ {\widetilde{\mathrm{\&eta;}}}_1{b}_{\mathrm{road}}\\ {\widetilde{\mathrm{\&eta;}}}_2{b}_{\mathrm{dist}}\\ {\widetilde{\mathrm{\&eta;}}}_5{b}_{\mathrm{temp}1}\\ {\widetilde{\mathrm{\&eta;}}}_5{b}_{\mathrm{temp}2}\\ {\widetilde{\mathrm{\&eta;}}}_3{q}_{h1}\\ {\widetilde{\mathrm{\&eta;}}}_3{q}_{h1}\\ {\widetilde{\mathrm{\&eta;}}}_5\\ {\widetilde{\mathrm{\&eta;}}}_4\\ {\widetilde{\mathrm{\&eta;}}}_5\end{array}\right)\left|\right|}^2.---\left(57\right)$

其中A_(·)和B_(·)分别对应各能量项与二维位置x和建筑高度h高度相关的矩阵，而b_(·)和q_(·)则是与对应各能量项与二维位置x和建筑高度h高度相关的向量，而表示原权重系数η_(·)的平方根。进一步将(47)写作正定二次形式：

$E_{\mathrm{dev}}(x,h,H)={\stackrel{\&OverBar;}{x}}^{T}H\stackrel{\&OverBar;}{x}-2{\stackrel{\&OverBar;}{x}}^T\stackrel{\&OverBar;}{A}d+d^{T}d.---\left(58\right)$

其中 $H={\stackrel{\&OverBar;}{A}}^T\stackrel{\&OverBar;}{A},\stackrel{\&OverBar;}{x}={(x,h,H)}^T.$

最小二乘形式(57)式具有良好的特性。一是矩阵H只与原始形状有关，因此可以预计算矩阵的Cholesky分解，可以很快解出x、h和H。二是与x、h和H相关的系数都分布在对角块上，按各个分块在小规模下分别求解甚至是并行求解，避免了求解整个大规模的线性方程组。

(2)数值求解

1)能量函数的线性化

对于E_occ(x，h，H)，求解能最小化E_occ(x，h，H)消除遮挡。约束可为令各遮挡对的可见性函数 $D_{\{p_r,B_i\}}(x,h,H)\&GreaterEqual;0,$ 有：

$\stackrel{\&OverBar;}{x}=\arg \underset{\stackrel{\&OverBar;}{x}\&Element;P}{\min }{\stackrel{\&OverBar;}{x}}^{T}H\stackrel{\&OverBar;}{x}-2{\stackrel{\&OverBar;}{x}}^T{\stackrel{\&OverBar;}{A}}^{T}d+d^{T}d$

(59)

一般的，道路扩展保持在一定宽度上，因此，令(46)式中的p_r变为常量，使得h_i(p_v-p_r(x))成为线性项。令视点高度H或是建筑位置成为常量，于是可见性函数(48)式变为线性函数

$\begin{array}{c}{D}_{\{p_r,B_i\}}\left(x\right)={l_i}^T(H(p_{B_i}\left(x\right)-p_r)-h_i(p_v-p_r))\\ ={l_i}^T(H({v_{B_i}}^{\&prime;}+t-p_r)-h_i(p_v-p_r))\end{array}.---\left(60\right)$

约束条件变为

${{\mathrm{Hl}}_i}^T{v_{B_i}}^{\&prime;}-h_i{l_i}^T(p_v-p_{r_i})\&GreaterEqual;{{\mathrm{Hl}}_i}^T(p_{r_i}-t),i=\mathrm{1,2}...m.---\left(61\right)$

由此将(59)式根据(57)式重写作

$s.t.G\left(\begin{array}{c}x\\ h\end{array}\right)\&GreaterEqual;c---\left(62\right)$

$A=\left(\begin{array}{c}{A}_{\mathrm{shp}}\\ A_{\mathrm{rela}}\\ A_{\mathrm{ges}}\\ {\widetilde{\mathrm{\&eta;}}}_1{A}_{\mathrm{road}}\\ {\widetilde{\mathrm{\&eta;}}}_2{A}_{\mathrm{dist}}\\ {\widetilde{\mathrm{\&eta;}}}_5{A}_{\mathrm{temp}1}\\ {\widetilde{\mathrm{\&eta;}}}_5{A}_{\mathrm{temp}2}\end{array}\right),b=\left(\begin{array}{c}{b}_{\mathrm{shp}}\\ b_{\mathrm{rela}}\\ b_{\mathrm{ges}}\\ {\widetilde{\mathrm{\&eta;}}}_1{b}_{\mathrm{road}}\\ {\widetilde{\mathrm{\&eta;}}}_2{b}_{\mathrm{dist}}\\ {\widetilde{\mathrm{\&eta;}}}_5{b}_{\mathrm{temp}1}\\ {\widetilde{\mathrm{\&eta;}}}_5{b}_{\mathrm{temp}2}\end{array}\right),B=\left(\begin{array}{c}{\widetilde{\mathrm{\&eta;}}}_3{B}_{h1}\\ {\widetilde{\mathrm{\&eta;}}}_3{B}_{h2}\\ {\widetilde{\mathrm{\&eta;}}}_{5}I\end{array}\right),q=\left(\begin{array}{c}{\widetilde{\mathrm{\&eta;}}}_3{q}_{h1}\\ {\widetilde{\mathrm{\&eta;}}}_3{q}_{h1}\\ {\widetilde{\mathrm{\&eta;}}}_4\end{array}\right).$

(62)式可以用有效集(ActiveSet)法迭代地解出。

2)用Kalman滤波方法实现视点高度的平滑变化

将视点高度H下消除遮挡所需的建筑位移及建筑高度变化能量表示为H的函数，(62)式中G和c都为H的线性函数。(62)式的解为

$\begin{array}{c}\stackrel{\&OverBar;}{x}=H^{-1}({\stackrel{\&OverBar;}{A}}^{T}d+\mathrm{G\&lambda;})\\ \mathrm{\&lambda;}={\left(G{H}^{-1}{G}^T\right)}^{-1}(c-{\mathrm{GH}}^{-1}{\stackrel{\&OverBar;}{A}}^{T}d)\end{array}.---\left(63\right)$

是H的非线性函数。给出一种启发式方法求解视点高度H。设造成遮挡的建筑集合为优化得到包含中建筑位置的向量与高度的向量(59)式可写作：

因为(64)中各建筑的v_i′和h′_i相互独立的，所以可快速求解和视点高度变化，为此需要将变化||v_i′-v_i||和|h′_i-h_i|分别转化为视点高度H的函数。由于

$|{h_i}^{\&prime;}-h_i|=t\left|\frac{\left|\right|P_B-{P_R}^{*}\left|\right|}{\left|\right|P_V-{P_R}^{*}\left|\right|}({H^{\max }}_i-H)\right|=t\left|k_i({H^{\max }}_i-H)\right|.---\left(65\right)$

$||{v_i}^{\&prime;}-v_i||=t|\frac{\left|\right|P_B-{P_R}^{*}\left|\right|}{\left|\right|H\left|\right|}({H^{\max }}_i-H)|=s|l_i({H^{\max }}_i-H)|.---\left(66\right)$

可得到最终视点高度的能量函数，有：

其中l_i＝||P_B-P_R ^*||/H。反复迭代求解(67)式，使得E_H(H)最小化。迭代开始之前先计算与H无关的l_i和k_i，每次迭代更新集合以及参数t_i和s_i，直到H收敛为止。为了满足视点的时间一致性，可使用Kalman滤波方法令视点高度平滑变化。

实施例1：

采用北京市海淀区城区的一块地形与建筑数据(共有1608幢建筑和1096个道路顶点)对本发明的方法进行了实施。为了充分说明本发明的效率，只绘制了未叠加纹理建筑的平顶模型，而对于道路数据则采用线状条带表示。在未打开全景视图模式下，系统可以达到近110fps的帧速率。

在全景视图模式下，本发明的方法依然能达到了交互级速率，三维城市场景能在约为30-40fps的交互级速率下运行。预处理过程主要计算包括城区的空间认知分析、建立稀疏正定矩阵及计算其Cholesky分解、矩阵求逆。建立DT子图及提取格式塔聚类共用时1.86秒。根据场景产生的矩阵大小为分别为2176×2176和608×608，建立两个矩阵并计算其Cholesky分解总共耗时0.37秒。预处理中主要计算开销在于求两个矩阵的逆矩阵，用时12.76秒，但这里仅用了一个线程，而实际上逆矩阵计算可以通过多线程按列并行地计算来加速。在各帧中，主要计算包括可见性光线跟踪时间和迭代求解二次规划时间。基于排序表这种良好的动态加速结构，平均每100次光线跟踪仅用时不到1毫秒。因为预计算了正定矩阵的Cholesky分解及其逆矩阵，平均求解一次二次规划的时间仅为约0.015秒，从而可以很好地保证系统运行时的交互性能。此外对于各帧的视点高度计算，实施中发现视点高度H很快达到收敛。内存消耗约为90MByte，主要用于存储了预计算的矩阵Cholesky分解和逆矩阵。

选定一段道路与视点作为实验对象，对原视图与全景地图进行了比较。道路扩宽大小设定为原道路宽度的2倍，原始视点高度被设定为100m。

图4(a)原始图中拐弯处及拐弯之前的路径都被建筑群遮挡，而且拐弯后的路径大部分被遮挡。图4(b)中，本发明保持视图相似的前提下，有效维持了道路的可见性。图5中道路的形状更为复杂。图5(a)原始图中，感兴趣道路被附近建筑遮挡，而遮挡在图5(b)中沿街的道路被朝视点方向平移，并被有限度地降低高度。图6(a)-(b)中的道路比较弯曲，图6(b)可以在保持大部分建筑高度的前提下消除遮挡。图7(a)-(b)给出了一条曲折的Z字形道路作为感兴趣特征。对于Z字形道路，作为遮挡物的建筑在平移过程中可能会被临近道路阻挡，在图7(b)中，大多建筑能被平移到不造成遮挡的位置，而城市的几何形态保持不变。表1列出了全景图中视点高度、最大建筑位移和最大建筑降低高度。

表1.全景图中的视点高度、最大位移与最大降低高度(原视点高度为100米，单位：米)

图	全景视图视点高度	最大位移	最大降低高度
				图4	138.7	25.4	12.0
图5	286.2	60.4	12.9
				图6	230.7	53.5	23.1
图7	196.5	68.2	15.2

可以发现，本发明不仅较好保持视点高度，建筑位置与建筑高度也服从了空间认知规律，并且有效消除了道路的遮挡，所产生的视图达到了交互级的帧速率。

Claims (1)

1.一种交互式三维城市全景地图的自动生成方法，其特征在于，包括如下步骤：

步骤一：基于ε距离邻域和β-骨架规则识别街区内建筑空间分布模式

所要构建DT图的结点为建筑足迹(footprint)多边形的重心，结点之间的距离为建筑足迹多边形间的最小距离；先构造结点的Delaunay三角网，从Delaunay三角网内筛选连接边，建立结点的最小生成树MST；下一步，基于ε距离邻域规则，在Delaunay三角网内选择距离小于阈值ε，且满足β＜1时，把β-骨架规则的非MST边加入子图中，得到满足ε距离邻域条件的DT子图；对于DT子图中度为1的结点，在Delaunay三角网中寻找其满足β＜1时，β-骨架规则最短的非MST边，并将该边加入子图中；

步骤二：提取城市建筑格式塔聚类

设定合理的阈值以找出潜在的格式塔聚类，进而用经训练的支持向量机分类器加以判别；整个建筑的格式塔聚类提取分为2步：(1)确定潜在的建筑格式塔聚类，(2)提取具有格式塔特征的建筑；

步骤三：基于光线跟踪的城市感兴趣特征可见性计算

利用光线跟踪实现城区道路可见性判断，将一条光线和与其相交的遮挡物作为一个遮挡对；每个遮挡对记录有光线、遮挡物及光线击中点参照于遮挡物重心的相对坐标；对于光线与建筑的相交测试，先判断光线在二维上是否相交，快速排除不相交的建筑，再在三维上判断光线是否与建筑的各个面相交；如有交点，则建筑为遮挡物，取距离参考点最近的交点为击中点，否则建筑不造成遮挡；采用动态的KD树结构加快二维平面上光线跟踪；

步骤四：建立城市全景图的约束条件和能量项

(1)建立城市全景图的能量函数偏离项

为了保证格式塔聚类的共同命运特征，建立街区内格式塔聚类位移约束；对道路形状、道路与建筑冲突进行约束，以保证街区与道路的相对位置，此外，约束建筑高度与视点高度，最后针对漫游时出现的视图抖动与闪烁现象，建立时空一致性约束；

(2)建立城市全景图能量函数的遮挡项

定义x为包含建筑重心与道路条带顶点的二维位置的向量，h为建筑高度向量；对于一个遮挡对{p_r，B_i}，参考点p_r的可见性表示为x的函数；设造成遮挡的建筑B_i高度为h_i，足迹重心为v′_Bi，与其相邻的道路线段为v′_i0v′_(i+1)0，则视线在建筑上的击中点为t_{pr，Bi}为p_r、之间的向量，且p_r由其相邻道路顶点表示p_r＝(1-s)v_i0′+sv_(i+1)0′，0＜s＜1，有可见性函数的二次形式：

$D_{\{p_r,B_i\}}(x,H)\&ap;{I_{\mathrm{LOS}}}^T\left(x\right)[H(p_{B_i}\left(x\right)-p_r\left(x\right))-h_i(p_v-p_r\left(x\right))]---\left(1\right)$

式中p_v为视点二维位置，H为其高度；

考虑所有遮挡对集合O，将(1)式离散化表示为x，h和H的函数(2)，使感兴趣道路可见，

由此能量函数的遮挡项表示为：

$E_{\mathrm{occ}}(x,h,H)=\left\{\begin{array}{cc}-D(x,h,H),& D(x,h,H)<0\\ 0,& D(x,h,H)\&GreaterEqual;0\end{array}\right.---\left(3\right)$

步骤五：城市全景图能量函数的最优化求解

(1)计算能量函数的最小二乘形式

综合各约束，得到最终的能量函数偏离项

E_dev(x，h，H)＝(E_shp(x)+E_rela(x)+E_ges(x))+η₁E_road(x)+η₂E_dist(x)(4)+η₃(E_h1(h)+E_h2(h))+η₄E_H(H)+η₅E_temp(x，h，H)

上述各项都具有正定二次形式，将各项最小二乘形式有：

$E_{\mathrm{dev}}(x,h,H)=\left|\right|\left(\begin{array}{ccc}{A}_{\mathrm{shp}}& & \\ A_{\mathrm{rela}}& & \\ A_{\mathrm{ges}}& & \\ {\widetilde{\mathrm{\&eta;}}}_1{A}_{\mathrm{road}}& & \\ {\widetilde{\mathrm{\&eta;}}}_2{A}_{\mathrm{dist}}& & \\ {\widetilde{\mathrm{\&eta;}}}_5{A}_{\mathrm{temp}1}& & \\ {\widetilde{\mathrm{\&eta;}}}_5{A}_{\mathrm{temp}2}& & \\ & {\widetilde{\mathrm{\&eta;}}}_3{B}_{h1}& \\ & {\widetilde{\mathrm{\&eta;}}}_3{B}_{h2}& \\ & {\widetilde{\mathrm{\&eta;}}}_{5}I& \\ & & {\widetilde{\mathrm{\&eta;}}}_4\\ & & {\widetilde{\mathrm{\&eta;}}}_5\end{array}\right)\left(\begin{array}{c}x\\ h\\ H\end{array}\right)-\left(\begin{array}{c}{b}_{\mathrm{shp}}\\ b_{\mathrm{rela}}\\ b_{\mathrm{ges}}\\ {\widetilde{\mathrm{\&eta;}}}_1{b}_{\mathrm{road}}\\ {\widetilde{\mathrm{\&eta;}}}_2{b}_{\mathrm{dist}}\\ {\widetilde{\mathrm{\&eta;}}}_5{b}_{\mathrm{temp}1}\\ {\widetilde{\mathrm{\&eta;}}}_5{b}_{\mathrm{temp}2}\\ {\widetilde{\mathrm{\&eta;}}}_3{q}_{h1}\\ {\widetilde{\mathrm{\&eta;}}}_3{q}_{h1}\\ {\widetilde{\mathrm{\&eta;}}}_5\\ {\widetilde{\mathrm{\&eta;}}}_4\\ {\widetilde{\mathrm{\&eta;}}}_5\end{array}\right){\left|\right|}^2--\left(5\right)$

其中A_(·)和B_(·)分别对应各能量项与二维位置x和建筑高度h高度相关的矩阵，而b_(·)和q_(·)则是与对应各能量项与二维位置x和建筑高度h高度相关的向量，而表示原权重系数η_(·)的平方根；

(2)数值求解

1)能量函数的线性化

对于E_occ(x，h，H)，求解能最小化E_occ(x，h，H)消除遮挡，约束为令各遮挡对的可见性函数 $D_{\{p_r,B_i\}}(x,h,H)\&GreaterEqual;0,$ 有：

$\stackrel{\&OverBar;}{x}=\arg \underset{\stackrel{\&OverBar;}{x}\&Element;P}{\min }{\stackrel{\&OverBar;}{x}}^{T}H\stackrel{\&OverBar;}{x}-2{\stackrel{\&OverBar;}{x}}^T{\stackrel{\&OverBar;}{A}}^{T}d+d^{T}d---\left(6\right)$

(6)式线性化后

$s.t.G\left(\begin{array}{c}x\\ h\end{array}\right)\&GreaterEqual;c$

$A=\left(\begin{array}{c}{A}_{\mathrm{shp}}\\ A_{\mathrm{rela}}\\ A_{\mathrm{ges}}\\ {\widetilde{\mathrm{\&eta;}}}_1{A}_{\mathrm{road}}\\ {\widetilde{\mathrm{\&eta;}}}_2{A}_{\mathrm{dist}}\\ {\widetilde{\mathrm{\&eta;}}}_5{A}_{\mathrm{temp}1}\\ {\widetilde{\mathrm{\&eta;}}}_5{A}_{\mathrm{temp}2}\end{array}\right),b=\left(\begin{array}{c}{b}_{\mathrm{shp}}\\ b_{\mathrm{rela}}\\ b_{\mathrm{ges}}\\ {\widetilde{\mathrm{\&eta;}}}_1{b}_{\mathrm{road}}\\ {\widetilde{\mathrm{\&eta;}}}_2{b}_{\mathrm{dist}}\\ {\widetilde{\mathrm{\&eta;}}}_5{b}_{\mathrm{temp}1}\\ {\widetilde{\mathrm{\&eta;}}}_5{b}_{\mathrm{temp}2}\end{array}\right),B=\left(\begin{array}{c}{\widetilde{\mathrm{\&eta;}}}_3{B}_{h1}\\ {\widetilde{\mathrm{\&eta;}}}_3{B}_{h2}\\ {\widetilde{\mathrm{\&eta;}}}_{5}I\end{array}\right),q=\left(\begin{array}{c}{\widetilde{\mathrm{\&eta;}}}_3{q}_{h1}\\ {\widetilde{\mathrm{\&eta;}}}_3{q}_{h1}\\ {\widetilde{\mathrm{\&eta;}}}_4\end{array}\right).---\left(7\right)$

2)用Kalman滤波方法实现视点高度的平滑变化

为了满足视点的时间一致性，使用Kalman滤波方法令视点高度平滑变化。