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. }  
  • 相关阅读:
    忘记密码破解
    关于本地变量的理解
    MVC的请求过程（或者MVC三者的关系）
    static 静态 关键字
    博客搬家通知
    C#一个可以马上跑起来的反射例子Assembly的使用
    C#之DateTime日期格式解析
    AddressParsing在C#中好用的地址拆分地址结构化库Net5
    C#中获取本地IP地址方法
    ComdeDom生成对象Emit之引用其他成员类库
  • 原文地址：https://www.cnblogs.com/huty/p/8516981.html
Copyright © 2011-2022 走看看