Opencv Stitching_detailed algorithm introduction

54 %
46 %
Information about Opencv Stitching_detailed algorithm introduction

Published on January 13, 2016

Author: williamzhang6

Source: slideshare.net

1. 一、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. |compressedPlanePortraitA1.5B1|paniniA2B1|paniniA1.5B1|paniniPortraitA2B1|paniniPor traitA1.5B1|mercator|transverseMercator) 选择融合的平面,默认是球形 --seam_megapix <float> 拼接缝像素的大小 默认是 0.1 ------------ --seam (no|voronoi|gc_color|gc_colorgrad) 拼接缝隙估计方法 默认是 gc_color --compose_megapix <float> 拼接分辨率,默认为-1 --expos_comp (no|gain|gain_blocks) 光照补偿方法,默认是 gain_blocks --blend (no|feather|multiband) 融合方法,默认是多频段融合 --blend_strength <float> 融合强度,0-100.默认是 5. --output <result_img> 输出图像的文件名,默认是 result,jpg 命令使用实例,以及程序运行时的提示: 三、程序代码分析 1.参数读入

3. 程序参数读入分析,将程序运行是输入的参数以及需要拼接的图像读入内存,检查图像 是否多于 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) [cpp] view plaincopy 1. work_scale = min(1.0, sqrt(work_megapix * 1e6 / full_img.size().area())); [cpp] view plaincopy 1. resize(full_img, img, Size(), work_scale, work_scale); 图像大小是 740*830,面积大于 6*10^5,所以计算出 work_scale = 0.98,然后对图像 resize。

4. 2.2 计算 seam_scale,也是根据图像的面积小于 seam_megapix*10^6, (seam_megapix 默认是 0.1),seam_work_aspect 目前还没用到 [cpp] view plaincopy 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。 [cpp] view plaincopy 1. (*finder)(img, features[i]);//matcher.cpp 348 2. features[i].img_idx = i; 特征点描述的结构体定义如下: [cpp] view plaincopy 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[]中 [cpp] view plaincopy

5. 1. resize(full_img, img, Size(), seam_scale, seam_scale); 2. images[i] = img.clone(); 3.图像匹配 对任意两副图形进行特征点匹配,然后使用查并集法,将图片的匹配关系找出,并删 除那些不属于同一全景图的图片。 3.1 使用最近邻和次近邻匹配,对任意两幅图进行特征点匹配。 [cpp] view plaincopy 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 函数: [cpp] view plaincopy 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 函数的定义(只是设置一下参数,属于构造函数): [cpp] view plaincopy 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

6. 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 的结构体定义: [cpp] view plaincopy 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 为准,:

7. 3.2 根据任意两幅图匹配的置信度,将所有置信度高于 conf_thresh(默认是 1.0)的 图像合并到一个全集中。 我们通过函数的参数 save_graph 打印匹配结果如下(我稍微修改了一下输出): [cpp] view plaincopy 1. "outimage101.jpg" -- "outimage102.jpg"[label="Nm=866, Ni=637, C=2.37864"]; 2. "outimage101.jpg" -- "outimage103.jpg"[label="Nm=165, Ni=59, C=1.02609"]; 3. "outimage102.jpg" -- "outimage103.jpg"[label="Nm=460, Ni=260, C=1.78082"]; Nm 代表匹配的数量,NI 代表正确匹配的数量,C 表示置信度 通过查并集方法,查并集介绍请参见博文: http://blog.csdn.net/skeeee/article/details/20471949 [cpp] view plaincopy 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。 [cpp] view plaincopy 1. estimateFocal(features, pairwise_matches, focals); 函数的主要源码如下: [cpp] view plaincopy 1. for (int i = 0; i < num_images; ++i)

8. 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 的定义如下: [cpp] view plaincopy 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. 函数的源码: [cpp] view plaincopy 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. 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: [cpp] view plaincopy 1. cameras.assign(num_images, CameraParams()); 2. for (int i = 0; i < num_images; ++i) 3. cameras[i].focal = focals[i];

