核线影像生成

??《航摄影像密集匹配的研究进展与展望》一文中介绍：航摄影像的密集匹配可以对单个立体影像对、也可以对多度重叠的多视影像序列进行。在获得影像的内外方位元素或者影像的相对方位元素的前提下，实施过程大体分为核线影像生成、匹配代价计算、匹配代价聚合、视差计算与精华、三维点云生成几个主要步骤。其一般流程可描述为下图：

??由上面的描述，我们可以知道，核线影像生成是航摄影像密集匹配的第一步,目的是利用核线相关的该概念将沿x、y方向搜索同名像点的二维相关问题转化为沿核线搜索同名像点。为了生成核线影像，需要影像的内外方位元素或者相对方位元素。这里我们可以分别用相对方位元素或外方位元素来生成核线影像。

>使用相对方位元素或外方位元素生成核线影像(以两幅影像为例)

需要的数据：相对方位元素(可以通过两幅影像特征匹配得到的匹配点对计算得到)、两幅影像的各四个框标点的扫描坐标和框标坐标(数字影像四角的像素坐标即为扫描坐标)。

过程：1.求解内定向参数；2.确定核线影像的范围；3.建立“水平”核线影像的像素坐标与像平面坐标系之间的变换关系；4.生成核线影像。

代码实现思路：
1.求解内定参数：这个过程主要是建立倾斜影像的像平面坐标系与其扫描坐标系之间的关系。原始倾斜影像的扫描坐标系与其像平面坐标系之间的变换关系可采用仿射变换。设影像的扫描坐标系为(c,r),影像像平面坐标系为 (x,y),则根据仿射变换公式 ,有如下关系 :

其中 a₀, a₁, a₂, b₀,b₁ , b₂ 为影像内定向的6个内定向参数。一般利用影像四个框标点坐标(x,y)组成误差方程,根据最小二乘平差求解这6个内定向参数。

2.建立“水平”核线影像的像素坐标与像平面坐标系之间的变换关系:通过建立两者之间的变换关系可以确定核线影像的范围，以及可以生成核线影像。

像空间坐标转换到像空间辅助坐标系的关系为：

像空间辅助坐标系转换核线影像像空间的关系：

原始像空间坐标系到核线影像的像空间的转换关系则为：

反过来，核线影像像空间坐标系到原始影像的关系则为：

代码实现：

