狠狠撸

一、stitching_detail 程序运行流程
1.命令行调用程序，输入源图像以及程序的参数
2.特征点检测，判断是使用 surf 还是 orb，默认是 surf。
3.对图像的特征点进行匹配，使用最近邻和次近邻方法，将两个最优的匹配的置信度
保存下来。
4.对图像进行排序以及将置信度高的图像保存到同一个集合中，删除置信度比较低的
图像间的匹配，得到能正确匹配的图像序列。这样将置信度高于门限的所有匹配合并到一个
集合中。
5.对所有图像进行相机参数粗略估计，然后求出旋转矩阵
6.使用光束平均法进一步精准的估计出旋转矩阵。
7.波形校正，水平或者垂直
8.拼接
9.融合，多频段融合，光照补偿，
二、stitching_detail 程序接口介绍
img1 img2 img3 输入图像
--preview 以预览模式运行程序，比正常模式要快，但输出图像分辨率低，拼接的分辨
率 compose_megapix 设置为 0.6
--try_gpu (yes|no) 是否使用 GPU(图形处理器)，默认为 no
/* 运动估计参数 */
--work_megapix <--work_megapix <float>> 图像匹配的分辨率大小，图像的面积尺寸
变为 work_megapix*100000，默认为 0.6
--features (surf|orb) 选择 surf 或者 orb 算法进行特征点计算，默认为 surf
--match_conf <float> 特征点检测置信等级，最近邻匹配距离与次近邻匹配距离的比值，
surf 默认为 0.65，orb 默认为 0.3
--conf_thresh <float> 两幅图来自同一全景图的置信度，默认为 1.0
--ba (reproj|ray) 光束平均法的误差函数选择，默认是 ray 方法
--ba_refine_mask (mask) ---------------
--wave_correct (no|horiz|vert) 波形校验水平，垂直或者没有默认是 horiz
--save_graph <file_name> 将匹配的图形以点的形式保存到文件中， Nm 代表匹配的
数量，NI 代表正确匹配的数量，C 表示置信度
/*图像融合参数:*/
--warp
(plane|cylindrical|spherical|fisheye|stereographic|compressedPlaneA2B1|compressedPla
neA1.5B1|compressedPlanePortraitA2B1

程序参数读入分析，将程序运行是输入的参数以及需要拼接的图像读入内存，检查图像
是否多于 2 张。
[cpp] view plaincopy
1. int retval = parseCmdArgs(argc, argv);
2. if (retval)
3. return retval;
4.
5. // Check if have enough images
6. int num_images = static_cast<int>(img_names.size());
7. if (num_images < 2)
8. {
9. LOGLN("Need more images");
10. return -1;
11. }
2.特征点检测
判断选择是 surf 还是 orb 特征点检测（默认是 surf）以及对图像进行预处理（尺寸缩
放），然后计算每幅图形的特征点，以及特征点描述子
2.1 计算 work_scale，将图像 resize 到面积在 work_megapix*10^6 以下，
（work_megapix 默认是 0.6）
1. work_scale = min(1.0, sqrt(work_megapix * 1e6 / full_img.size().area()));
1. resize(full_img, img, Size(), work_scale, work_scale);
图像大小是 740*830，面积大于 6*10^5，所以计算出 work_scale = 0.98,然后对图像
resize。

2.2 计算 seam_scale，也是根据图像的面积小于 seam_megapix*10^6，
（seam_megapix 默认是 0.1），seam_work_aspect 目前还没用到
1. seam_scale = min(1.0, sqrt(seam_megapix * 1e6 / full_img.size().area()));
2. seam_work_aspect = seam_scale / work_scale; //seam_megapix = 0.1 seam_work_a
spect = 0.69
2.3 计算图像特征点，以及计算特征点描述子，并将 img_idx 设置为 i。
1. (*finder)(img, features[i]);//matcher.cpp 348
2. features[i].img_idx = i;
特征点描述的结构体定义如下：
1. struct detail::ImageFeatures
2. Structure containing image keypoints and descriptors.
3. struct CV_EXPORTS ImageFeatures
4. {
5. int img_idx;//
6. Size img_size;
7. std::vector<KeyPoint> keypoints;
8. Mat descriptors;
9. };
2.4 将源图像 resize 到 seam_megapix*10^6，并存入 image[]中

1. resize(full_img, img, Size(), seam_scale, seam_scale);
2. images[i] = img.clone();
3.图像匹配
对任意两副图形进行特征点匹配，然后使用查并集法，将图片的匹配关系找出，并删
除那些不属于同一全景图的图片。
3.1 使用最近邻和次近邻匹配，对任意两幅图进行特征点匹配。
1. vector<MatchesInfo> pairwise_matches;//Structure containing information abou
t matches between two images.
2. BestOf2NearestMatcher matcher(try_gpu, match_conf);//最近邻和次近邻法
3. matcher(features, pairwise_matches);//对每两个图片进行 matcher，20-》
400 matchers.cpp 502
介绍一下 BestOf2NearestMatcher 函数：
1. //Features matcher which finds two best matches for each feature and leav
es the best one only if the ratio between descriptor distances is greater th
an the threshold match_conf.
2. detail::BestOf2NearestMatcher::BestOf2NearestMatcher(bool try_use_gpu=fals
e,float match_conf=0.3f,
3. intnum_matches_thresh1=6,
int num_matches_thresh2=6)
4. Parameters: try_use_gpu – Should try to use GPU or not
5. match_conf – Match distances ration threshold
6. num_matches_thresh1 – Minimum number of matches required for the 2D projecti
ve
7. transform estimation used in the inliers classification step
8. num_matches_thresh2 – Minimum number of matches required for the 2D projecti
ve
9. transform re-estimation on inliers
函数的定义（只是设置一下参数，属于构造函数）：
1. BestOf2NearestMatcher::BestOf2NearestMatcher(bool try_use_gpu, float match_c
onf, int num_matches_thresh1, int num_matches_thresh2)
2. {
3. #ifdef HAVE_OPENCV_GPU
4. if (try_use_gpu && getCudaEnabledDeviceCount() > 0)
5. impl_ = new GpuMatcher(match_conf);
6. else

7. #else
8. (void)try_use_gpu;
9. #endif
10. impl_ = new CpuMatcher(match_conf);
11.
12. is_thread_safe_ = impl_->isThreadSafe();
13. num_matches_thresh1_ = num_matches_thresh1;
14. num_matches_thresh2_ = num_matches_thresh2;
15. }
以及 MatchesInfo 的结构体定义：
1. Structure containing information about matches between two images. It’s assu
med that there is a homography between those images.
2. struct CV_EXPORTS MatchesInfo
3. {
4. MatchesInfo();
5. MatchesInfo(const MatchesInfo &other);
6. const MatchesInfo& operator =(const MatchesInfo &other);
7. int src_img_idx, dst_img_idx; // Images indices (optional)
8. std::vector<DMatch> matches;
9. std::vector<uchar> inliers_mask; // Geometrically consistent matches
mask
10. int num_inliers; // Number of geometrically consistent matches
11. Mat H; // Estimated homography
12. double confidence; // Confidence two images are from the same panora
ma
13. };
求出图像匹配的结果如下（具体匹配参见 sift 特征点匹配），任意两幅图都进行匹配
（3*3=9），其中 1-》2 和 2-》1 只计算了一次，以 1-》2 为准，：

3.2 根据任意两幅图匹配的置信度，将所有置信度高于 conf_thresh（默认是 1.0）的
图像合并到一个全集中。
我们通过函数的参数 save_graph 打印匹配结果如下（我稍微修改了一下输出）：
1. "outimage101.jpg" -- "outimage102.jpg"[label="Nm=866, Ni=637, C=2.37864"];
Nm 代表匹配的数量，NI 代表正确匹配的数量，C 表示置信度
通过查并集方法，查并集介绍请参见博文：
http://blog.csdn.net/skeeee/article/details/20471949
1. vector<int> indices = leaveBiggestComponent(features, pairwise_matches, conf
_thresh);//将置信度高于门限的所有匹配合并到一个集合中
2. vector<Mat> img_subset;
3. vector<string> img_names_subset;
4. vector<Size> full_img_sizes_subset;
5. for (size_t i = 0; i < indices.size(); ++i)
6. {
7. img_names_subset.push_back(img_names[indices[i]]);
8. img_subset.push_back(images[indices[i]]);
9. full_img_sizes_subset.push_back(full_img_sizes[indices[i]]);
10. }
11.
12. images = img_subset;
13. img_names = img_names_subset;
14. full_img_sizes = full_img_sizes_subset;
4.根据单应性矩阵粗略估计出相机参数（焦距）
4.1 焦距参数的估计
根据前面求出的任意两幅图的匹配，我们根据两幅图的单应性矩阵 H，求出符合条
件的 f，（4 副图，16 个匹配，求出 8 个符合条件的 f），然后求出 f 的均值或者中值当成
所有图形的粗略估计的 f。
1. estimateFocal(features, pairwise_matches, focals);
函数的主要源码如下：
1. for (int i = 0; i < num_images; ++i)

2. {
3. for (int j = 0; j < num_images; ++j)
4. {
5. const MatchesInfo &m = pairwise_matches[i*num_images + j];
6. if (m.H.empty())
7. continue;
8. double f0, f1;
9. bool f0ok, f1ok;
10. focalsFromHomography(m.H, f0, f1, f0ok, f1ok);//Tries to estimate focal leng
ths from the given homography 79
11. //under the assumption that the camera undergoes rotations around its centre
only.
12. if (f0ok && f1ok)
13. all_focals.push_back(sqrt(f0 * f1));
14. }
15. }
其中函数 focalsFromHomography 的定义如下：
1. Tries to estimate focal lengths from the given homography
2. under the assumption that the camera undergoes rotations around its cent
re only.
3. Parameters
4. H – Homography.
5. f0 – Estimated focal length along X axis.
6. f1 – Estimated focal length along Y axis.
7. f0_ok – True, if f0 was estimated successfully, false otherwise.
8. f1_ok – True, if f1 was estimated successfully, false otherwise.
函数的源码：
1. void focalsFromHomography(const Mat& H, double &f0, double &f1, bool &f0_ok,
bool &f1_ok)
2. {
3. CV_Assert(H.type() == CV_64F && H.size() == Size(3, 3));//Checks a condi
tion at runtime and throws exception if it fails
4.
5. const double* h = reinterpret_cast<const double*>(H.data);//强制类型转
换
6.
7. double d1, d2; // Denominators
8. double v1, v2; // Focal squares value candidates

9.
10. //具体的计算过程有点看不懂啊
11. f1_ok = true;
12. d1 = h[6] * h[7];
13. d2 = (h[7] - h[6]) * (h[7] + h[6]);
14. v1 = -(h[0] * h[1] + h[3] * h[4]) / d1;
15. v2 = (h[0] * h[0] + h[3] * h[3] - h[1] * h[1] - h[4] * h[4]) / d2;
16. if (v1 < v2) std::swap(v1, v2);
17. if (v1 > 0 && v2 > 0) f1 = sqrt(std::abs(d1) > std::abs(d2) ? v1 : v2);
18. else if (v1 > 0) f1 = sqrt(v1);
19. else f1_ok = false;
20.
21. f0_ok = true;
22. d1 = h[0] * h[3] + h[1] * h[4];
23. d2 = h[0] * h[0] + h[1] * h[1] - h[3] * h[3] - h[4] * h[4];
24. v1 = -h[2] * h[5] / d1;
25. v2 = (h[5] * h[5] - h[2] * h[2]) / d2;
26. if (v1 < v2) std::swap(v1, v2);
27. if (v1 > 0 && v2 > 0) f0 = sqrt(std::abs(d1) > std::abs(d2) ? v1 : v2);
28. else if (v1 > 0) f0 = sqrt(v1);
29. else f0_ok = false;
30. }
求出的焦距有 8 个
求出的焦距取中值或者平均值，然后就是所有图片的焦距。
并构建 camera 参数，将矩阵写入 camera：
1. cameras.assign(num_images, CameraParams());
3. cameras[i].focal = focals[i];

4.2 根据匹配的内点构建最大生成树，然后广度优先搜索求出根节点，并求出 camera
的 R 矩阵，K 矩阵以及光轴中心
camera 其他参数：
aspect = 1.0，ppx，ppy 目前等于 0，最后会赋值成图像中心点的。
K 矩阵的值：
1. Mat CameraParams::K() const
2. {
3. Mat_<double> k = Mat::eye(3, 3, CV_64F);
4. k(0,0) = focal; k(0,2) = ppx;
5. k(1,1) = focal * aspect; k(1,2) = ppy;
6. return k;
7. }
R 矩阵的值：
1. void operator ()(const GraphEdge &edge)
2. {
3. int pair_idx = edge.from * num_images + edge.to;
4.
5. Mat_<double> K_from = Mat::eye(3, 3, CV_64F);
6. K_from(0,0) = cameras[edge.from].focal;
7. K_from(1,1) = cameras[edge.from].focal * cameras[edge.from].aspect;
8. K_from(0,2) = cameras[edge.from].ppx;
9. K_from(0,2) = cameras[edge.from].ppy;
10.
11. Mat_<double> K_to = Mat::eye(3, 3, CV_64F);
12. K_to(0,0) = cameras[edge.to].focal;
13. K_to(1,1) = cameras[edge.to].focal * cameras[edge.to].aspect;
14. K_to(0,2) = cameras[edge.to].ppx;

15. K_to(0,2) = cameras[edge.to].ppy;
16.
17. Mat R = K_from.inv() * pairwise_matches[pair_idx].H.inv() * K_to;
18. cameras[edge.to].R = cameras[edge.from].R * R;
19. }
光轴中心的值
2. {
3. cameras[i].ppx += 0.5 * features[i].img_size.width;
4. cameras[i].ppy += 0.5 * features[i].img_size.height;
5. }
5.使用 Bundle Adjustment 方法对所有图片进行相机参数校正
Bundle Adjustment（光束法平差）算法主要是解决所有相机参数的联合。这是全景拼
接必须的一步，因为多个成对的单应性矩阵合成全景图时，会忽略全局的限制，造成累积误
差。因此每一个图像都要加上光束法平差值，使图像被初始化成相同的旋转和焦距长度。
光束法平差的目标函数是一个具有鲁棒性的映射误差的平方和函数。即每一个特征点
都要映射到其他的图像中，计算出使误差的平方和最小的相机参数。具体的推导过程可以参
见 Automatic Panoramic Image Stitching using Invariant Features.pdf 的第五章，本文只介
绍一下 opencv 实现的过程（完整的理论和公式暂时还没看懂，希望有人能一起交流）
opencv 中误差描述函数有两种如下：（opencv 默认是 BundleAdjusterRay ）：
1. BundleAdjusterReproj // BundleAdjusterBase(7, 2)//最小投影误差平方和
2. Implementation of the camera parameters refinement algorithm which minimizes
sum of the reprojection error squares
3.
4. BundleAdjusterRay // BundleAdjusterBase(4, 3)//最小特征点与相机中心点的距离
和
5. Implementation of the camera parameters refinement algorithm which minimizes
sum of the distances between the
6. rays passing through the camera center and a feature.
5.1 首先计算 cam_params_的值：
1. setUpInitialCameraParams(cameras);

函数主要源码：
1. cam_params_.create(num_images_ * 4, 1, CV_64F);
2. SVD svd;//奇异值分解
3. for (int i = 0; i < num_images_; ++i)
4. {
5. cam_params_.at<double>(i * 4, 0) = cameras[i].focal;
6.
7. svd(cameras[i].R, SVD::FULL_UV);
8. Mat R = svd.u * svd.vt;
9. if (determinant(R) < 0)
10. R *= -1;
11.
12. Mat rvec;
13. Rodrigues(R, rvec);
14. CV_Assert(rvec.type() == CV_32F);
15. cam_params_.at<double>(i * 4 + 1, 0) = rvec.at<float>(0, 0);
18. }
计算 cam_params_的值，先初始化 cam_params(num_images_*4,1,CV_64F);
cam_params_[i*4+0] = cameras[i].focal;
cam_params_后面 3 个值，是 cameras[i].R 先经过奇异值分解，然后对 u*vt 进行
Rodrigues 运算，得到的 rvec 第一行 3 个值赋给 cam_params_。
奇异值分解的定义：
在矩阵 M 的奇异值分解中 M = UΣV*
·U 的列(columns)组成一套对 M 的正交"输入"或"分析"的基向量。这些向量是 MM*的特
征向量。
·V 的列(columns)组成一套对 M的正交"输出"的基向量。这些向量是 M*M的特征向量。
·Σ 对角线上的元素是奇异值，可视为是在输入与输出间进行的标量的"膨胀控制"。这
些是 M*M 及 MM*的奇异值，并与 U 和 V 的行向量相对应。
5.2 删除置信度小于门限的匹配对
1. // Leave only consistent image pairs
2. edges_.clear();
3. for (int i = 0; i < num_images_ - 1; ++i)
4. {
5. for (int j = i + 1; j < num_images_; ++j)
6. {

7. const MatchesInfo& matches_info = pairwise_matches_[i * num_images_
+ j];
8. if (matches_info.confidence > conf_thresh_)
9. edges_.push_back(make_pair(i, j));
10. }
11. }
5.3 使用 LM 算法计算 camera 参数。
首先初始化 LM 的参数（具体理论还没有看懂）
1. //计算所有内点之和
2. for (size_t i = 0; i < edges_.size(); ++i)
3. total_num_matches_ += static_cast<int>(pairwise_matches[edges_[i].fi
rst * num_images_ +
4. edges_[i].se
cond].num_inliers);
5.
6. CvLevMarq solver(num_images_ * num_params_per_cam_,
7. total_num_matches_ * num_errs_per_measurement_,
8. term_criteria_);
9.
10. Mat err, jac;
11. CvMat matParams = cam_params_;
12. cvCopy(&matParams, solver.param);
13.
14. int iter = 0;
15. for(;;)//类似于 while（1）,但是比 while（1）效率高
16. {
17. const CvMat* _param = 0;
18. CvMat* _jac = 0;
19. CvMat* _err = 0;
20.
21. bool proceed = solver.update(_param, _jac, _err);
22.
23. cvCopy(_param, &matParams);
24.
25. if (!proceed || !_err)
26. break;
27.
28. if (_jac)
29. {
30. calcJacobian(jac); //构造雅阁比行列式
31. CvMat tmp = jac;

32. cvCopy(&tmp, _jac);
33. }
34.
35. if (_err)
36. {
37. calcError(err);//计算 err
38. LOG_CHAT(".");
39. iter++;
40. CvMat tmp = err;
41. cvCopy(&tmp, _err);
42. }
43. }
计算 camera
1. obtainRefinedCameraParams(cameras);//Gets the refined camera parameters.
函数源代码：
1. void BundleAdjusterRay::obtainRefinedCameraParams(vector<CameraParams> &came
ras) const
2. {
4. {
5. cameras[i].focal = cam_params_.at<double>(i * 4, 0);
6.
7. Mat rvec(3, 1, CV_64F);
8. rvec.at<double>(0, 0) = cam_params_.at<double>(i * 4 + 1, 0);
11. Rodrigues(rvec, cameras[i].R);
12.
13. Mat tmp;
14. cameras[i].R.convertTo(tmp, CV_32F);
15. cameras[i].R = tmp;
16. }
17. }
求出根节点，然后归一化旋转矩阵 R
1. // Normalize motion to center image
2. Graph span_tree;

3. vector<int> span_tree_centers;
4. findMaxSpanningTree(num_images_, pairwise_matches, span_tree, span_tree_cent
ers);
5. Mat R_inv = cameras[span_tree_centers[0]].R.inv();
7. cameras[i].R = R_inv * cameras[i].R;
6 波形校正
前面几节把相机旋转矩阵计算出来，但是还有一个因素需要考虑，就是由于拍摄者拍摄
图片的时候不一定是水平的，轻微的倾斜会导致全景图像出现飞机曲线，因此我们要对图像
进行波形校正，主要是寻找每幅图形的“上升向量”（up_vector），使他校正成
。
波形校正的效果图
opencv 实现的源码如下（也是暂时没看懂，囧）

1. <span style="white-space:pre"> </span>vector<Mat> rmats;
2. for (size_t i = 0; i < cameras.size(); ++i)
3. rmats.push_back(cameras[i].R);
4. waveCorrect(rmats, wave_correct);//Tries to make panorama more horizo
ntal (or vertical).
5. for (size_t i = 0; i < cameras.size(); ++i)
6. cameras[i].R = rmats[i];
其中 waveCorrect(rmats, wave_correct)源码如下：
1. void waveCorrect(vector<Mat> &rmats, WaveCorrectKind kind)
2. {
3. LOGLN("Wave correcting...");
4. #if ENABLE_LOG
5. int64 t = getTickCount();
6. #endif
7.
8. Mat moment = Mat::zeros(3, 3, CV_32F);
9. for (size_t i = 0; i < rmats.size(); ++i)
10. {
11. Mat col = rmats[i].col(0);
12. moment += col * col.t();//相机 R 矩阵第一列转置相乘然后相加
13. }
14. Mat eigen_vals, eigen_vecs;
15. eigen(moment, eigen_vals, eigen_vecs);//Calculates eigenvalues and eigen
vectors of a symmetric matrix.
16.
17. Mat rg1;
18. if (kind == WAVE_CORRECT_HORIZ)
19. rg1 = eigen_vecs.row(2).t();//如果是水平校正，去特征向量的第三行
20. else if (kind == WAVE_CORRECT_VERT)
21. rg1 = eigen_vecs.row(0).t();//如果是垂直校正，特征向量第一行
22. else
23. CV_Error(CV_StsBadArg, "unsupported kind of wave correction");
24.
25. Mat img_k = Mat::zeros(3, 1, CV_32F);
27. img_k += rmats[i].col(2);//R 函数第 3 列相加
28. Mat rg0 = rg1.cross(img_k);//rg1 与 img_k 向量积
29. rg0 /= norm(rg0);//归一化？
30.
31. Mat rg2 = rg0.cross(rg1);
32.
33. double conf = 0;

34. if (kind == WAVE_CORRECT_HORIZ)
35. {
37. conf += rg0.dot(rmats[i].col(0));//Computes a dot-product of two
vectors.数量积
38. if (conf < 0)
39. {
40. rg0 *= -1;
41. rg1 *= -1;
42. }
43. }
44. else if (kind == WAVE_CORRECT_VERT)
45. {
47. conf -= rg1.dot(rmats[i].col(0));
48. if (conf < 0)
49. {
50. rg0 *= -1;
51. rg1 *= -1;
52. }
53. }
54.
55. Mat R = Mat::zeros(3, 3, CV_32F);
56. Mat tmp = R.row(0);
57. Mat(rg0.t()).copyTo(tmp);
58. tmp = R.row(1);
60. tmp = R.row(2);
62.
64. rmats[i] = R * rmats[i];
65.
66. LOGLN("Wave correcting, time: " << ((getTickCount() - t) / getTickFreque
ncy()) << " sec");
67. }
7.单应性矩阵变换
由图像匹配，Bundle Adjustment 算法以及波形校验，求出了图像的相机参数以及旋转
矩阵，接下来就对图形进行单应性矩阵变换，亮度的增量补偿以及多波段融合（图像金字塔）。
首先介绍的就是单应性矩阵变换：
源图像的点（x，y，z=1），图像的旋转矩阵 R，图像的相机参数矩阵 K，经过变换后
的同一坐标（x_,y_,z_），然后映射到球形坐标（u，v，w），他们之间的关系如下：

1. void SphericalProjector::mapForward(float x, float y, float &u, float &v)
2. {
3. float x_ = r_kinv[0] * x + r_kinv[1] * y + r_kinv[2];
4. float y_ = r_kinv[3] * x + r_kinv[4] * y + r_kinv[5];
5. float z_ = r_kinv[6] * x + r_kinv[7] * y + r_kinv[8];
6.
7. u = scale * atan2f(x_, z_);
8. float w = y_ / sqrtf(x_ * x_ + y_ * y_ + z_ * z_);
9. v = scale * (static_cast<float>(CV_PI) - acosf(w == w ? w : 0));
10. }
根据映射公式，对图像的上下左右四个边界求映射后的坐标，然后确定变换后图像的
左上角和右上角的坐标，
如果是球形拼接，则需要再加上判断（暂时还没研究透）：
1. float tl_uf = static_cast<float>(dst_tl.x);
2. float tl_vf = static_cast<float>(dst_tl.y);
3. float br_uf = static_cast<float>(dst_br.x);
4. float br_vf = static_cast<float>(dst_br.y);
5.
6. float x = projector_.rinv[1];
7. float y = projector_.rinv[4];
8. float z = projector_.rinv[7];
9. if (y > 0.f)
10. {
11. float x_ = (projector_.k[0] * x + projector_.k[1] * y) / z + projector_.
k[2];
12. float y_ = projector_.k[4] * y / z + projector_.k[5];
13. if (x_ > 0.f && x_ < src_size.width && y_ > 0.f && y_ < src_size.height)
14. {

15. tl_uf = min(tl_uf, 0.f); tl_vf = min(tl_vf, static_cast<float>(CV_PI
* projector_.scale));
16. br_uf = max(br_uf, 0.f); br_vf = max(br_vf, static_cast<float>(CV_PI
* projector_.scale));
17. }
18. }
19.
20. x = projector_.rinv[1];
21. y = -projector_.rinv[4];
22. z = projector_.rinv[7];
23. if (y > 0.f)
24. {
25. float x_ = (projector_.k[0] * x + projector_.k[1] * y) / z + projector_.
k[2];
26. float y_ = projector_.k[4] * y / z + projector_.k[5];
27. if (x_ > 0.f && x_ < src_size.width && y_ > 0.f && y_ < src_size.height)
28. {
29. tl_uf = min(tl_uf, 0.f); tl_vf = min(tl_vf, static_cast<float>(0));
30. br_uf = max(br_uf, 0.f); br_vf = max(br_vf, static_cast<float>(0));
31. }
32. }
然后利用反投影将图形反投影到变换的图像上，像素计算默认是二维线性插值。
反投影的公式：
1. void SphericalProjector::mapBackward(float u, float v, float &x, float &y)
2. {
3. u /= scale;
4. v /= scale;
5.
6. float sinv = sinf(static_cast<float>(CV_PI) - v);
7. float x_ = sinv * sinf(u);
8. float y_ = cosf(static_cast<float>(CV_PI) - v);
9. float z_ = sinv * cosf(u);
10.
11. float z;
12. x = k_rinv[0] * x_ + k_rinv[1] * y_ + k_rinv[2] * z_;
13. y = k_rinv[3] * x_ + k_rinv[4] * y_ + k_rinv[5] * z_;
14. z = k_rinv[6] * x_ + k_rinv[7] * y_ + k_rinv[8] * z_;
15.

16. if (z > 0) { x /= z; y /= z; }
17. else x = y = -1;
18. }
然后将反投影求出的 x，y 值写入 Mat 矩阵 xmap 和 ymap 中
1. xmap.create(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, CV_32F);
2. ymap.create(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, CV_32F);
3.
4. float x, y;
5. for (int v = dst_tl.y; v <= dst_br.y; ++v)
6. {
7. for (int u = dst_tl.x; u <= dst_br.x; ++u)
8. {
9. projector_.mapBackward(static_cast<float>(u), static_cast<float>(
v), x, y);
10. xmap.at<float>(v - dst_tl.y, u - dst_tl.x) = x;
11. ymap.at<float>(v - dst_tl.y, u - dst_tl.x) = y;
12. }
13. }
最后使用 opencv 自带的 remap 函数将图像重新绘制:
1. remap(src, dst, xmap, ymap, interp_mode, border_mode);//重映射，xmap，yamp 分别
是 u，v 反投影对应的 x，y 值，默认是双线性插值
8.光照补偿
图像拼接中，由于拍摄的图片有可能因为光圈或者光线的问题，导致相邻图片重叠区
域出现亮度差，所以在拼接时就需要对图像进行亮度补偿，（opencv 只对重叠区域进行了
亮度补偿，这样会导致图像融合处虽然光照渐变，但是图像整体的光强没有柔和的过渡）。
首先，将所有图像，mask 矩阵进行分块（大小在 32*32 附近）。
1. for (int img_idx = 0; img_idx < num_images; ++img_idx)
2. {
3. Size bl_per_img((images[img_idx].cols + bl_width_ - 1) / bl_width_,
4. (images[img_idx].rows + bl_height_ - 1) / bl_height_);
5. int bl_width = (images[img_idx].cols + bl_per_img.width - 1) / bl_per_
img.width;

6. int bl_height = (images[img_idx].rows + bl_per_img.height - 1) / bl_pe
r_img.height;
7. bl_per_imgs[img_idx] = bl_per_img;
8. for (int by = 0; by < bl_per_img.height; ++by)
9. {
10. for (int bx = 0; bx < bl_per_img.width; ++bx)
11. {
12. Point bl_tl(bx * bl_width, by * bl_height);
13. Point bl_br(min(bl_tl.x + bl_width, images[img_idx].cols),
14. min(bl_tl.y + bl_height, images[img_idx].rows));
15.
16. block_corners.push_back(corners[img_idx] + bl_tl);
17. block_images.push_back(images[img_idx](Rect(bl_tl, bl_br)));
18. block_masks.push_back(make_pair(masks[img_idx].first(Rect(bl_t
l, bl_br)),
19. masks[img_idx].second));
20. }
21. }
22. }
然后，求出任意两块图像的重叠区域的平均光强，
1. //计算每一块区域的光照均值 sqrt(sqrt（R）+sqrt(G)+sqrt(B))
2. //光照均值是对称矩阵，所以一次循环计算两个光照值，（i，j），与（j，i）
4. {
5. for (int j = i; j < num_images; ++j)
6. {
7. Rect roi;
8. //判断 image[i]与 image[j]是否有重叠部分
9. if (overlapRoi(corners[i], corners[j], images[i].size(), images[
j].size(), roi))
10. {
11. subimg1 = images[i](Rect(roi.tl() - corners[i], roi.br() - c
orners[i]));
12. subimg2 = images[j](Rect(roi.tl() - corners[j], roi.br() - c
orners[j]));
13.
14. submask1 = masks[i].first(Rect(roi.tl() - corners[i], roi.br
() - corners[i]));
15. submask2 = masks[j].first(Rect(roi.tl() - corners[j], roi.br
() - corners[j]));

16. intersect = (submask1 == masks[i].second) & (submask2 == mas
ks[j].second);
17.
18. N(i, j) = N(j, i) = max(1, countNonZero(intersect));
19.
20. double Isum1 = 0, Isum2 = 0;
21. for (int y = 0; y < roi.height; ++y)
22. {
23. const Point3_<uchar>* r1 = subimg1.ptr<Point3_<uchar> >(
y);
24. const Point3_<uchar>* r2 = subimg2.ptr<Point3_<uchar> >(
y);
25. for (int x = 0; x < roi.width; ++x)
26. {
27. if (intersect(y, x))
28. {
29. Isum1 += sqrt(static_cast<double>(sqr(r1[x].x) +
sqr(r1[x].y) + sqr(r1[x].z)));
30. Isum2 += sqrt(static_cast<double>(sqr(r2[x].x) +
sqr(r2[x].y) + sqr(r2[x].z)));
31. }
32. }
33. }
34. I(i, j) = Isum1 / N(i, j);
35. I(j, i) = Isum2 / N(i, j);
36. }
37. }
38. }
建立方程，求出每个光强的调整系数
1. Mat_<double> A(num_images, num_images); A.setTo(0);
2. Mat_<double> b(num_images, 1); b.setTo(0);//beta*N(i,j)
4. {
5. for (int j = 0; j < num_images; ++j)
6. {
7. b(i, 0) += beta * N(i, j);
8. A(i, i) += beta * N(i, j);
9. if (j == i) continue;
10. A(i, i) += 2 * alpha * I(i, j) * I(i, j) * N(i, j);
11. A(i, j) -= 2 * alpha * I(i, j) * I(j, i) * N(i, j);
12. }
13. }

14.
15. solve(A, b, gains_);//求方程的解 A*gains = B
gains_原理分析:
num_images ：表示图像分块的个数，零 num_images = n
N 矩阵，大小 n*n，N(i,j)表示第 i 幅图像与第 j 幅图像重合的像素点数，N(i,j)=N(j,i)
I 矩阵，大小 n*n，I(i,j)与 I(j,i)表示第 i，j 块区域重合部分的像素平均值，I(i，j)是重合区域
中 i 快的平均亮度值，
参数 alpha 和 beta，默认值是 alpha=0.01，beta=100.
A 矩阵，大小 n*n，公式图片不全
b 矩阵，大小 n*1，
然后根据求解矩阵
gains_矩阵，大小 1*n，A*gains = B
然后将 gains_进行线性滤波
1. Mat_<float> ker(1, 3);
2. ker(0,0) = 0.25; ker(0,1) = 0.5; ker(0,2) = 0.25;
3.
4. int bl_idx = 0;
5. for (int img_idx = 0; img_idx < num_images; ++img_idx)
6. {
7. Size bl_per_img = bl_per_imgs[img_idx];
8. gain_maps_[img_idx].create(bl_per_img);
9.
10. for (int by = 0; by < bl_per_img.height; ++by)

11. for (int bx = 0; bx < bl_per_img.width; ++bx, ++bl_idx)
12. gain_maps_[img_idx](by, bx) = static_cast<float>(gains[bl_idx]
);
13. //用分解的核函数对图像做卷积。首先，图像的每一行与一维的核 kernelX 做卷积；然后，运算
结果的每一列与一维的核 kernelY 做卷积
14. sepFilter2D(gain_maps_[img_idx], gain_maps_[img_idx], CV_32F, ker, ker
);
15. sepFilter2D(gain_maps_[img_idx], gain_maps_[img_idx], CV_32F, ker, ker
);
16. }
然后构建一个 gain_maps 的三维矩阵，gain_main[图像的个数][图像分块的行数][图像
分块的列数]，然后对没一副图像的 gain 进行滤波。
9.Seam Estimation
缝隙估计有 6 种方法，默认就是第三种方法，seam_find_type == "gc_color"，该方法
是利用最大流方法检测。
1. if (seam_find_type == "no")
2. seam_finder = new detail::NoSeamFinder();//Stub seam estimator which
does nothing.
3. else if (seam_find_type == "voronoi")
4. seam_finder = new detail::VoronoiSeamFinder();//Voronoi diagram-base
d seam estimator. 泰森多边形缝隙估计
5. else if (seam_find_type == "gc_color")
6. {
8. if (try_gpu && gpu::getCudaEnabledDeviceCount() > 0)
9. seam_finder = new detail::GraphCutSeamFinderGpu(GraphCutSeamFind
erBase::COST_COLOR);
10. else
11. #endif
12. seam_finder = new detail::GraphCutSeamFinder(GraphCutSeamFinderB
ase::COST_COLOR);//Minimum graph cut-based seam estimator
13. }
14. else if (seam_find_type == "gc_colorgrad")
15. {
17. if (try_gpu && gpu::getCudaEnabledDeviceCount() > 0)
18. seam_finder = new detail::GraphCutSeamFinderGpu(GraphCutSeamFind
erBase::COST_COLOR_GRAD);
19. else

20. #endif
21. seam_finder = new detail::GraphCutSeamFinder(GraphCutSeamFinderB
ase::COST_COLOR_GRAD);
22. }
23. else if (seam_find_type == "dp_color")
24. seam_finder = new detail::DpSeamFinder(DpSeamFinder::COLOR);
25. else if (seam_find_type == "dp_colorgrad")
26. seam_finder = new detail::DpSeamFinder(DpSeamFinder::COLOR_GRAD);
27. if (seam_finder.empty())
28. {
29. cout << "Can't create the following seam finder '" << seam_find_type
<< "'n";
30. return 1;
31. }
程序解读可参见：
http://blog.csdn.net/chlele0105/article/details/12624541
http://blog.csdn.net/zouxy09/article/details/8534954
http://blog.csdn.net/yangtrees/article/details/7930640
10.多波段融合
由于以前进行处理的图片都是以 work_scale（3.1 节有介绍）进行缩放的，所以图像
的内参，corner（同一坐标后的顶点），mask（融合的掩码）都需要重新计算。以及将之
前计算的光照增强的 gain 也要计算进去。
1. // Read image and resize it if necessary
2. full_img = imread(img_names[img_idx]);
3. if (!is_compose_scale_set)
4. {
5. if (compose_megapix > 0)
6. compose_scale = min(1.0, sqrt(compose_megapix * 1e6 / full_img
.size().area()));
7. is_compose_scale_set = true;
8.
9. // Compute relative scales
10. //compose_seam_aspect = compose_scale / seam_scale;
11. compose_work_aspect = compose_scale / work_scale;
12.
13. // Update warped image scale
14. warped_image_scale *= static_cast<float>(compose_work_aspect);
15. warper = warper_creator->create(warped_image_scale);
16.

17. // Update corners and sizes
19. {
20. // Update intrinsics
21. cameras[i].focal *= compose_work_aspect;
22. cameras[i].ppx *= compose_work_aspect;
23. cameras[i].ppy *= compose_work_aspect;
24.
25. // Update corner and size
26. Size sz = full_img_sizes[i];
27. if (std::abs(compose_scale - 1) > 1e-1)
28. {
29. sz.width = cvRound(full_img_sizes[i].width * compose_scale
);//取整
30. sz.height = cvRound(full_img_sizes[i].height * compose_sca
le);
31. }
32.
33. Mat K;
34. cameras[i].K().convertTo(K, CV_32F);
35. Rect roi = warper->warpRoi(sz, K, cameras[i].R);//Returns Proj
ected image minimum bounding box
36. corners[i] = roi.tl();//! the top-left corner
37. sizes[i] = roi.size();//! size of the real buffer
38. }
39. }
40. if (abs(compose_scale - 1) > 1e-1)
41. resize(full_img, img, Size(), compose_scale, compose_scale);
42. else
43. img = full_img;
44. full_img.release();
45. Size img_size = img.size();
46.
47. Mat K;
48. cameras[img_idx].K().convertTo(K, CV_32F);
49.
50. // Warp the current image
51. warper->warp(img, K, cameras[img_idx].R, INTER_LINEAR, BORDER_REFLECT,
img_warped);
52. // Warp the current image mask
53. mask.create(img_size, CV_8U);
54. mask.setTo(Scalar::all(255));
55. warper->warp(mask, K, cameras[img_idx].R, INTER_NEAREST, BORDER_CONSTA
NT, mask_warped);

56. // Compensate exposure
57. compensator->apply(img_idx, corners[img_idx], img_warped, mask_warped)
;//光照补偿
58. img_warped.convertTo(img_warped_s, CV_16S);
59. img_warped.release();
60. img.release();
61. mask.release();
62.
63. dilate(masks_warped[img_idx], dilated_mask, Mat());
64. resize(dilated_mask, seam_mask, mask_warped.size());
65. mask_warped = seam_mask & mask_warped;
对图像进行光照补偿，就是将对应区域乘以 gain：
1. //块光照补偿
2. void BlocksGainCompensator::apply(int index, Point /*corner*/, Mat &image, c
onst Mat &/*mask*/)
3. {
4. CV_Assert(image.type() == CV_8UC3);
5.
6. Mat_<float> gain_map;
7. if (gain_maps_[index].size() == image.size())
8. gain_map = gain_maps_[index];
9. else
10. resize(gain_maps_[index], gain_map, image.size(), 0, 0, INTER_LINEAR
);
11.
12. for (int y = 0; y < image.rows; ++y)
13. {
14. const float* gain_row = gain_map.ptr<float>(y);
15. Point3_<uchar>* row = image.ptr<Point3_<uchar> >(y);
16. for (int x = 0; x < image.cols; ++x)
17. {
18. row[x].x = saturate_cast<uchar>(row[x].x * gain_row[x]);
19. row[x].y = saturate_cast<uchar>(row[x].y * gain_row[x]);
20. row[x].z = saturate_cast<uchar>(row[x].z * gain_row[x]);
21. }
22. }
23. }
进行多波段融合，首先初始化 blend，确定 blender 的融合的方式，默认是多波段融合
MULTI_BAND，以及根据 corners 顶点和图像的大小确定最终全景图的尺寸。

1. <span> </span>//初始化 blender
2. if (blender.empty())
3. {
4. blender = Blender::createDefault(blend_type, try_gpu);
5. Size dst_sz = resultRoi(corners, sizes).size();//计算最后图像的大
小
6. float blend_width = sqrt(static_cast<float>(dst_sz.area())) * bl
end_strength / 100.f;
7. if (blend_width < 1.f)
8. blender = Blender::createDefault(Blender::NO, try_gpu);
9. else if (blend_type == Blender::MULTI_BAND)
10. {
11. MultiBandBlender* mb = dynamic_cast<MultiBandBlender*>(stati
c_cast<Blender*>(blender));
12. mb->setNumBands(static_cast<int>(ceil(log(blend_width)/log(2
.)) - 1.));
13. LOGLN("Multi-band blender, number of bands: " << mb->numBand
s());
14. }
15. else if (blend_type == Blender::FEATHER)
16. {
17. FeatherBlender* fb = dynamic_cast<FeatherBlender*>(static_ca
st<Blender*>(blender));
18. fb->setSharpness(1.f/blend_width);
19. LOGLN("Feather blender, sharpness: " << fb->sharpness());
20. }
21. blender->prepare(corners, sizes);//根据 corners 顶点和图像的大小确定
最终全景图的尺寸
22. }
然后对每幅图图形构建金字塔，层数可以由输入的系数确定，默认是 5 层。
先对顶点以及图像的宽和高做处理，使其能被 2^num_bands 除尽，这样才能将进行
向下采样 num_bands 次，首先从源图像 pyrDown 向下采样，在由最底部的图像 pyrUp 向
上采样，把对应层得上采样和下采样的相减，就得到了图像的高频分量，存储到每一个金字
塔中。然后根据 mask，将每幅图像的各层金字塔分别写入最终的金字塔层 src_pyr_laplace
中。
最后将各层的金字塔叠加，就得到了最终完整的全景图。
1. blender->feed(img_warped_s, mask_warped, corners[img_idx]);//将图像写入金字塔
中
源码：

1. void MultiBandBlender::feed(const Mat &img, const Mat &mask, Point tl)
2. {
3. CV_Assert(img.type() == CV_16SC3 || img.type() == CV_8UC3);
4. CV_Assert(mask.type() == CV_8U);
5.
6. // Keep source image in memory with small border
7. int gap = 3 * (1 << num_bands_);
8. Point tl_new(max(dst_roi_.x, tl.x - gap),
9. max(dst_roi_.y, tl.y - gap));
10. Point br_new(min(dst_roi_.br().x, tl.x + img.cols + gap),
11. min(dst_roi_.br().y, tl.y + img.rows + gap));
12.
13. // Ensure coordinates of top-left, bottom-right corners are divided by (
1 << num_bands_).
14. // After that scale between layers is exactly 2.
15. //
16. // We do it to avoid interpolation problems when keeping sub-images only
. There is no such problem when
17. // image is bordered to have size equal to the final image size, but thi
s is too memory hungry approach.
18. //将顶点调整成能被 2^num_bank 次方除尽的值
19. tl_new.x = dst_roi_.x + (((tl_new.x - dst_roi_.x) >> num_bands_) << num_
bands_);
20. tl_new.y = dst_roi_.y + (((tl_new.y - dst_roi_.y) >> num_bands_) << num_
bands_);
21. int width = br_new.x - tl_new.x;
22. int height = br_new.y - tl_new.y;
23. width += ((1 << num_bands_) - width % (1 << num_bands_)) % (1 << num_ban
ds_);
24. height += ((1 << num_bands_) - height % (1 << num_bands_)) % (1 << num_b
ands_);
25. br_new.x = tl_new.x + width;
26. br_new.y = tl_new.y + height;
27. int dy = max(br_new.y - dst_roi_.br().y, 0);
28. int dx = max(br_new.x - dst_roi_.br().x, 0);
29. tl_new.x -= dx; br_new.x -= dx;
30. tl_new.y -= dy; br_new.y -= dy;
31.
32. int top = tl.y - tl_new.y;
33. int left = tl.x - tl_new.x;
34. int bottom = br_new.y - tl.y - img.rows;
35. int right = br_new.x - tl.x - img.cols;
36.

37. // Create the source image Laplacian pyramid
38. Mat img_with_border;
39. copyMakeBorder(img, img_with_border, top, bottom, left, right,
40. BORDER_REFLECT);//给图像设置一个边界，BORDER_REFLECT 边界颜色
任意
41. vector<Mat> src_pyr_laplace;
42. if (can_use_gpu_ && img_with_border.depth() == CV_16S)
43. createLaplacePyrGpu(img_with_border, num_bands_, src_pyr_laplace);
44. else
45. createLaplacePyr(img_with_border, num_bands_, src_pyr_laplace);//创建
高斯金字塔，每一层保存的全是高频信息
46.
47. // Create the weight map Gaussian pyramid
48. Mat weight_map;
49. vector<Mat> weight_pyr_gauss(num_bands_ + 1);
50.
51. if(weight_type_ == CV_32F)
52. {
53. mask.convertTo(weight_map, CV_32F, 1./255.);//将 mask 的 0，255 归一化成
0,1
54. }
55. else// weight_type_ == CV_16S
56. {
57. mask.convertTo(weight_map, CV_16S);
58. add(weight_map, 1, weight_map, mask != 0);
59. }
60.
61. copyMakeBorder(weight_map, weight_pyr_gauss[0], top, bottom, left, right
, BORDER_CONSTANT);
62.
63. for (int i = 0; i < num_bands_; ++i)
64. pyrDown(weight_pyr_gauss[i], weight_pyr_gauss[i + 1]);
65.
66. int y_tl = tl_new.y - dst_roi_.y;
67. int y_br = br_new.y - dst_roi_.y;
68. int x_tl = tl_new.x - dst_roi_.x;
69. int x_br = br_new.x - dst_roi_.x;
70.
71. // Add weighted layer of the source image to the final Laplacian pyramid
layer
72. if(weight_type_ == CV_32F)
73. {
74. for (int i = 0; i <= num_bands_; ++i)
75. {

76. for (int y = y_tl; y < y_br; ++y)
77. {
78. int y_ = y - y_tl;
79. const Point3_<short>* src_row = src_pyr_laplace[i].ptr<Point
3_<short> >(y_);
80. Point3_<short>* dst_row = dst_pyr_laplace_[i].ptr<Point3_<sh
ort> >(y);
81. const float* weight_row = weight_pyr_gauss[i].ptr<float>(y_)
;
82. float* dst_weight_row = dst_band_weights_[i].ptr<float>(y);
83.
84. for (int x = x_tl; x < x_br; ++x)
85. {
86. int x_ = x - x_tl;
87. dst_row[x].x += static_cast<short>(src_row[x_].x * weigh
t_row[x_]);
88. dst_row[x].y += static_cast<short>(src_row[x_].y * weigh
t_row[x_]);
89. dst_row[x].z += static_cast<short>(src_row[x_].z * weigh
t_row[x_]);
90. dst_weight_row[x] += weight_row[x_];
91. }
92. }
93. x_tl /= 2; y_tl /= 2;
94. x_br /= 2; y_br /= 2;
95. }
96. }
97. else// weight_type_ == CV_16S
98. {
100. {
101. for (int y = y_tl; y < y_br; ++y)
102. {
103. int y_ = y - y_tl;
104. const Point3_<short>* src_row = src_pyr_laplace[i].ptr<Poin
t3_<short> >(y_);
105. Point3_<short>* dst_row = dst_pyr_laplace_[i].ptr<Point3_<s
hort> >(y);
106. const short* weight_row = weight_pyr_gauss[i].ptr<short>(y_
);
107. short* dst_weight_row = dst_band_weights_[i].ptr<short>(y);
108.

109. for (int x = x_tl; x < x_br; ++x)
110. {
111. int x_ = x - x_tl;
112. dst_row[x].x += short((src_row[x_].x * weight_row[x_])
>> 8);
113. dst_row[x].y += short((src_row[x_].y * weight_row[x_])
>> 8);
114. dst_row[x].z += short((src_row[x_].z * weight_row[x_])
>> 8);
115. dst_weight_row[x] += weight_row[x_];
116. }
117. }
118. x_tl /= 2; y_tl /= 2;
119. x_br /= 2; y_br /= 2;
120. }
121. }
122. }
其中，金字塔构建的源码：
1. void createLaplacePyr(const Mat &img, int num_levels, vector<Mat> &pyr)
2. {
3. #ifdef HAVE_TEGRA_OPTIMIZATION
4. if(tegra::createLaplacePyr(img, num_levels, pyr))
5. return;
6. #endif
7.
8. pyr.resize(num_levels + 1);
9.
10. if(img.depth() == CV_8U)
11. {
12. if(num_levels == 0)
13. {
14. img.convertTo(pyr[0], CV_16S);
15. return;
16. }
17.
18. Mat downNext;
19. Mat current = img;
20. pyrDown(img, downNext);
21.
22. for(int i = 1; i < num_levels; ++i)
23. {
24. Mat lvl_up;

25. Mat lvl_down;
26.
27. pyrDown(downNext, lvl_down);
28. pyrUp(downNext, lvl_up, current.size());
29. subtract(current, lvl_up, pyr[i-1], noArray(), CV_16S);
30.
31. current = downNext;
32. downNext = lvl_down;
33. }
34.
35. {
36. Mat lvl_up;
37. pyrUp(downNext, lvl_up, current.size());
38. subtract(current, lvl_up, pyr[num_levels-1], noArray(), CV_16S);
39.
40. downNext.convertTo(pyr[num_levels], CV_16S);
41. }
42. }
43. else
44. {
45. pyr[0] = img;
46. //构建高斯金字塔
47. for (int i = 0; i < num_levels; ++i)
48. pyrDown(pyr[i], pyr[i + 1]);//先高斯滤波，在亚采样，得到比 pyr【i】缩
小一半的图像
49. Mat tmp;
50. for (int i = 0; i < num_levels; ++i)
51. {
52. pyrUp(pyr[i + 1], tmp, pyr[i].size());//插值（偶数行，偶数列赋值为 0），
然后高斯滤波，核是 5*5。
53. subtract(pyr[i], tmp, pyr[i]);//pyr[i] = pyr[i]-tmp，得到的全是高频
信息
54. }
55. }
56. }
最终把所有层得金字塔叠加的程序：
1. Mat result, result_mask;
2. blender->blend(result, result_mask);//将多层金字塔图形叠加
源码：

1. void MultiBandBlender::blend(Mat &dst, Mat &dst_mask)
2. {
4. normalizeUsingWeightMap(dst_band_weights_[i], dst_pyr_laplace_[i]);
5.
6. if (can_use_gpu_)
7. restoreImageFromLaplacePyrGpu(dst_pyr_laplace_);
8. else
9. restoreImageFromLaplacePyr(dst_pyr_laplace_);
10.
11. dst_ = dst_pyr_laplace_[0];
12. dst_ = dst_(Range(0, dst_roi_final_.height), Range(0, dst_roi_final_.wid
th));
13. dst_mask_ = dst_band_weights_[0] > WEIGHT_EPS;
14. dst_mask_ = dst_mask_(Range(0, dst_roi_final_.height), Range(0, dst_roi_
final_.width));
15. dst_pyr_laplace_.clear();
16. dst_band_weights_.clear();
17.
18. Blender::blend(dst, dst_mask);
19. }

狠狠撸

Opencv Stitching_detailed algorithm introduction

Recommended

More Related Content

What's hot (20)

Similar to Opencv Stitching_detailed algorithm introduction (20)

Opencv Stitching_detailed algorithm introduction