一、项目名称：全景图像拼接
二、语言、库和开发环境：C++、opencv249、VS2013
三、参考
2：计算机视觉中的多视图几何
四、包含原理：
1、特征点
2、2D射影变换
3、相机参数
4、拼接缝的查找
4.3、相机参数评估
从图6可以看出，真实物体通过小孔映射到成像平面上，小孔到成像平面的距离称为焦距f。在成像平面上的图像是镜像倒立的，所以为了研究方便，在小孔和物体之间定义一个虚拟成像平面（在后面，我们把该平面也称为成像平面），它与小孔的距离也为焦距，则两个成像平面的图像大小是相同的，但虚拟成像平面上的图像与原物体的方向是一样的。
我们以小孔为坐标原点建立一个三维直角坐标系XYZ，坐标原点C称为相机的光心。成像平面xy平行于XY，并且距离光心C为f，其中Z轴定义为光轴，它与成像平面xy的交点为P，因此CP＝f。设空间中的任一点Q的坐标为（X, Y, Z），该点映射到成像平面的点q的坐标为(x, y,f)。
我们以小孔为坐标原点建立一个三维直角坐标系XYZ（如图7所示），坐标原点C称为相机的光心。成像平面xy平行于XY，并且距离光心C为f，其中Z轴定义为光轴，它与成像平面xy的交点为P，因此CP＝f。设空间中的任一点Q的坐标为（X, Y, Z），该点映射到成像平面的点q的坐标为(x, y,f)。
考虑人站在原地，没有发生位移，只是旋转相机角度
这部分结果对前两行、前两列调整至系数为1，便于计算。
当有多幅图像需要拼接为一幅图像时，是要以其中一幅图像为基准，其他图像都要旋转到该基准图像平面上的，所以就需要找到基准平面。这里用到的算法为最大生成树算法。
待拼接图像的排列是无序的，而且我们事先是不知道它们之间的关系的，我们只知道它们之间的单应矩阵，而单应矩阵是由图像间的内点计算评估得到的。由此我们可以构造一个无向图G，G的节点为图像，G的边为内点数，然后利用并查集在该G中得到一棵最大生成树。
我们把树的中心节点作为基准图像。中心节点的确定方法为：计算每一个节点到所有叶节点的距离，把其中的最大值作为该节点的值；然后选择这些值中最小者作为中心节点。这里的距离指的是节点间的节点数。如图8(c)所示，节点A和C为中心节点。中心节点可能是1个，也可能是2个，如果是2个，则选择其中一个即可。
基于以上方法，我们得到了相机的内外参数，但这样得到的参数忽略了多个图像间的约束，而且会产生累计误差。这时，我们就需要用到光束平差法（Bundle Adjustment）来精确化相机参数。光束指的是相机光心“发射”出来的光束（或射线），它透过相片达到物点，因此相片中的点应该和物点处于一个光束线上，但当两者不共线时，我们就需要对结构和视角参数进行调节，以达到最优解甚至共线的目的。最优化一般采用前文介绍过的LM算法。
应用于光束平差法的LM算法，误差指标函数可以有两个，一个是重映射误差，另一个是射线发散误差。
代码：相机参数评估
HomographyBasedEstimator estimator; //相机参数评估
vector<CameraParams> cameras;
estimator(features, pairwise_matches, cameras);
基本公式：
uv1=f0ppx0afppy001R|tXYZ1
结果：当前两张图像对应两个相机，每一相机由10个参数组成
具体代码：
1、输入数据：两张图像得到的features和匹配后的pairwise_matches特征关系，
1）Features的格式为：
图像0中有976个特征点
图像1中有2325个特征点
其中，第0张图片的特征keypoints
其中，第0张图片的第0个特征点内容：
其中，第0张图片的特征描述子descriptors：为976*64的矩阵，代表976个特征，每个特征占64个字节。
- pairwise_matches的格式
共两张图像，可以得到4中对应关系，图像i和j的匹配为：
pairwise_matches[i*num_images + j]
其中，以i=0，j=1为例，每一匹配信息为：
其中，matches为280表示图像0到图像1之间的映射有280组对应特征点，以第0个match为例：
queryIdx : 查询点的索引，指0图片中特征库的索引
trainIdx : 被查询到的点的索引，指1图片中特征库的索引
2、参数初始化
1）初始化 f0，f1，由H中得到F0、F1
原理：
在计算f时，为减少参数量，不考虑ppx、ppy和a三个参数，
这部分结果对前两行、前两列调整至系数为1，便于计算。
Step 1: 计算和确定f1
d1 = h[6] * h[7]; //式48的分母
d2 = (h[7] - h[6]) * (h[7] + h[6]); //式47的分母
v1 = -(h[0] * h[1] + h[3] * h[4]) / d1; //式48
v2 = (h[0] * h[0] + h[3] * h[3] - h[1] * h[1] - h[4] * h[4]) / d2; //式47
if (v1 > 0 && v2 > 0) f1 = sqrt(std::abs(d1) > std::abs(d2) ? v1 : v2);
在v1和v2中选择分母较大的作为f1
Step 2：计算和确定f0
d1 = h[0] * h[3] + h[1] * h[4]; //式44的分母
d2 = h[0] * h[0] + h[1] * h[1] - h[3] * h[3] - h[4] * h[4]; //式43的分母
v1 = -h[2] * h[5] / d1; //式44
v2 = (h[5] * h[5] - h[2] * h[2]) / d2; //式43
if (v1 < v2) std::swap(v1, v2); //使v1不小于v2
if (v1 > 0 && v2 > 0) f0 = sqrt(std::abs(d1) > std::abs(d2) ? v1 : v2);
在v1和v2中选择分母较大的作为f0
Step 3：计算确定f
取中值
all_focals.push_back(sqrt(f0 * f1)); //式49，把焦距存入all_focals内
median = all_focals[all_focals.size() / 2];
Step 4：计算最大生成树，确定中心点图像
Step 5：计算旋转矩阵
原理：
cameras[edge.from].R=I
cameras[edge.to].aspect 为1
cameras[edge.from].ppx; //表示式33的cx 为0
cameras[edge.from].ppy; //表示式33的cy 为0
int pair_idx = edge.from * num_images + edge.to; //表示匹配点对的索引
//构造式51中的参数K0
Mat_<double> K_from = Mat::eye(3, 3, CV_64F); //初始化
K_from(0, 0) = cameras[edge.from].focal; //表示式33的fu
//表示式33的fv
K_from(1, 1) = cameras[edge.from].focal * cameras[edge.from].aspect;
K_from(0, 2) = cameras[edge.from].ppx; //表示式33的cx 为0
K_from(1, 2) = cameras[edge.from].ppy; //表示式33的cy 为0
//构造式51中的参数K1
Mat_<double> K_to = Mat::eye(3, 3, CV_64F); //初始化
K_to(0, 0) = cameras[edge.to].focal; //表示式33的fu
K_to(1, 1) = cameras[edge.to].focal * cameras[edge.to].aspect; //表示式33的fv
K_to(0, 2) = cameras[edge.to].ppx; //表示式33的cx
K_to(1, 2) = cameras[edge.to].ppy; //表示式33的cy
Mat R = K_from.inv() * pairwise_matches[pair_idx].H.inv() * K_to; //式51
//式52，可见CameraParams变量中R实际存储的是相机旋转矩阵变量的逆
cameras[edge.to].R = cameras[edge.from].R * R;
Step 6：初始化参数 ppx，ppy
cameras[i].ppx += 0.5 * features[i].img_size.width;
cameras[i].ppy += 0.5 * features[i].img_size.height;
Step 7：BA
原理：
损失：
- 参数初始化
一个相机共7个参数：f、PPX、PPY、a、R（3个）
cam_params_.at<double>(i * 7, 0) = cameras[i].focal; //初始化fu
cam_params_.at<double>(i * 7 + 1, 0) = cameras[i].ppx; //初始化cx
cam_params_.at<double>(i * 7 + 2, 0) = cameras[i].ppy; //初始化cy
cam_params_.at<double>(i * 7 + 3, 0) = cameras[i].aspect; //初始化α
svd(cameras[i].R, SVD::FULL_UV);
Mat R = svd.u * svd.vt;
if (determinant(R) < 0)
R *= -1;
Mat rvec;
Rodrigues(R, rvec); //旋转矩阵R由Rodrigues算法得到旋转向量rvec
CV_Assert(rvec.type() == CV_32F);
cam_params_.at<double>(i * 7 + 4, 0) = rvec.at<float>(0, 0); //初始化rx
cam_params_.at<double>(i * 7 + 5, 0) = rvec.at<float>(1, 0); //初始化ry
cam_params_.at<double>(i * 7 + 6, 0) = rvec.at<float>(2, 0); //初始化rz
2) 计算迭代损失函数：
3)得到单应，计算重投影误差：
4)LM优化
PS：在 opencv实现的BA的 LM算法迭代优化中存在疑惑？？
opencv中的原代码：
每次迭代中_jac 和 _err 赋值为0，在算法中是否起到更新和迭代的作用？？？
改正后的代码
扫一扫
3404
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
