zoukankan      html  css  js  c++  java

  • 【计算机视觉】stitching_detail算法介绍

    已经不负责图像拼接相关工作，有技术问题请自己解决，谢谢。

    一、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|compressedPlaneA1.5B1|compressedPlanePortraitA2B1

    |compressedPlanePortraitA1.5B1|paniniA2B1|paniniA1.5B1|paniniPortraitA2B1|paniniPortraitA1.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.参数读入

         程序参数读入分析，将程序运行是输入的参数以及需要拼接的图像读入内存，检查图像是否多于2张。

    [cpp] view plain copy
    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 plain copy
    1. work_scale = min(1.0, sqrt(work_megapix * 1e6 / full_img.size().area()));  
    [cpp] view plain copy
    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目前还没用到

    [cpp] view plain copy
    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_aspect = 0.69  
         
        2.3 计算图像特征点，以及计算特征点描述子，并将img_idx设置为i。

    [cpp] view plain copy
    1. (*finder)(img, features[i]);//matcher.cpp 348  
    2. features[i].img_idx = i;  
        特征点描述的结构体定义如下：

    [cpp] view plain copy
    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 plain copy
    1. resize(full_img, img, Size(), seam_scale, seam_scale);  
    2. images[i] = img.clone();  
        3.图像匹配

           对任意两副图形进行特征点匹配，然后使用查并集法，将图片的匹配关系找出，并删除那些不属于同一全景图的图片。

        3.1 使用最近邻和次近邻匹配，对任意两幅图进行特征点匹配。

    [cpp] view plain copy
    1. vector<MatchesInfo> pairwise_matches;//Structure containing information about 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 plain copy
    1.    //Features matcher which finds two best matches for each feature and leaves the best one only if the ratio between descriptor distances is greater than the threshold match_conf.  
    2.   detail::BestOf2NearestMatcher::BestOf2NearestMatcher(bool try_use_gpu=false,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 projective  
    7. transform estimation used in the inliers classification step  
    8. num_matches_thresh2 – Minimum number of matches required for the 2D projective  
    9. transform re-estimation on inliers  
         函数的定义（只是设置一下参数，属于构造函数）：

    [cpp] view plain copy
    1. BestOf2NearestMatcher::BestOf2NearestMatcher(bool try_use_gpu, float match_conf, 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的结构体定义：

    [cpp] view plain copy
    1. Structure containing information about matches between two images. It’s assumed 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 panorama  
    13.     };  
          求出图像匹配的结果如下（具体匹配参见sift特征点匹配），任意两幅图都进行匹配（3*3=9），其中1-》2和2-》1只计算了一次，以1-》2为准，：

         

           3.2 根据任意两幅图匹配的置信度，将所有置信度高于conf_thresh（默认是1.0）的图像合并到一个全集中。

           我们通过函数的参数 save_graph打印匹配结果如下（我稍微修改了一下输出）：

    [cpp] view plain copy
    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 plain copy
    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 plain copy
    1. estimateFocal(features, pairwise_matches, focals);  
           函数的主要源码如下：

    [cpp] view plain copy
    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 lengths 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 plain copy
    1. Tries to estimate focal lengths from the given homography  
    2.     under the assumption that the camera undergoes rotations around its centre 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 plain copy
    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 condition 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：

    [cpp] view plain copy
    1. cameras.assign(num_images, CameraParams());  
    2. for (int i = 0; i < num_images; ++i)  
    3.     cameras[i].focal = focals[i];  

         4.2 根据匹配的内点构建最大生成树，然后广度优先搜索求出根节点，并求出camera的R矩阵，K矩阵以及光轴中心

          camera其他参数：

         aspect = 1.0，ppx，ppy目前等于0，最后会赋值成图像中心点的。

          K矩阵的值：


    [cpp] view plain copy
    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 plain copy
    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. }  
             光轴中心的值

    [cpp] view plain copy
    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 plain copy
    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 plain copy
    1. setUpInitialCameraParams(cameras);  

          函数主要源码：

    [cpp] view plain copy
    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 plain copy
    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的参数（具体理论还没有看懂）

    [cpp] view plain copy
    1. //计算所有内点之和  
    2.     for (size_t i = 0; i < edges_.size(); ++i)  
    3.         total_num_matches_ += static_cast<int>(pairwise_matches[edges_[i].first * num_images_ +  
    4.                                                                 edges_[i].second].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

    [cpp] view plain copy
    1. obtainRefinedCameraParams(cameras);//Gets the refined camera parameters.  
           函数源代码：

    [cpp] view plain copy
    1. void BundleAdjusterRay::obtainRefinedCameraParams(vector<CameraParams> &cameras) 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 plain copy
    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_centers);  
    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 plain copy
    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 horizontal (or vertical).  
    5.        for (size_t i = 0; i < cameras.size(); ++i)  
    6.            cameras[i].R = rmats[i];  
           其中waveCorrect(rmats, wave_correct)源码如下：

    [cpp] view plain copy
    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 eigenvectors 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;  
    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) / getTickFrequency()) << " sec");  
    67. }  

         7.单应性矩阵变换

          由图像匹配，Bundle Adjustment算法以及波形校验，求出了图像的相机参数以及旋转矩阵，接下来就对图形进行单应性矩阵变换，亮度的增量补偿以及多波段融合（图像金字塔）。首先介绍的就是单应性矩阵变换：

          源图像的点（x，y，z=1），图像的旋转矩阵R，图像的相机参数矩阵K，经过变换后的同一坐标（x_,y_,z_），然后映射到球形坐标（u，v，w），他们之间的关系如下：

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

         反投影的公式：

    [cpp] view plain copy
    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中

    [cpp] view plain copy
    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 plain copy
    1. remap(src, dst, xmap, ymap, interp_mode, border_mode);//重映射，xmap，yamp分别是u，v反投影对应的x，y值，默认是双线性插值  
           
          8.光照补偿

          图像拼接中，由于拍摄的图片有可能因为光圈或者光线的问题，导致相邻图片重叠区域出现亮度差，所以在拼接时就需要对图像进行亮度补偿，（opencv只对重叠区域进行了亮度补偿，这样会导致图像融合处虽然光照渐变，但是图像整体的光强没有柔和的过渡）。

          首先，将所有图像，mask矩阵进行分块（大小在32*32附近）。

    [cpp] view plain copy
    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_per_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_tl, bl_br)),  
    19.                                               masks[img_idx].second));  
    20.           }  
    21.       }  
    22.   }  

          然后，求出任意两块图像的重叠区域的平均光强，

    [cpp] view plain copy
    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() - corners[i]));  
    12.                 subimg2 = images[j](Rect(roi.tl() - corners[j], roi.br() - corners[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 == masks[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 plain copy
    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. }  
    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 plain copy
    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"，该方法是利用最大流方法检测。

    [cpp] view plain copy
    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-based 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(GraphCutSeamFinderBase::COST_COLOR);  
    10.         else  
    11. #endif  
    12.             seam_finder = new detail::GraphCutSeamFinder(GraphCutSeamFinderBase::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(GraphCutSeamFinderBase::COST_COLOR_GRAD);  
    19.         else  
    20. #endif  
    21.             seam_finder = new detail::GraphCutSeamFinder(GraphCutSeamFinderBase::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 << "' ";  
    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 plain copy
    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  
    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_scale);  
    31.               }  
    32.   
    33.               Mat K;  
    34.               cameras[i].K().convertTo(K, CV_32F);  
    35.               Rect roi = warper->warpRoi(sz, K, cameras[i].R);//Returns Projected 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_CONSTANT, 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：

    [cpp] view plain copy
    1. //块光照补偿  
    2. void BlocksGainCompensator::apply(int index, Point /*corner*/, Mat &image, const 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 plain copy
    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())) * blend_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*>(static_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->numBands());  
    14.             }  
    15.             else if (blend_type == Blender::FEATHER)  
    16.             {  
    17.                 FeatherBlender* fb = dynamic_cast<FeatherBlender*>(static_cast<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 plain copy
    1. blender->feed(img_warped_s, mask_warped, corners[img_idx]);//将图像写入金字塔中  
          源码：

    [cpp] view plain copy
    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 this 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_bands_);  
    24.     height += ((1 << num_bands_) - height % (1 << num_bands_)) % (1 << num_bands_);  
    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<Point3_<short> >(y_);  
    80.                 Point3_<short>* dst_row = dst_pyr_laplace_[i].ptr<Point3_<short> >(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 * weight_row[x_]);  
    88.                     dst_row[x].y += static_cast<short>(src_row[x_].y * weight_row[x_]);  
    89.                     dst_row[x].z += static_cast<short>(src_row[x_].z * weight_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<Point3_<short> >(y_);  
    105.                 Point3_<short>* dst_row = dst_pyr_laplace_[i].ptr<Point3_<short> >(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. }  
            其中，金字塔构建的源码：
    [cpp] view plain copy
    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. }  
          最终把所有层得金字塔叠加的程序：

    [cpp] view plain copy
    1. Mat result, result_mask;  
    2. blender->blend(result, result_mask);//将多层金字塔图形叠加  
         源码：

    [cpp] view plain copy
    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_.width));  
    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. }  
  • 相关阅读:
    算法----(1)冒泡排序
    淘宝爬虫
    爬虫_豆瓣电影top250 （正则表达式）
    爬虫_猫眼电影top100（正则表达式）
    Android 简单调用摄像头
    Android 简单天气预报
    思维模型
    This view is not constrained, it only has designtime positions, so it will jump to (0,0) unless you
    Android studio preview界面无法预览，报错render problem
    Android studio 3.1.2报错，no target device found
  • 原文地址：https://www.cnblogs.com/huty/p/8516981.html
Copyright © 2011-2022 走看看