一种基于MVS稠密点云的场景正射影像生成方法

一种基于mvs稠密点云的场景正射影像生成方法
技术领域
1.本发明属于地景建模与测绘技术领域，涉及一种利用多视角立体三维重建(mvs)算法得到的地面场景稠密点云生成正射影像的方法。


背景技术：

2.随着计算机视觉技术的不断的发展，对地面场景进行建模并且生成目标场景的正射影像也处于越来越重要的地位。对地面场景建模常见的一种方法是利用倾斜影像采集地面场景的图像数据，再通过图像进一步完成地面场景三维模型的重建。具体而言，首先利用无人机倾斜摄影获得建模区域的图像，其次通过运动推断结构(structrue-from-motion，sfm)算法得到相机位置和姿势，然后利用多视图立体三维重建算法(mulit-view steros，mvs)得到地面场景的三维模型。其中sfm算法是可以从多幅二维图像中推断出相机的参数和场景三维结构的技术，其中相机的参数包括相机的内参和外参；而场景的三维结构是指图像中特征点的三维信息，称为稀疏点云。mvs三维重建是利用一组标定的二维多视角图像重建三维几何形状的算法；实现方法是利用sfm获取二维图像的相机参数和稀疏点云，重建出目标区域的详细信息，常见的输出结果是目标模型的稠密点云。
3.本发明在对目标地面场景建模的过程中，使用已生成的稠密点云生成目标地面场景的正射影像。
4.正射影像是具有正射投影性质的图像，其具有现势性和完整性。正射影像的快速准确的获取有着至关重要的作用，在实际应用中正射影像可作为背景控制信息，可评价其他数据的精度，可从中提取自然和人文信息，可用做地图的更新，也可用于数字城市建设和城市规划设计。对于地面场景建模而言，在建模的过程中，获得正射影像可以用来进一步完善建模结果，对目标场景进行分类，提高效率，优化建模的效果。
5.传统的正射影像生成方法是利用中心投影影像通过数字纠正的方法生成。但是在获取地面场景稠密点云的前提下，利用本发明所提出的方法可以直观准确的生成正射影像，可以避免在由中心投影影像制作正射影像过程中，逐像元进行辐射纠正、微分纠正和镶嵌所产生的运算量和误差。


技术实现要素：

