【計算機視覺】stitching_detail算法介紹

時間  2019-11-06
標籤 計算機視覺 stitching 算法 介紹 欄目 快樂工作 简体版
原文   原文鏈接

已經不負責圖像拼接相關工做，有技術問題請本身解決，謝謝。css

1、stitching_detail程序運行流程html

      1.命令行調用程序，輸入源圖像以及程序的參數算法

      2.特徵點檢測，判斷是使用surf仍是orb，默認是surf。api

      3.對圖像的特徵點進行匹配，使用最近鄰和次近鄰方法，將兩個最優的匹配的置信度保存下來。app

      4.對圖像進行排序以及將置信度高的圖像保存到同一個集合中，刪除置信度比較低的圖像間的匹配，獲得能正確匹配的圖像序列。這樣將置信度高於門限的全部匹配合併到一個集合中。ide

     5.對全部圖像進行相機參數粗略估計，而後求出旋轉矩陣函數

     6.使用光束平均法進一步精準的估計出旋轉矩陣。ui

     7.波形校訂，水平或者垂直this

     8.拼接spa

     9.融合，多頻段融合，光照補償，


2、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

    命令使用實例，以及程序運行時的提示：


3、程序代碼分析

    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 << "'\n";  
  30.         return 1;  
  31.     }  
      程序解讀可參見：

http://blog.csdn.net/chlele0105/article/details/12624541

http://blog.csdn.net/zouxy09/article/details/8534954

http://blog.csdn.net/yangtrees/article/details/7930640

     

     10.多波段融合

      因爲之前進行處理的圖片都是以work_scale（3.1節有介紹）進行縮放的，因此圖像的內參，corner（同一座標後的頂點），mask（融合的掩碼）都須要從新計算。以及將以前計算的光照加強的gain也要計算進去。

[cpp]  view 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. }  
相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程