近景摄影测量编程与实习

前言

记录编程中的问题

影像获取

摄影方式的选择

摄影方式主要包括两种：正直摄影方式和交向摄影方式。正直摄影方式中，像片对两像片的主光轴彼此平行，且垂直于摄影基线；交向摄影方式中，像片对两像片的主光轴大体位于同一平面但彼此不平行，且不垂直于摄影基线。

受限于检校场的大小，为了能拍摄到更多的标志，获得更大的重叠度，可以选择交向摄影方式。同时，为了保证交向摄影的精度，还应该注意交会角的大小。

焦距的选择

检校场留有的摄影空间较小，摄影距离受限，为了拍摄的影像尽可能包含较多的标志，应当选择较小的焦距。但是，根据交向摄影的精度估算式，为了拍摄出比例尺较大的影像以提高精度，焦距也不应过小。

调焦方式的选择

为了保持畸变参数和内方位元素不变，以简化后续计算，拍摄过程中，不应进行调焦。因而相机的调焦方式应当选择手动调焦。

曝光程序的选择

曝光程序主要有快门优先和光圈优先两种。快门优先即先确定曝光时间，再测定光圈号数；光圈优先即先确定光圈号数，再测定曝光时间。

由于控制景深需要设置光圈的大小，因而要选择光圈优先曝光程序。

景深的控制

近景摄影测量中，获取清晰的高质量影像至关重要，特别需要关注的是景深的正确选择。景深 $\Delta D$ 是在给定光圈和模糊圈大小条件下被摄影空间获得清晰构象的深度范围，也就是后景距 $D_2$ 与前景距 $D_1$ 的差值。

$\begin{equation} D_1=\cfrac{HD}{H+D} \end{equation}$

$\begin{equation} D_2=\cfrac{HD}{H-D} \end{equation}$

$\begin{equation} \Delta D=D_2-D_1 \end{equation}$

其中

$\begin{equation} H=\cfrac{F^2}{kE} \end{equation}$

$D$ 代表调焦距， $H$ 代表超焦点距离， $F$ 为摄影机焦距， $k$ 为光圈号数， $E$ 为模糊圆直径。

要保证拍摄物体清晰，就需要控制景深，或者说，使被摄物体位于前景距和后景距之间。由公式(1)~(4)可知，调焦距 $D$ 、摄影机焦距 $F$ 、光圈号数 $k$ 、模糊圆直径 $E$ 是影响前景距与后景距的因素。而由于检校场的摄影空间较小， $D$ 与 $F$ 调节范围受限，而 $E$ 又是基本确定的值。因此，为了使标志能够成像清晰，需要调节光圈号数 $k$ 。这也是曝光程序选择光圈优先的原因。

分析公式(1)~(4)，当 $H>D$ 时， $k$ 越大， $H$ 越小，则前景距 $D_1$ 越小且后景距 $D_2$ 越大。所以，若检校场的标志不能清晰成像，可以选择加大光圈号数，从而使标志点位于前景距与后景距之间，获得清晰成像。

曝光时间的控制

为了确定曝光时间，可以使用经验法、曝光表测定法，试片法，或推算比较法。在曝光程序为光圈优先的情况下，实验使用的相机可以通过测光功能，自动安置曝光时间。但是，在实际的近景摄影测量工作中，有时需要有人工的略微修正。

相机方位元素的记录

拍摄影像对的目的之一就是为了通过后方交会，计算相机的内外方位元素和畸变参数。由于后方交会解求的是非线性方程，需要提供方位元素的近似值，所以在摄影过程中，应当记录方位元素的近似值。

像点坐标量测

圆形识别

采用的是边缘检测，并拟合椭圆的方法。这种方法简单有效，易于操作，但是对参数的设置要求较高。

主要函数：

GaussianBlur：进行高斯模糊，减小噪声。噪声会影响Canny算子提取边缘的效果。