6.本发明针对目标区域地面场景正射影像的生成问题，提供一种基于mvs稠密点云生成正射影像的方法。以得到的正射影像为基础，可以进一步对目标区域进行分析，也可以完善地面场景的重建结果。
7.本发明提供的所述一种基于mvs稠密点云的场景正射影像生成方法，包括以下步骤：
8.s1：利用倾斜摄影的方法获得目标地面场景区域的航空影像，利用sfm算法和mvs算法获取目标地面场景区域稠密点云，然后利用滤波算法去除稠密点云中的离群点；
9.s2：利用地面标记点、sfm算法所得相机位置计算出稠密点云中对应的水平面表达
式，对稠密点云进行旋转得到投影平面；
10.s3：计算正射影像的图像分辨率，利用稠密点云中三维点的颜色初步生成正射影像；
11.s4：利用相邻像素的颜色rgb信息处理空像素，优化正射影像最终效果。
12.进一步地，所述步骤s1中生成目标地面场景区域稠密点云的算法流程为：
13.s1-1：在无人机采集目标区域地景航拍图像的过程前，地面设置t个标记点；
14.s1-2：操作无人机采集的目标建模区域的图像，并获飞行过程中无人机地理位置信息；
15.s1-3：经过sfm获取相机外参与图像中特征点的三维信息，其中外参是指相机的位置与姿势；
16.s1-4：经过mvs重建获取的目标区域稠密点云，稠密点云中包含地面标记点对应的三维点集；稠密点云要求基本可以覆盖目标建模场景，没有大面积的缺失信息；
17.s1-5：利用pcl(point cloud library)中提供的高斯滤波、条件滤波、基于随机采样一致性滤波等滤波方法去除稠密点云中的离群点。
18.至此，已经获得了目标区域地面场景的稠密点云，稠密点云基本覆盖目标区域，没有大面积的缺失信息，点云基本没有离群点。
19.进一步地，所属步骤s2中利用地面标记点、相机位置计算出稠密点云中对应的水平面，对目标场景稠密点云进行旋转得到投影平面的算法流程为：
20.s2-1：计算出一个参考向量用于判断相对水平面法向量的方向；
21.s2-2：当在s1-1步骤中设置地面标记点时，利用地面标记点信息求解目标场景稠密点云中的相对水平面；
22.s2-3：若没有设置地面标记点，利用sfm得到的相机位置求解目标场景稠密点云中的相对水平面；
23.s2-4：利用相对水平面表达式计算出旋转矩阵r，通过旋转矩阵r将稠密进行旋转，使得场景中的水平面方向向量为(0,0,1)
t
，将旋转后的水平面视为投影平面。
24.至此，已经计算出稠密点云中相对水平面的方程表达式，并将其作为投影平面，在后续步骤中得到正射影像。
25.进一步地，所述步骤s3计算正射影像的图像分辨率，利用稠密点云中三维点的颜色初步生成正射影像的算法流程为：
26.s3-1：从完成旋转的稠密点云中，选择生成正射影像的区域，记录其包围盒；
27.s3-2：计算正射影像图像分辨率m
×
n，其中，若设置地面标记点，利用地面标记点和设置的地面分辨率，计算出正射影像图像分辨率；若没有设置地面标记点，则利用生成正射影像的区域范围和点云密度，计算出正射影像图像分辨率；
28.s3-3：计算出点云中的三维点对应的像素坐标(i,j)，利用点云中三维点的颜色rgb值，初步确认每个像素的rgb颜色值。如果像素对应三维点数目为0，则记录空像素(i,j)，在后续的步骤中进行进一步处理空像素；如果对应三维点数目为1，像素(i,j)的rgb颜色值等于三维点rgb颜色值；如果对应三维点数目大于1，像素(i,j)的rgb颜色值等于多个三维点中z值最大的三维点的rgb颜色值。由于本发明要求稠密点云可以覆盖目标场景，没有大面积的缺失。所以空像素的数量较少，可以通过修复，进一步优化生成正射投影图的显
示结果。
29.进一步地，所述步骤s4处理空像素，优化正射影像最终显示效果的算法流程为：
30.s4-1：若(i,j)为空像素，8-邻域像素全部不为空，则利用8-邻域像素的rgb颜色值计算(i,j)的rgb颜色值；
31.s4-2：若空像素(i,j)的8-邻域像素内的像素全部为空，则从相邻像素不全为空的像素开始计算；
32.s4-3：若空像素(i,j)的8-邻域内的像素部分为空，则只计算不为空的像素的颜色rgb的均值；
33.s4-4：根据每个空像素(i,j)的8-邻域内的像素的情况，计算出所有像素的rgb值。
34.本发明实施提供的一种利用目标场景中的稠密点云生成正射影像的方法。利用此正射影像对场景模型进行完善，也可以对目标场景进行分析。与传统生成正射影像的方法相比，在获得目标区域场景稠密点云的前提下，本发明的优点在于：
35.(1)可以通过稠密点云方便快捷的获得目标区域的正射影像，计算方式直接快捷。
36.(2)所生成的正射影像是直接从三维信息中生成的，所以避免传统方法中从倾斜影像中生成正射影像时出现的遮蔽现象，进而避免影像校准的过程，提高正射影像的准确率。
附图说明
37.为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，通过参考附图会更加清楚的理解本发明的特征和优点，附图是示意性的而不应理解为对本发明进行任何限制，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，可以根据这些附图获得其他的附图。其中：
38.图1为本发明实施例提供的基于mvs算法重建地面场景稠密点云的正射影像方法流程图；
39.图2为获得目标区域地面场景稠密点云流程图
40.图3为目标区域地面场景稠密点云局部示意图；
41.图4为旋转目标地面场景稠密点云得到投影平面流程图；
42.图5为利用目标地面场景稠密点云生成正射影像流程图；
43.图6为利用稠密点云中三维点的颜色初步生成正射影像局部示意图；
44.图7为经过优化消除空像素局部示意图
45.图8为目标区域地面场景正射影像示意图
具体实施方式
46.为使本发明的技术方案更加清楚，下面结合附图对本发明做进一步阐述。整体流程图如图1所示。
47.具体地，一种基于mvs稠密点云的场景正射影像生成方法，具体包括以下步骤：
48.s1：利用倾斜摄影的方法获得目标地面场景区域的航空影像，利用sfm算法和mvs算法重建出所述目标区域的稠密点云，并通过滤波算法去除所述稠密点云中的离群点。
49.s2：利用地面标记点和sfm算法所得相机位置计算出所述稠密点云中对应的水平
面表达式，旋转目标场景稠密点云得到投影平面。由于经过mvs重建出的地面场景稠密点云的水平面并不与xyz直角坐标系中的某个坐标系平面平行。所以在此步骤中，旋转稠密点云使得稠密点云中的相对水平面与xoy平面平行。
50.s3：计算正射影像图像分辨率，利用稠密点云中三维点的颜色信息初步生成正射影像。
51.s4：利用相邻像素的颜色rgb值消除s3步骤生成的正射影像中存在的空像素。由于本发明要求稠密点云可以覆盖目标场景，没有大面积的缺失。所以空像素的数量较少，可以进行修复，进一步优化生成正射投影图的显示结果。
52.如图2流程图所示，步骤s1中获得目标区域地面场景稠密点云具体算法流程为：
53.s1-1：在无人机采集目标区域地景航拍图像的过程前，地面设置t个标记点，其中t＞2；
54.s1-2：对搭载倾斜相机的无人机进行巡航任务规划和手动操作，采集建模出目标区域的航拍图片，在此过程中记录无人机地理位置信息(pos点)，pos点主要包括gps数据和imu数据，即倾斜摄影测量中的外方位元素，其中gps数据包括相机拍摄时经度、维度和高程数据；对于在地面设置的标记点，航拍图像中包括标记点信息；
55.s1-3：利用sfm得到图像对应的相机在稠密点云三维坐标系中的位置和姿势，并得到图像中特征点的三维信息；开源sfm算法包括：opensfm，openmvg，colmap等；
56.s1-4：利用sfm得到的相机位置姿势和特征点三维信息，经过mvs算法，重建得到目标区域的稠密点云。开源mvs算法包括：pmvs，openmvs，colmap等；此步骤中，会生成一组包含与稠密点云中的三维点集p
mark
，p
mark
与航拍图像中的地面标记点p
real
一一对应；为了最终的正射影像的完整性，本发明要求稠密点云可以覆盖目标区域的场景，没有大面积的缺失信息；
57.s1-5：利用pcl库提供的基于随机采样一致性滤波的方法，删除稠密点云中的离群点；
58.最终目标区域建模稠密点云局部如图3所示，至此，已经获得了目标区域地面场景的稠密点云，稠密点云基本覆盖目标区域，没有大面积的缺失信息，点云基本没有离群点。
59.稠密点云的相对水平面并不与xyz直角坐标系中的某个坐标系平面平行，在步骤s2中旋转稠密点云，将模型中的相对水平面旋转到与xoy平面平行，将相对水平面作为投影平面。如图4流程图所示，步骤s2具体算法流程为：
60.s2-1：在稠密点云中，随机选取一个三维点，记为p
t
(x
t
,y
t
,z
t
)
t
；随机选择相机位置记为pf(xf,yf,zf)
t
，利用p
t
和pf计算出一个参考向量如公式(1所示，在后续步骤中用于判断模型中相对水平面法向量的方向；
[0061][0062]
s2-2：当在s1-1步骤中地面设置标记点时，求解目标场景稠密点云中的相对水平面，水平面方程表达式记为a1x+b1y+z+c1＝0；真实场景中的地面标记为p
real
＝{p1,p2,
…
,pn}，其中pi＝(xi,yi,zi)
t
，而在稠密点云三维坐标中与p
real
一一对应的三维点集为p
mark
＝{p1,p2,
…
,pn}，其中pi＝(xi,yi,zi)
t
；
[0063]
计算p
real
中z值最大值z
max
与最小值z
min
，若(z
max-z
min
)≤0.05m，则可以判断目标场
景内地面基本与水平面平行，即与p
real
一一对应的点集p
mark
拟合出的平面为模型场景中的相对水平面；采用最小二乘法利用三维点集合p
mark
＝{p1,p2,
…
,pn}拟合点云中的水平面；水平面方程表达式记为a1x+b1y+z+c1＝0，最小二乘法矩阵表示如公式(2)所示：
[0064]
ax＝b
ꢀꢀ
(2)
[0065]
其中
[0066]
其中(x1,y1,z1)
t
,(x2,y2,z2)
t
,
…
,(xn,yn,zn)
t
表示p
mark
中三维点坐标，使用正规方程x＝(a
t
a)-1at
b求得x中a1,b1,c1的值，此时相对水平面法向量记为若向量(a1,b
1,1t
与参考向量之间的夹角为锐角，则平面法向量记为否则平面法向量记为
[0067]
若(z
max-z
min
)＞0.05m，表示目标场景内地面与水平面不平行，将地面标记点p
real
中z值代表的高度设置为同一常数，即可表示为真实场景中同一水平面的一组点，记为p
horizontal_real
＝{p
h,1
,p
h,2
,
…
,p
h,n
}，其中所有的z值均设置为t，p
horizontal_real
集合中的p
h,i
与p
real
中的pi有相同x值和y值，即若pi＝(xi,yi,zi)
t
，则有p
h,i
＝(xi,yi,t)
t
；
[0068]
任选两地面标记点pj和pk之间的实际距离，记为ln(单位：米m)；在稠密点云中计算地面标记点pj和pk在模型中的距离，记为lm，则模型中的单位长度对应的真实场景中的实际距离记为l，计算方法如公式(3)所示：
[0069][0070]
由于在实际运算过程中pj和pk可能会存在较大误差，影响最终结果，因此l的计算如公式(4)所示：
[0071][0072]
其中pi表示p
mark
中第i个点，pi表示p
real
中第i个点，||p
i+1-pi||计算pi与p
i+1
之间的距离；
[0073]
点集合p
horizontal_real
中的一点p
h,i
在稠密点云中的坐标p
h,i
＝(x
h,i
,y
h,i
,z
h,i
)
t
可以通过公式(5)计算：
[0074][0075]
其中zi表示与p
h,i
有相同x值和y值的地面标记点pi的z值，(xi,yi,zi)
t
表示地面标记点pi对应的稠密点云中的三维点pi的三维坐标。(a,b,c)
t
表示目标场景点云中相对水平面的单位法向量，即a2+b2+c2＝1；相对水平面的法向量可等价表示为其中
[0076]
若向量(a1,b1,1)
t
与参考向量之间的夹角为锐角，则平面法向量记为之间的夹角为锐角，则平面法向量记为否则平面法向量记为
[0077]
s2-3：若在采集图像的过程中，若地面没有设置标记点，则用如下方法计算模型中相对水平面的表达式a1x+b1y+z+c1＝0；在无人机搭载倾斜相机采集数据时，获取的无人机影像会携带配套的pos点数据。使用pos点数据筛选出无人机保持在同一高度h进行巡航时获取的图像，结合这些图像经过sfm算法得到的相机位置，利用最小二乘法可以拟合出稠密点云场景坐标中高度为h的平面方程，此平面方程与水平面平行；若向量(a1,b1,1
t
与参考向量之间的夹角为锐角，则平面法向量记为否则平面法向量记为否则平面法向量记为
[0078]
s2-4：使用rodrigues公式，计算出向量与向量(0,0,1)的旋转矩阵r，当时，旋转矩阵r如公式(6)所示：
[0079][0080]
其中，表示向量(a1,b1,1
t
与向量(0,0,1)
t
夹角的余弦值，此时平面法向量可以通过旋转矩阵r进行旋转，旋转后结果为
[0081]
通过旋转矩阵r将稠密点云进行旋转，使得场景中的水平面平行于xoy平面，并将旋转后的水平面视为投影平面。对于稠密模型中的任意三维点p
before
，可得旋转后的点p
after
：
[0082]
p
after
＝rp
before
ꢀꢀ
(7)
[0083]
至此，已经计算出稠密点云中相对水平面的方程表达式，并将稠密点云进行旋转，水平面作为投影平面在后续的步骤中生成正射影像。
[0084]
如图5流程图所示，步骤s3和s4利用旋转后的目标区域地面场景稠密点云生成正射影像具体算法流程为：
[0085]
s3-1：从完成旋转的稠密点云中选择生成正射影像的区域。该区域用包围盒进行表示，x方向上最小值x
min
和最大值x
max
，y方向最小值y
min
和最大值y
max
，z方向上的最小值z
min
和最大值z
max
；
[0086]
s3-2：计算正射影像图像分辨率。正射影像的图像分辨率记为m
×
n，其中m和n分别
1,j-1),(i-1,j+1),(i+1,j-1),(i+1,j+1)如果不全为空，则可以根据8-邻域内不为空的像素的rgb值计算出空像素(i,j)的rgb值。
[0102]
s4-1：若(i,j)为空像素，8-邻域像素全部不为空，则计算(i,j)的rgb颜色值如公式(15)、(16)、(17)所示：
[0103][0104][0105][0106]
s4-2：若空像素(i,j)的8-邻域像素内的像素全部为空，则从相邻像素不全为空的像素开始计算。例如，空像素(i,j)的8邻域中，只有(i-1,j)、(i-1,j-1)、(i,j-1)、(i+1,j-1)四个像素有不为空，其他的四个像素都为空时，计算(i,j)的rgb颜色值如公式(18)、(19)、(20)所示：
[0107][0108][0109][0110]
s4-3：若空像素(i,j)的8-邻域内的像素部分为空，则只计算不为空的像素的颜色rgb的均值；
[0111]
s4-4：根据每个空像素(i,j)的8-邻域内的像素的情况，计算出所有像素的rgb值。
[0112]
初步生成正射影像没有进行优化的正射影像效果局部示意图如图6所示，图中白色像素为空像素。经过优化的正射影像局部示意图如图7所示。
[0113]
最终生成的目标区域的正射影像如图8所示。