10. 4.2 根据匹配的内点构建最大生成树,然后广度优先搜索求出根节点,并求出 camera 的 R 矩阵,K 矩阵以及光轴中心 camera 其他参数: aspect = 1.0,ppx,ppy 目前等于 0,最后会赋值成图像中心点的。 K 矩阵的值: [cpp] view plaincopy 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 矩阵的值: [cpp] view plaincopy 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;

11. 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. } 光轴中心的值 [cpp] view plaincopy 1. for (int i = 0; i < num_images; ++i) 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 ): [cpp] view plaincopy 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_的值: [cpp] view plaincopy 1. setUpInitialCameraParams(cameras);

12. 函数主要源码: [cpp] view plaincopy 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); 16. cam_params_.at<double>(i * 4 + 2, 0) = rvec.at<float>(1, 0); 17. cam_params_.at<double>(i * 4 + 3, 0) = rvec.at<float>(2, 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 删除置信度小于门限的匹配对 [cpp] view plaincopy 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. {

13. 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 的参数(具体理论还没有看懂) [cpp] view plaincopy 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;

14. 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 [cpp] view plaincopy 1. obtainRefinedCameraParams(cameras);//Gets the refined camera parameters. 函数源代码: [cpp] view plaincopy 1. void BundleAdjusterRay::obtainRefinedCameraParams(vector<CameraParams> &came ras) const 2. { 3. for (int i = 0; i < num_images_; ++i) 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); 9. rvec.at<double>(1, 0) = cam_params_.at<double>(i * 4 + 2, 0); 10. rvec.at<double>(2, 0) = cam_params_.at<double>(i * 4 + 3, 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 [cpp] view plaincopy 1. // Normalize motion to center image 2. Graph span_tree;

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

16. 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)源码如下: [cpp] view plaincopy 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); 26. for (size_t i = 0; i < rmats.size(); ++i) 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;

17. 34. if (kind == WAVE_CORRECT_HORIZ) 35. { 36. for (size_t i = 0; i < rmats.size(); ++i) 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. { 46. for (size_t i = 0; i < rmats.size(); ++i) 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); 59. Mat(rg1.t()).copyTo(tmp); 60. tmp = R.row(2); 61. Mat(rg2.t()).copyTo(tmp); 62. 63. for (size_t i = 0; i < rmats.size(); ++i) 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),他们之间的关系如下:

18. [cpp] view plaincopy 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. } 根据映射公式,对图像的上下左右四个边界求映射后的坐标,然后确定变换后图像的 左上角和右上角的坐标, 如果是球形拼接,则需要再加上判断(暂时还没研究透): [cpp] view plaincopy 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. {

19. 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. } 然后利用反投影将图形反投影到变换的图像上,像素计算默认是二维线性插值。 反投影的公式: [cpp] view plaincopy 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.

20. 16. if (z > 0) { x /= z; y /= z; } 17. else x = y = -1; 18. } 然后将反投影求出的 x,y 值写入 Mat 矩阵 xmap 和 ymap 中 [cpp] view plaincopy 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 函数将图像重新绘制: [cpp] view plaincopy 1. remap(src, dst, xmap, ymap, interp_mode, border_mode);//重映射,xmap,yamp 分别 是 u,v 反投影对应的 x,y 值,默认是双线性插值 8.光照补偿 图像拼接中,由于拍摄的图片有可能因为光圈或者光线的问题,导致相邻图片重叠区 域出现亮度差,所以在拼接时就需要对图像进行亮度补偿,(opencv 只对重叠区域进行了 亮度补偿,这样会导致图像融合处虽然光照渐变,但是图像整体的光强没有柔和的过渡)。 首先,将所有图像,mask 矩阵进行分块(大小在 32*32 附近)。 [cpp] view plaincopy 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;

21. 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. } 然后,求出任意两块图像的重叠区域的平均光强, [cpp] view plaincopy 1. //计算每一块区域的光照均值 sqrt(sqrt(R)+sqrt(G)+sqrt(B)) 2. //光照均值是对称矩阵,所以一次循环计算两个光照值,(i,j),与(j,i) 3. for (int i = 0; i < num_images; ++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]));