确实影响很大，不做高斯模糊有的很明显的圆都识别不出来，还困扰了我很久。
Canny：提取边缘。
findContours：边缘追踪，得到轮廓。
contourArea：计算轮廓包围的面积，可以控制识别圆形的大小。

但是要注意contourArea算的面积是有误差的。
arcLength：计算轮廓的周长，结合面积可以计算圆形度。
fitEllipse：根据轮廓拟合椭圆。

cvtColor(focusImg, grayImg, CV_BGR2GRAY);

//边缘提取
GaussianBlur(grayImg, grayImg, Size(5, 5), 0);

Mat edgeImage;
Canny(grayImg, edgeImage, 30, 70);

//边缘追踪
vector<vector<Point>> contours;		// 轮廓

findContours(edgeImage, contours, CV_RETR_LIST, CV_CHAIN_APPROX_NONE);


double area = 0;
double perimeter = 0;
Mat cirImg = Mat::zeros(edgeImage.size(), CV_8UC3);
size_t count = 0;
Point center = Point(-999, -999);

//找出并筛选圆形
for (auto it1 = contours.begin(); it1 != contours.end(); )
{
    //拟合点个数筛选
    count = it1->size();
    if (count < 6)
    {
        it1 = contours.erase(it1);
        continue;
    }

    //面积筛选与圆度筛选
    //这里面积的阈值如何确定？（可以用户输入）
    area = contourArea(*it1);
    perimeter = arcLength(*it1, true);
    if (area <minArea || area>maxArea)
    {
        it1 = contours.erase(it1);
        continue;
    }
    if (4 * 3.1415926 * area / pow(perimeter, 2) < roundness)
    {
        it1 = contours.erase(it1);
        continue;
    }


    Point temp = fitEllipse(*it1).center;
    if (sqrtf(powf((center.x - temp.x), 2) + powf((center.y - temp.y), 2)) < distance)
    {
        it1 = contours.erase(it1);
        continue;
    }
    else
    {
        center = temp;
    }


    //迭代器后移
    it1++;
}

区域框选

需要的主要函数是鼠标事件函数setMouseCallback。比较麻烦的就是鼠标事件获取坐标之后如何传递的问题，我现在还没有理清。~~（我用的是全局变量，不得不说全局变量是真的丑）~~。

Rect select;				//框选的矩形
bool select_flag = false;	//框选状态
Point origin;				//框选的起点

//鼠标事件的回调函数，控制框选区域的是select(全局变量)
void onMouse(int event, int x, int y, int, void*)
{
	if (event == EVENT_LBUTTONDOWN)
	{
		select_flag = true;			//左键按下赋真值
		origin = Point(x, y);		//保存左键单击的坐标
		select = Rect(x, y, 0, 0);	

	}
	else if (event == EVENT_LBUTTONUP)
	{
		select_flag = false;
	}
	if (select_flag)
	{
		select.x = MIN(origin.x, x);					//鼠标按下开始到弹起这段时间实时计算所选矩形框左上角点坐标
		select.y = MIN(origin.y, y);
		select.width = abs(x - origin.x);				//算矩形宽度和高度
		select.height = abs(y - origin.y);
		select &= Rect(0, 0, preImg.cols, preImg.rows);	//保证所选矩形框在图片区域之内
	}
}

代码可以参考：基础学习笔记之opencv(5)：实现鼠标选定矩形框 - tornadomeet - 博客园 (cnblogs.com)

键盘操作

是很简单的输入流与waitKey()控制。输入输出流用于在控制台进行交互。waitKey()则读取按键信息。

奇怪的地方在于，waitKey()没办法识别delete、shift等等的按键，不知道是什么原因。

此外waitKey()有一个重要的功能就是让imshow()能显示出图像，特别是在与控制台交互输入点号信息的时候。如果不用waitKey(30)，窗口就会显示灰色而不是图像。

单像空间后方交会

单像空间后方交会是根据一张像片覆盖的一定数量控制点的物方空间坐标及其像点坐标，以像点坐标作为观测值，按共线条件方程解算该像片的内、外方位元素以及其它附加参数的过程。

