学生:邬代杰

院系:遥感信息工程学院

学号:2019302130329

班级:19D9

目录

• 一. 基本原理
  • 1. 寻找控制点坐标
  • 2. 对控制点进行编号
  • 3. 计算相机外部参数
  • 4. 计算控制点三维坐标
• 二. 实现过程
  • 1. 寻找控制点
  • 2. 寻找圆盘中心点
  • 3. 对控制点分组
  • 4. 对控制点编号
  • 5. 优化外部参数
• 三. 操作说明

一. 基本原理

1. 寻找控制点坐标

控制点就是圆盘上黑色区域的中心坐标，体现在图像上就是黑色像素集合区域的中心，那我们首先要做的就是提取出黑色像素。由于原图像为RGB彩色图像，佛像与背景、控制点黑色区域与圆盘其他部分对比度较高，故考虑将图像二值化后再进行操作，具体操作应用MATLAB自带函数就可以解决，不再赘述。

2. 对控制点进行编号

提取出控制点坐标后，接下来的任务就是利用已知的控制点编号规则在图像上对控制点进行编号。

注意到控制点是在一个圆盘之上的，即是说：
任意过同一编号组控制点的直线一定都相交于圆盘圆心且过圆心的直线至多与4个控制点中心相交

又有：投影（相机摄像）不会改变物体的相对位置，并且图片已经消除了畸变，不会出现同一组控制点不在一条直线上的情况。

基于以上事实，我们可以得到将控制点编号的方法：

1. 选取两组控制点中各两点，列直线方程；

2. 计算两直线相交点作为圆盘中心点；

3. 逐点进行遍历并分组：首先任意选取一点作为起始点，然后应用直线拟合对其他点进行判定（即是否为同一组控制点）并记录同一组的控制点的像素坐标，然后不断选取未分组的点进行分组操作，直到全部点分组完毕；

4. 将分组后的控制点从外到内进行排序：在分好组的控制点间计算其与中心点之间的棋盘距离，采用冒泡排序以便于同时更新控制点坐标；

5. 编号：根据控制点区域大小赋给控制点不同的值（0或1），4个点得到一个二进制数（不满4个点的控制点组舍弃），依此序列计算组号。以0100为例，组号为：\(0\times 2^3+1\times2^2+0\times2^1+0\times2^0=4\). 若组号为\(n\)，从外到内的编号依次为：

\[N=4(n-1)+i-1\quad(i是从外到内点的位置) \]

至此，编号结束。

3. 计算相机外部参数

采取后方交会计算左右相片各自的外部参数。

后方交会原理在给的参考资料中已经很详细了，不再赘述。

4. 计算控制点三维坐标

原理：前方交会

已知：

\[\begin{align} 相机内参：&x_0,y_0,f,width,height\\ 相片外参：&X_{S1},Y_{S1},Z_{S1},\varphi_1,\omega_1,\kappa_1\\ &X_{S2},Y_{S2},Z_{S2},\varphi_2,\omega_2,\kappa_2\\ 同名像点坐标：&(x_1,y_1),(x_2,y_2) \end{align} \]

待求：

\[控制点三维坐标：X_p,Y_p,Z_p \]

过程：

1. 计算左右片在地辅坐标系中旋转矩阵的方向余弦：该步骤与后方交会时计算旋转矩阵一致；

\[R= \begin{bmatrix} a_1 \quad a_2\quad a_3\\ b_1 \quad b_2\quad b_3\\ c_1 \quad c_2\quad c_3 \end{bmatrix} \\ R'=\begin{bmatrix} a_1' \quad a_2'\quad a_3'\\ b_1' \quad b_2'\quad b_3'\\ c_1' \quad c_2'\quad c_3' \end{bmatrix} \]

2. 计算基线分量

\[ B _ { X } = X _ { s2 } - X _ { s1 } \\ B _ { Y } = Y _ { s 2 } - Y _ { s1 } \\ B _ { Z } = Z _ { s 2 } - Z _ { s 1 } \]

3. 计算像点的像空间辅助坐标

\[ \begin{bmatrix} X\\Y\\Z \end{bmatrix} =R\begin{bmatrix} x_1\\y_1\\-f \end{bmatrix}, \begin{bmatrix} X'\\Y'\\Z' \end{bmatrix} =R\begin{bmatrix} x_2\\y_2\\-f \end{bmatrix} \]

4. 计算投影系数

\[ N=\frac{B_X Z'-B_Z X'}{XZ'-ZX'}\\ N'=\frac{B_X Z-B_Z X}{XZ'-ZX'} \]

5. 计算地面点的左像辅系坐标

\[ \begin{align} \Delta X&=NX\\ \Delta Y&=(NY+N'Y'+B_Y)/2\\ \Delta Z&=NZ \end{align} \]

6. 计算地面点的地面坐标

\[ \begin{align} X_p&=X_{s1}+\Delta X\\ Y_p&=Y_{s1}+\Delta Y\\ Z_p&=Z_{s1}+\Delta Z\\ \end{align} \]

二. 实现过程

以下均以左边图像为例

1. 寻找控制点

main.m第52~74行:

RGB转灰度图-->灰度图转二值图-->对二值图取质心及面积-->储存于结构体ldftInfo.coodInfo中

2. 寻找圆盘中心点

main.m第76~85行

选定第1,13,28,40点计算中心点坐标-->计算交点(中心点坐标)-->存于ldftInfo.centreInfo

3. 对控制点分组

main.m第87~119行

flag数组用于标记是否分组-->将选择出来的未分组的控制点与中心点进行一次幂拟合,得到系数-->遍历所有控制点,寻找在一条直线上的点-->对有的不足4个元素的数组进行补充-->重复迭代直至遍历所有控制点

4. 对控制点编号

距离中心点的距离从小到大进行排序-->求二进制编码组-->编号-->存于leftInfo.coodInfo第4列(Inf为图中有但未进行编号的控制点,其原因是一组中控制点由于佛像遮挡有丢失,无法编号)

5. 优化外部参数

三. 操作说明