22. 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. } 建立方程,求出每个光强的调整系数 [cpp] view plaincopy 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) 3. for (int i = 0; i < num_images; ++i) 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. }

23. 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_进行线性滤波 [cpp] view plaincopy 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)

24. 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",该方法 是利用最大流方法检测。 [cpp] view plaincopy 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. { 7. #ifdef HAVE_OPENCV_GPU 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. { 16. #ifdef HAVE_OPENCV_GPU 17. if (try_gpu && gpu::getCudaEnabledDeviceCount() > 0) 18. seam_finder = new detail::GraphCutSeamFinderGpu(GraphCutSeamFind erBase::COST_COLOR_GRAD); 19. else

25. 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 也要计算进去。 [cpp] view plaincopy 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.

26. 17. // Update corners and sizes 18. for (int i = 0; i < num_images; ++i) 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);

27. 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: [cpp] view plaincopy 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 顶点和图像的大小确定最终全景图的尺寸。 [cpp] view plaincopy

28. 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 中。 最后将各层的金字塔叠加,就得到了最终完整的全景图。 [cpp] view plaincopy 1. blender->feed(img_warped_s, mask_warped, corners[img_idx]);//将图像写入金字塔 中 源码:

29. [cpp] view plaincopy 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.

30. 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. {

31. 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. { 99. for (int i = 0; i <= num_bands_; ++i) 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.

32. 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. } 其中,金字塔构建的源码: [cpp] view plaincopy 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;

33. 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. } 最终把所有层得金字塔叠加的程序: [cpp] view plaincopy 1. Mat result, result_mask; 2. blender->blend(result, result_mask);//将多层金字塔图形叠加 源码: [cpp] view plaincopy

34. 1. void MultiBandBlender::blend(Mat &dst, Mat &dst_mask) 2. { 3. for (int i = 0; i <= num_bands_; ++i) 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. }

#ifdef presentations

Add a comment

Related pages

Fast Panorama Stitching - Intel® Software

Introduction Taking panoramic ... OpenCV* and PanoTools* ... We will be using OpenCV’s sample application “cpp-example-stitching_detailed” as a ...
Read more

remove cuda module · opencv/opencv@124ac15 · GitHub

opencv - Open Source ... -CUDA Module Introduction ... 4 samples/cpp/stitching_detailed.cpp @@ -673,7 +673,7 @@ int main(int argc, char* argv[]) seam ...
Read more

stitching opencv example - Findeen

The source code is very easy. I think if we use stitching algorithm in realtime, we should be programing by GPU. ///// #include < stdio.h > #include ...
Read more

opencv: File List - doxygen documentation | Fossies Dox

... OpenCV (Open Source ... algorithm.cpp alloc.cpp arithm ... windows_visual_studio_Opencv introduction_windows_vs.cpp
Read more

[OpenCV] 影像拼接2 (Image Stitching) | 逍遙文工作室

1>stitching_detailed.obj : ... 重裝純的2.4.0後會發生 opencv_core240.lib ... algorithm; coffee; article; sex; line; girl; culture;
Read more

OpenCV - CMakeCache.txt - OpenCV DevZone

//This is the directory where this CMakeCache.txt was created 1012: CMAKE_CACHEFILE_DIR: ... /OpenCV/opencv/modules/core/src/algorithm.cpp;C: ...
Read more

code.opencv.org

code.opencv.org
Read more

opencv_cuda : sort unclassified routines from the module ...

@@ -6,4 +6,5 @@ set(the_description "CUDA-accelerated Computer Vision (legacy)") ocv_warnings_disable(CMAKE_CXX_FLAGS /wd4127 /wd4130 /wd4324 /wd4512 ...
Read more

Panorama Stitching - scribd.com

Hand-picked favorites from our editors. Editors' Picks Sheet Music. Hand-picked favorites from our editors
Read more