int main()
{
	//读取原始影像对
	Mat img2 = imread("2.tif");
	Mat img1 = imread("1.tif");
	

	if (img1.empty() || img2.empty())
	{
		cout << "image open error!!" << endl;
		return -1;
	}
	//获取原始影像对宽高
	// 获取影像大小
	int n_Width_Left = img1.cols;
	int n_Height_Left = img1.rows;
	int n_Width_Right = img2.cols;
	int n_Height_Right = img2.rows;
	/* 以下是内定向部分 */
	//即通过左右倾斜像片的各自四个框标点的扫描坐标和框标坐标完成内定向6个参数的确定。
	double delta = xxx;
	double f = xxx;
	double x0 = xxx;
	double y0 = xxx;
	//扫描坐标系下框表坐标的行列号
	double kuangbiao[8] = { 0,0,n_Width_Left,0,n_Width_Left,n_Height_Left,0,n_Height_Left};
	//以框表点在像平面坐标系下的坐标
	double scanpos[8] = { xxx,xxx,xxx,xxx,xxx,xxx,xxx,xxx};
	//左右像片内定向参数
	double np[6] = { 0.0 };
	Internal_ori(scanpos, kuangbiao, np, delta);

	
	//已有数据获取核线的像空间坐标转换到原始影像的像空间坐标
	//相对方位元素
	double BX = xxx;
	double BY = xxx;
	double BZ = xxx;
	double B = sqrt(BX * BX + BY * BY + BZ * BZ);
	double m_T = atan2(BY, BX);
	double m_V = asin(BZ / B);

	Mat Rtemp(3, 3, CV_64FC1);
	Rtemp.at<double>(0, 0) = cos(m_T) * cos(m_V);
	Rtemp.at<double>(0, 1) = sin(m_T) * cos(m_V);
	Rtemp.at<double>(0, 2) = sin(m_V);
	Rtemp.at<double>(1, 0) = -sin(m_T);
	Rtemp.at<double>(1, 1) = cos(m_T);
	Rtemp.at<double>(1, 2) = 0;
	Rtemp.at<double>(2, 0) = -cos(m_T) * sin(m_V);
	Rtemp.at<double>(2, 1) = -sin(m_T) * sin(m_V);
	Rtemp.at<double>(2, 2) = cos(m_V);

	double fai = xxx;
	double w = xxx;
	double k = xxx;
	Mat Rotationr(3, 3, CV_64FC1);
	Rotationr.at<double>(0, 0) = cos(fai) * cos(k) - sin(fai) * sin(w) * sin(k);
	Rotationr.at<double>(0, 1) = -cos(fai) * sin(k) - sin(fai) * sin(w) * cos(k);
	Rotationr.at<double>(0, 2) = -sin(fai) * cos(w);
	Rotationr.at<double>(1, 0) = cos(w) * sin(k);
	Rotationr.at<double>(1, 1) = cos(w) * cos(k);
	Rotationr.at<double>(1, 2) = -sin(w);
	Rotationr.at<double>(2, 0) = sin(fai) * cos(k) + cos(fai) * sin(w) * sin(k);
	Rotationr.at<double>(2, 1) = -sin(fai) * sin(k) + cos(fai) * sin(w) * cos(k);
	Rotationr.at<double>(2, 2) = cos(fai) * cos(w);

	Mat Rleft = Rtemp;
	Mat Rright = Rtemp * Rotationr;

	int HxWidth_left, HxHeight_left;
	int HxWidth_right, HxHeight_right;
	//获取核线影像大小
	getHxImageSize(n_Width_Left, n_Height_Left, f, delta, Rleft, HxWidth_left, HxHeight_left, np, x0, y0);
	getHxImageSize(n_Width_Right, n_Height_Right, f, delta, Rright, HxWidth_right, HxHeight_right, np, x0, y0);

	getQxImageSize(HxWidth_left, HxHeight_left, f, delta, Rleft, HxWidth_left, HxHeight_left, np, x0, y0);
	getQxImageSize(HxWidth_right, HxHeight_right, f, delta, Rright, HxWidth_right, HxHeight_right, np, x0, y0);

	int HxImageHeight = HxHeight_left > HxHeight_right ? HxHeight_left : HxHeight_right;
	int HxImageWidth = HxWidth_left > HxWidth_right ? HxWidth_left : HxWidth_right;

	/* 开始制作核线影像 */
	Mat HxImageLeft = Mat(HxImageWidth, HxImageHeight, img1.type());
	Mat HxImageRight = Mat(HxImageWidth, HxImageHeight, img2.type());
	double u1, v1, x1, y1;
	double u2, v2, x2, y2;
	double c1, r1 ,c2, r2;
	for (int i = 0;i < HxImageHeight;i++)
	{
		for (int j = 0;j < HxImageWidth;j++)
		{
			scanning2plane(i, j, u1, v1, np, x0, y0, delta);
			hx2qx(u1, v1, x1, y1, Rleft, f);
			plane2scanning(c1, r1, x1, y1, np, x0, y0, delta);
			
			scanning2plane(i, j, u2, v2, np, x0, y0, delta);
			hx2qx(u2, v2, x2, y2, Rright, f);
			plane2scanning(c2, r2, x2, y2, np, x0, y0, delta);

			if ((int(c1) < 0.0) || (int(r1) < 0.0) || (r1 > n_Height_Left-1) || (c1 > n_Width_Left-1))
			{
				if (img1.channels() == 3)
				{
					HxImageLeft.at<cv::Vec3b>(j, i)[0] = 0;
					HxImageLeft.at<cv::Vec3b>(j, i)[1] = 0;
					HxImageLeft.at<cv::Vec3b>(j, i)[2] = 0;
				}
			}
			else// 双线性内插
			{
				if (img1.channels() == 3)
				{
					HxImageLeft.at<cv::Vec3b>(j, i)[0] = BilinearInterpolation(&img1, r1, c1, 0);
					HxImageLeft.at<cv::Vec3b>(j, i)[1] = BilinearInterpolation(&img1, r1, c1, 1);
					HxImageLeft.at<cv::Vec3b>(j, i)[2] = BilinearInterpolation(&img1, r1, c1, 2);
				}
			}
			if ((int(c2) < 0.0) || (int(r2) < 0.0) || (r2 > n_Height_Right - 1) || (c2 > n_Width_Right - 1))
			{
	            if (img2.channels() == 3)
				{
					HxImageRight.at<cv::Vec3b>(j, i)[0] = 0;
					HxImageRight.at<cv::Vec3b>(j, i)[1] = 0;
					HxImageRight.at<cv::Vec3b>(j, i)[2] = 0;
				}
			}
			else// 双线性内插
			{
				if (img2.channels() == 3)
				{
					HxImageRight.at<cv::Vec3b>(j, i)[0] = BilinearInterpolation(&img2, r2, c2, 0);
					HxImageRight.at<cv::Vec3b>(j, i)[1] = BilinearInterpolation(&img2, r2, c2, 1);
					HxImageRight.at<cv::Vec3b>(j, i)[2] = BilinearInterpolation(&img2, r2, c2, 2);
				}
			}

		}
		//cout << i << endl;
	}

	imwrite("hx1.tif", HxImageLeft);
	imwrite("hx2.tif", HxImageRight);


	return 0;
}

