经过这几天的测试,完成无人机航拍图像几何校正,也叫做航片的正射影像处理。 主要步骤如下:
现在对每一个步骤分析,谈论为什么要这样做。 获取航片信息
获取航片信息是为了计算影像的分辨率,相关操作可以查看这里。 根据偏航角旋转图片 这步骤是必要的。详细讲解如何实现这一步骤。
扩充图像的意思是,用0值在图像的边缘填充,这一步骤的目的是,为了 旋转后的图像还在指定尺寸范围内。 扩充图像效果如下:
下一步是对图像进行旋转,假设偏航角为-45°,则此时要对图像进行逆时针旋转45°。 这涉及到另一个问题,偏航角的定义。示意图如下:
简单来说,航迹线偏在航线右边,偏航角为正,航迹线偏在航线左边,偏航角为负。请看下面的灵魂画图。
这里假设偏航角是-45度,旋转中心为图像的中心,旋转后的图像如下:
为什么要先扩充图像,再旋转?因为这里调用opencv的函数,如果直接对原始图像进行旋转,会丢失部分信息。 直接对原始图像进行旋转,效果如下:
如上图所示,部分信息被省略了,所以要先用0值扩充图像边缘,使图像的尺寸增加一倍,然后再对扩充后的图像进行旋转。 裁剪图像 这一步骤的目的是,为了减少旋转图像的无效像素,缩小图像的尺寸的同时,仅保留有用信息。 裁剪后的旋转图像的效果如下:
计算左上角坐标 现在得到的结果,是不包含经纬度信息的。我们在第一步的时候,获取了图像的中心的的经纬度、高度、镜头参数。所以我们可以计算此时无人机拍摄照片的空间分辨率。
空间分辨率 = (像元尺寸 x 高度)/ 焦距
以上的单位都要换算为米。
像元尺寸要根据你的无人机硬件参数决定。现在我们的无人机镜头的像元尺寸为4.4微米。 高度、焦距自动从影像文件中获取。
此时的空间分辨率为米,我们需要进行单位换算,把米换算为度。 假设地球的半径为6371004米,则有以下的换算公式:
x_res = x_res_meter / (2 * math.pi * 6371004) * 360
左上角经纬度换算公式如下: 左上角经度 = 中心点经度 - 裁剪后的旋转图像的宽度 x 0.5 x 经度分辨率
左上角纬度 = 中心点纬度 + 裁剪后的旋转图像的宽度 x 0.5 x 纬度分辨率
仿射六参数拟合
通过以上的步骤,我们获得了 1.左上角经纬度 2.空间分辨率 所以我们可以得到该影像的仿射六参数。其格式如下: [左上角经度,经度分辨率,0,左上角纬度,0,纬度分辨率 ] 由于我们在北半球,所以,上面公式中的经纬度分辨率换算公式如下: 经度分辨率 = 空间分辨率 纬度分辨率 = -1 x 空间分辨率 合成为TIF格式 这里跳用GDAL库,轻松地把裁剪后的旋转图像与仿射六参数结合,生成TIF格式的影像。 最终结果展示 google earth上的测试区范围如下:
把结果拖到google earth上,如下图所示:
讨论
以上可以认为是航片的粗几何校正。存在有一定的偏差,在2米左右。
这是因为我们在计算空间分辨率的时候,是假设地球半径。实际上当地的地球半径不是我们假设的通用地球半径。优化思路可以根据当地的经纬度,进行UTM坐标转换为WGS84坐标。
如果要精细化校正,则再需要进行特征点匹配,进行航片与底图的配准、航片与航片直接的配置。 现在的商业软件基本都集成了航片的正射校正,相关的开源资料却很少,以上的步骤是慢慢思考,总结出来的。
所以,这篇文章可能在存在一些错误,欢迎在评论区留言指出,一起学习!
代码已开源,地址在: https://github.com/ytkz11/AerialPictureCorrection
欢迎给我一个免费的star,这次代码和样本数据就不放在百度云了,直接放在github。
样本数据是一个小伙伴提供的,现在只分享一景,作为测试数据。 之后会对代码进行讲解。
明天先把代码整合,使其可以批量对航片处理,最后打包为exe,会用云盘的方式分享出来。 吭吭哧哧花费了好几个小时才完成这篇干货技术文章,你们随一个赞是不过分的吧?!
往期文章回顾