坐标系转换

根据实习提供的控制点物方坐标数据，可以推断出控制点物方坐标系各坐标轴的朝向。如图1左图所示，物方坐标数据的坐标系是一个左手坐标系。一方面，左手坐标系不符合通常情况下共线条件方程式的推导结果。另一方面，不恰当的坐标系，会使外方位角元素的初始值（记录的水平角与竖直角）难以转换到选择的坐标系中。采用如图 1 右图所示的物方坐标系，使 Z 轴与摄影方向一致。既能保证右手坐标系，又能将其视为从空中向地面的竖直摄影，方便地使用推导的公式，还可以简单地将水平角与竖直角转换为外方位角元素。

此外，由于像点坐标量测使用了 OpenCV，量测坐标的平面坐标系与常用的平面坐标系定义不同，如图 2 左图所示。也需要将其转换为如图 2 右图所示的坐标系并配合转换后的物方空间坐标系使用。同理，在输出结果时，要将坐标变换回原始像平面坐标系中的坐标。

[{"url":"https://cdn.zerolacqua.top/images/blogs/近景摄影测量编程与实习/1.png","alt":"图1 物方空间坐标系的变换"},{"url":"https://cdn.zerolacqua.top/images/blogs/近景摄影测量编程与实习/2.png","alt":"图2 像平面坐标系的转换"}]

误差方程式

误差方程式参见教材。

需要注意以下几点：

畸变差系数的导数前的符号，这是由共线条件方程中中畸变差的形式决定的。
误差方程式解算的是各未知数的改正数

精度计算

单像空间后方交会解算中计算的是内精度。

根据平差原理知识，单位权中误差如下式所示。

$\begin{equation} \sigma_0=\sqrt{\cfrac{V^\mathrm{T}V}{r}} \end{equation}$

其中 $r$ 为多余观测数。

则由单位权中误差，可以计算出各未知数的中误差

$\begin{equation} m_X=\sigma_0\sqrt{Diag(Q_{XX})} \end{equation}$

直接线性变换

直接线性变换解法是建立像点的“坐标仪坐标”和相应物点的物方空间坐标直接的线性关系的解法。

坐标系转换

与单像空间后方交会中所述一致。

近似值的计算

近似值的计算参见教材。

计算 $l$ 系数的近似值时，还需要计算像主点坐标 $x_0,y_0$ ，这是因为在直接线性变换中使用了畸变系数，需要使用像主点坐标 $x_0,y_0$ 。

误差方程式

误差方程式参见教材。

需要注意以下几点：

误差方程式解算的是各未知数本身

精度计算

直接线性变换解算解算 $l$ 系数的过程，计算的是内精度。原理与单像空间后方交会中所述一致。

直接线性变换解算待定点物方坐标的过程计算的是内精度和外精度。内精度计算原理与单像空间后方交会中所述一致。外精度则通过下式计算：

$\begin{equation} \begin{aligned} m_X&=\sqrt{\cfrac{1}{n}\sum(X-\hat{X})^2}\\ m_Y&=\sqrt{\cfrac{1}{n}\sum(Y-\hat{Y})^2}\\ m_Z&=\sqrt{\cfrac{1}{n}\sum(Z-\hat{Z})^2}\\ m_p&=\sqrt{m_X^2+m_Y^2+m_Z^2}\\ \end{aligned} \end{equation}$

方位元素解算

方位元素的计算参见教材。

对于 $b_1,b_2$ 的计算方法，教材上似乎没有明确给出，个人认为其计算方法如下：

$\begin{equation} \begin{aligned} b_1&=\cfrac{\gamma_3l_2+b_3x_0+b_2f_x\tan{\mathrm{d}\beta}}{f_x}\\ b_2&=\cfrac{\gamma_3(l_6+l_{10}y_0)(1+\mathrm{d}s)\cos{\mathrm{d}\beta}}{f_x} \end{aligned} \end{equation}$