极线校正

??极线校正是通过对两个相机进行旋转，并重新定义新的像平面，让极线对共线且平行于像平面的某条坐标轴(通常是水平轴)，该操作同时建立了新的立体像对。校正目的是对两幅图像的二维搜索变成一维，节省计算量，排除虚假匹配点。校正过程是将相机在数学上对准到同一观察平面上，使得相机上像素行是严格对齐的。
在这里插入图片描述

>Fusiello校正法

为了将摄像机极点变换到无穷远使得极线水平对准，我们创建了一个旋转矩阵，首先是X轴的旋转，我们将X轴旋转到与基线相同的方向即旋转后新X轴就是O_l—O_r，旋转向量：

确定旋转后的新Y轴，只需满足与e₁正交即可：选择与主光轴e₁相交的方向：

新的Z轴，与e₁,e₂垂直且满足右手定则：

最终得到的旋转矩阵：

新的旋转矩阵：

校正后，baseline T计算：

下式中是指右相机光心在左相机坐标系下的坐标（也等于与右相机光心[0,0,0]的距离），求解公式的意义是左相机光心在右相机校正后的坐标系下的坐标（也等于与右相机校正后的坐标系下光心[0,0,0]的距离）。故baseline就是x轴方向的距离。

>Fusiello校正法伪代码

 function [A,R,t] = art(P)
  % ART: factorize a PPM as  P=A*[R;t]
  
  Q = inv(P(1:3, 1:3));
  [U,B] = qr(Q);
  
  R = inv(U);
  t = B*P(1:3,4);
  A = inv(B);
  A = A ./A(3,3);
  function [T1,T2,Pn1,Pn2] = rectify(Po1,Po2)
  
  % RECTIFY: 计算校正矩阵 
  
  % 分解原始的PPMs
  [A1,R1,t1] = art(Po1);
  [A2,R2,t2] = art(Po2);
  
  % 光学中心(不变)
  c1 = - inv(Po1(:,1:3))*Po1(:,4);
  c2 = - inv(Po2(:,1:3))*Po2(:,4);
  
  % 新的x轴（基线方向）
  v1 = (c1-c2);
  % 新的y轴（与新的x和旧的z正交） 
  v2 = cross(R1(3,:)',v1);
  % 新的z轴（与基线和y正交）
  v3 = cross(v1,v2);
  
  % 新的外参
  R = [v1'/norm(v1)
       v2'/norm(v2)
       v3'/norm(v3)];
  % translation is left unchanged
  
  % 新的内参（任意）
  A = (A1 + A2)./2;
  A(1,2)=0; % no skew
  
  % 新的投影矩阵 
  Pn1 = A * [R -R*c1 ];
  Pn2 = A * [R -R*c2 ];
  
  % 校正图像变换 
  T1 = Pn1(1:3,1:3)* inv(Po1(1:3,1:3));
  T2 = Pn2(1:3,1:3)* inv(Po2(1:3,1:3));