经过这几天的测试，完成无人机航拍图像几何校正，也叫做航片的正射影像处理。 主要步骤如下：
现在对每一个步骤分析，谈论为什么要这样做。 获取航片信息
获取航片信息是为了计算影像的分辨率，相关操作可以查看这里。 根据偏航角旋转图片 这步骤是必要的。详细讲解如何实现这一步骤。
扩充图像的意思是，用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，会用云盘的方式分享出来。 吭吭哧哧花费了好几个小时才完成这篇干货技术文章，你们随一个赞是不过分的吧？！
往期文章回顾
1334
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
