710官方网站以下部分情节。2.4图纸并接….3输入源图像以及程序的参数。

并不曾还念了,不过感觉还是有必要举行一个笔记的,毕竟这仅仅是随笔不是文章,所以可以发小写多少,也终于工作总了,最着重之是这好于足,完成具有有义文档的摸,比打自己之word来说高级很多~~~。

1.概述..2

曾经休担当图像拼接相关工作,有技艺问题要好解决,谢谢。

以下一些内容

2.重要步骤..2

一、stitching_detail程序运行流程

以下

2.1.图像获取….2

      1.命令行调用程序,输入源图像及程序的参数

全文翻译。不过我实在是有限零星当念。。。

2.2鱼类眼图像矫正….2

      2.特征点检测,判断是以surf还是orb,默认是surf。

http://www.docin.com/p-817312337.html

2.3图片匹配….3

     
3.针对图像的特征点进行匹配,使用以来邻和次近邻方法,将少独最美妙的匹配的购信度保存下来。

 

2.4图并接….3

     
4.对图像进行排序和将置信度高之图像保存到和一个集合中,删除置信度比较低之图像中的相当,得到能是匹配的图像序列。这样以贾信度高于门限的具有匹配合并至一个汇聚中。

一、简介:

段同样:我们的拼凑不欲手动校直。

段二:根据文献,图像自动对那和拼接有少数学方式:1、直接的 2、基于特征的

直接的主意:1、提供十分精确之一定,2、但是需要初始化处理。

冲特征的配准:1、不需初始化,2、但是可靠性上相差。(有或是Harris)

段三:

我们的配合有:

1、可靠匹配,在转悠,缩放和光照还无安宁的场面下。

2、给起无关图,也能找到有关图里的关系。【

来东西需要实现转了,刚刚的拼接,并没控制变量。

再试一次。

委发现了autostitch具有双重强的可靠性,stitch会跑起大,不安静。不亮 stitching_detail具有多胜过的康乐。

因此我们得先stitch 然后重新 detail 最后再 autostitch?

3、多波段融合,结果越可靠。

4、增益补偿

5、自动校直

6、捆绑调整, 展示了对于多独波段的重合图像进行多波段融合。

 

段四:

其余部分结构如下:

亚有些说明,所研究问题的几乎何学,以及选择未转移特征的原委。

老三片介绍了:图像匹配方法(RANSAC)和说明图像匹配的概率模型。

季部分:描述图像对准算法(捆绑调整),共同优化摄像头的参数。

五交七有:处理过程 自动校直,增益补偿,多波段融合

【八片去哪了?】

九片段:结论与未来展望

 

2.5图像融合….3

     5.对准所有图像进行相机参数粗略估算,然后要出旋转矩阵

其次、特征匹配:

段一:

全景识别算法的首先步是当具有图像里提取及匹配SITF特征点。

Sift特征监测子位于不同尺度空间高斯差值函数的极值点处。

Sift特征点使得特征尺度和可行性为确定。为测量提供了一个般不转换的组织。对于仿射变化有所一定之健壮性。【什么是放射变换?】

空中累积计算实现活动不变性。

亮度不变性由梯度(消除偏差)

描述子矢量归一化(消除增益)

段二:sift在转悠和原则变化时凡休移的,因此可处理缩放图像,旋转图像。

Harris角点图像修补的相关性特征匹配技术就是无可知落实。

因此我是未是也只要跃跃欲试一下旋转图像方向。

风土人情技术实现的相关性在图像旋转时候是变化之。

Harris不仅于转悠时候别,在标准化变换时,也是浮动之。

段三:

若果相机绕光学中心旋转,图像的变化群是一个应和矩阵的新鲜群。由一个旋转矢量710官方网站 1与焦距将每个摄像头参数化,就深受有了成对的附和矩阵,710官方网站 2,其中

710官方网站 3(1)

并且710官方网站 4710官方网站 5大凡咸匀的图像左边(710官方网站 6,其中710官方网站 7凡是二维的图像坐标)。4参数相机模型定义也:

710官方网站 8(2)

 

对旋转使用指数表示:

710官方网站 9710官方网站 10(3)

 

段四:

 

每当是变换群中,理想之规格下用会见使无转换的图像特点。可是,在图像左边被对此小的更换表示如下:

 

710官方网站 11(4)

要么当于710官方网站 12其中,

710官方网站 13(5)

举凡通过一个关于710官方网站 14的对应线性化得到的仿射变化。这象征每个微之图像修补经过逐一仿射变换,并且成立采取了当仿射变换下部分无转移的SIFT特征。

段五:

只要从持有n个图像中取特征点后(线性时间外),需对特征点进行匹配。由于大多独图像可能重叠在一个单一的光明上,在特点空间内每个特征点需和她最近的k个临域点配合(k=4),通过使用k-d树算法找到类似最近底临域点,时间复杂度为O(nlogn)。K-d树是一致种植到对那的二进制空间划分,它当平均高方差维递归划分特征空间。

 

2.6全景图像投射….3

     6.使用光束平均法更加精准的估算起旋转矩阵。

老三、图像匹配

段同样:图像匹配的对象是找到有匹配(例如重叠)图像,然后图像匹配联通集会成为全景图。由于每个图或和擅自其他一个配合,这个题目同样开始就是呈现出是图像数的2次方。【所以这个数字应该是几?对于我们曾输入的n张图片来讲。】为了取一个好之拼接结果,对于图像而言,每个图像只待同个别叠的图像来配合。【如果重叠那拼接结果从是最最好】

段同样总结:问题难度巨大,要是发生重合就哼惩治。

段二:从特征匹配步骤中,已经查找有大量配合的图像中有大量相当配点的图像。对于当前图像,我们将m幅图形作为或的相当图像(假设m=6),这m幅图像和时图像发极度充分数目之特点匹配点。首先,使用RANSAC算法选择一样多级及图像中对许矩阵兼容的内点,然后使用概率模型做更的求证。

段二总:用RANSAC算法。

【备选图像发m张,那么到底要无若入选,需要RANSAC算法来控制】

3.算法技巧点介绍..3

     7.波形校正,水平还是垂直

3.1运用RANSAC算法的矫健对承诺矩阵估计

段三:RANSAC(随机取样一致性算法)【RANSAC算法是使用最少的同等组自由采样匹配点的平等栽健康估计过程,这个看下原文那立句话是纪念说,这种算法采用的匹配点少,并且为死年轻力壮?】用来估计图像变换参数,并找到与数量有极其好一致性的化解方案。【使用RANSAC算法来估算图像变换参数,来找到最佳解决方案。】在全景图的情形下,选择r=4对匹配特征点,使用直接线性变换(DLT)方法计算图像里的应和矩阵H。【r是什么?R代表匹配特征点,此处数量为4.】使用直接线性变换方法算图像中的附和矩阵H。【什么事一直线性变换方法?图像里的相应矩阵又是凭什么?】重复试验500破,选择内点数量极其老之解决方案(在像从误差范围外,其预测与H是均等的)。假而同一对郎才女貌图像里的特征匹配点的正确率(内点概率)为,n次试验后找到正确变换概率为:

710官方网站 15 (6)

 

由此大量试后,找到正确对应矩阵的几率非常特别。例如,对于内点概率710官方网站 16=0.5,在500涂鸦试后,未找到科学对应矩阵的票房价值为710官方网站 17。【那这地方我们是免是可稍微弄得稍微一些。以加强运行的效率,抑或者因为某些用途,可以延续于这个数字转换多少。不过肯定都远非必要了,又或者我们可调动是比例参数,在一些重合区域的图像内点概率就会坏十分,那么实验次数虽得对应的稍,来弱化多少计算压力,对于某些内点概率小,地方我们可以要求尝试次数尽量的那个,所以啊是内接触吧?】

 

段三总结:当内点采用某平等数值的事态下,我们可以通过增加试行次数来加强相应矩阵的配合程度。

 

段四:RANSAC算法本质上是同等栽估计H的采样方法,如果因此对数似然和的最大化代替内点数量之最大化,结果是最最要命似然估计(MLE)。此外,如果换参数的先验值是行之有效之,可以算起最好深后验概率(MAP)。这些算法为分级成为MLESAC和MAPSAC。

 

段四总:RANSAC本质上是一样种量图像中对承诺矩阵的采样方法。如果因此对数似然和的最大化代替内点数量之最大化,结果是最好充分似然估计(MLE)。如果参数先验值有效,可以算起尽可怜后验概率MAP。对应之算法分别成为:MLESAC 和MAPSAC。不过为何对数似然和之最大化为什么可以代表内点数量之最大化。内点指得是呀?

 

3.1图像获取….3

     8.拼接

3.2图像匹配关系验证的概率模型

 

段五:对少少贪图如中是否在相当关系,我们利用同一雨后春笋几哪里一样的特色匹配点(RANSAC内点)和一致系列以重叠区域外,但切莫一样的特征点(RANSAC外点)来说明。【我接近想起来了,RANSAC貌似是同栽算法,界外的就是外点,属于无过关的触及,不相互配合的触发,而RANSAC内接触,就是特点匹配点。】验证模型通过比较这些对匹配有的一致名目繁多内点和左匹配有的等同系列外点的概率来拓展验证。

 

段五总结:是否匹配 通过
特征一致点和特色未一致点来判定。那下面将说判定方法了:

 

段六:对于同适合给定的图像,重叠区域外总的匹配特征点数量为710官方网站 18,内点数也710官方网站 19。图像是否可行匹配通过二进制变量710官方网站 20代表。第i个匹配特征点710官方网站 21是不是为内点被如为独立的伯努利分布,以便让内接触总数从伯努利分部:

 

 

 

710官方网站 22(7)

710官方网站 23(8)

 

其中,710官方网站 24举凡可正确匹配图像时特征点内点的几率,710官方网站 25凡勿可知兑现图像匹配时特征点为内点的概率;710官方网站 26

代表特征匹配点变量的会师710官方网站 27,内点数710官方网站 28710官方网站 29举凡伯努利分部,表示如下:710官方网站 30.(9)

 

选择710官方网站 31=0.6,710官方网站 32=0.1,则好下贝叶斯规则(式10、11)计算科学图像匹配的先验概率。

710官方网站 33(10)

       =710官方网站 34(11)

假设满足710官方网站 35

 

 

710官方网站 36(12)

咱俩可兑现图像拼接。假定710官方网站 37710官方网站 38,进一步得出正确图像匹配的判断条件:

710官方网站 39 (13)

其中,710官方网站 40=8.0,710官方网站 41=0.3。尽管当马上我们捎了710官方网站 42710官方网站 43710官方网站 44710官方网站 45710官方网站 46凡是咱若的值,这些数据我们得更改的,或者
是可以给他更加纯粹的。比如,我们可由此以那个之数据集中计算一部分匹配点和是的相应矩阵相平等,来打量710官方网站 47。最终独自是为着总结这样一个公式710官方网站 48,要想然匹配,必须满足内点数大于这么一个同总匹配特征点的一个函数。如果低于等于,就不可知实现对匹配。其中710官方网站 49。(内点数要低于总点数。)

段七:

倘若图像里的匹配点对
确定,我么你得找到全景序列作为连接匹配图像集,它可辨认图像集中之大都个全景,拒绝
不般配的噪音 图像。(见图2)

710官方网站 50

 

710官方网站 51

(g)根据对应矩阵的图像对准

希冀1,从具有图像中提取SIFT特征点。使用k-d数匹有特征点后,对于一个加图像,用最为多特点匹配点的m幅图像进行图像匹配。首先执行RANSAC算法计算起对承诺矩阵,然后采用概率模型验证基于内点数的图像匹配,在此例子中,输入图像是517*374,有247只不利特征匹配点。

710官方网站 52

710官方网站 53

图2,可识别全景图。考虑一个特色匹配点的噪声集,我们运用RANSAC算法好概率验证过程找到同样的图像匹配(a),每个图像对 间的 箭头表示在图像对
间找到同样的风味匹配点集,图像匹配连接分量被找到(b),拼接成全景图(c);注意到拖欠算法对 不属全景图的噪音图像
不灵活。

 

 

扎调整下一个号在拍卖。。。

 

3.2鲜鱼眼图像矫正….4

     9.融合,多频段融合,光照补偿,

3.3图形匹配….5

二、stitching_detail程序接口介绍

3.3.1和特性无关之相当方式….5

    img1 img2 img3 输入图像

3.3.2基于特征进行匹配的法门….5

     –preview
 以预览模式运行程序,比正规模式一旦尽快,但输出图像分辨率低,拼接的分辨率compose_megapix 设置为0.6

3.4图纸并接….6

     –try_gpu  (yes|no)  是否采用GPU(图形处理器),默认为no

3.5图像融合….7

/* 运动估计参数 */

3.5.1平均叠加法….7

    –work_megapix <–work_megapix <float>>
图像匹配的分辨率大小,图像的面积尺寸变为work_megapix*100000,默认为0.6

3.5.2线性法….7

    –features (surf|orb) 选择surf或者orb算法进行特征点计算,默认为surf

3.5.3加权函数法….7

    –match_conf <float>
特征点检测置信等级,最近邻匹配距离与不好走近邻匹配距离的比率,surf默认为0.65,orb默认为0.3

3.5.4基本上段融合法(多分辨率样条)….7

    –conf_thresh <float> 两轴图来源同全景图的进信度,默认为1.0

3.6全景图像投射….7

    –ba (reproj|ray) 光束平均法的误差函数选择,默认是ray方法

3.6.1柱面全景图….7

    –ba_refine_mask (mask)  —————

3.6.2球面全景图….8

    –wave_correct (no|horiz|vert) 波形校验 水平,垂直或者无
默认是horiz

3.6.3多面体全景图….8

     –save_graph <file_name>
将匹配的图形以点的花样保留至文件被,
Nm代表匹配的数码,NI代表对匹配的数额,C表示请信度

4.开头源图像算法库OpenCV拼接模块..9

/*图像融合参数:*/

4.1

    –warp
(plane|cylindrical|spherical|fisheye|stereographic|compressedPlaneA2B1|compressedPlaneA1.5B1|compressedPlanePortraitA2B1

stitching_detail程序运行流程….9

|compressedPlanePortraitA1.5B1|paniniA2B1|paniniA1.5B1|paniniPortraitA2B1|paniniPortraitA1.5B1|mercator|transverseMercator)

4.2

    选择融合的面,默认是球形

stitching_detail程序接口介绍….9

    –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

    命令下实例,以及程序运行时的提醒:

710官方网站 54

其三、程序代码分析

    1.参数读入

   
 程序参数读入分析,将程序运行是输入的参数和用拼接的图像读入内存,检查图像是否多于2张。

    int retval = parseCmdArgs(argc, argv);
    if (retval)
        return retval;

    // Check if have enough images
    int num_images = static_cast<int>(img_names.size());
    if (num_images < 2)
    {
        LOGLN("Need more images");
        return -1;
    }

      2.特征点检测

     
判断选择是surf还是orb特征点检测(默认是surf)以及对图像进行事先处理(尺寸缩放),然后计算各国幅图的特征点,以及特征点描述子

      2.1
计算work_scale,将图像resize到面积在work_megapix*10^6以下,(work_megapix
默认是0.6)

work_scale = min(1.0, sqrt(work_megapix * 1e6 / full_img.size().area()));

resize(full_img, img, Size(), work_scale, work_scale);

      图像大小是740*830,面积过6*10^5,所以算出work_scale =
0.98,然后对图像resize。 

     710官方网站 55

     2.2
计算seam_scale,也是依据图像的面积小于seam_megapix*10^6,(seam_megapix
默认是0.1),seam_work_aspect目前还未曾因此到

   seam_scale = min(1.0, sqrt(seam_megapix * 1e6 / full_img.size().area()));
   seam_work_aspect = seam_scale / work_scale; //seam_megapix = 0.1 seam_work_aspect = 0.69

     710官方网站 56
    2.3 计算图像特征点,以及计算特征点描述子,并将img_idx设置为i。

  (*finder)(img, features[i]);//matcher.cpp 348
  features[i].img_idx = i;

    特征点描述的组织体定义如下:

   struct detail::ImageFeatures
    Structure containing image keypoints and descriptors.
    struct CV_EXPORTS ImageFeatures
    {
        int img_idx;// 
        Size img_size;
        std::vector<KeyPoint> keypoints;
        Mat descriptors;
    };

    710官方网站 57
     2.4 将自图像resize到seam_megapix*10^6,并存入image[]中

        resize(full_img, img, Size(), seam_scale, seam_scale);
        images[i] = img.clone();

    3.图像匹配

     
 对擅自两合图形进行特征点匹配,然后用查并集法,将图片的配合关系找来,并去那些未属同一全景图的图。

    3.1 使用以来邻和次近邻匹配,对随意两轴图进行特征点匹配。

    vector<MatchesInfo> pairwise_matches;//Structure containing information about matches between two images. 
    BestOf2NearestMatcher matcher(try_gpu, match_conf);//最近邻和次近邻法
    matcher(features, pairwise_matches);//对每两个图片进行matcher,20-》400 matchers.cpp 502

    介绍一下BestOf2NearestMatcher 函数:

      //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.
     detail::BestOf2NearestMatcher::BestOf2NearestMatcher(bool try_use_gpu=false,float match_conf=0.3f,
                                                         intnum_matches_thresh1=6, int num_matches_thresh2=6)
     Parameters:    try_use_gpu – Should try to use GPU or not
            match_conf – Match distances ration threshold
            num_matches_thresh1 – Minimum number of matches required for the 2D projective
            transform estimation used in the inliers classification step
            num_matches_thresh2 – Minimum number of matches required for the 2D projective
            transform re-estimation on inliers

     函数的定义(只是设置一下参数,属于构造函数):

BestOf2NearestMatcher::BestOf2NearestMatcher(bool try_use_gpu, float match_conf, int num_matches_thresh1, int num_matches_thresh2)
{
#ifdef HAVE_OPENCV_GPU
    if (try_use_gpu && getCudaEnabledDeviceCount() > 0)
        impl_ = new GpuMatcher(match_conf);
    else
#else
    (void)try_use_gpu;
#endif
        impl_ = new CpuMatcher(match_conf);

    is_thread_safe_ = impl_->isThreadSafe();
    num_matches_thresh1_ = num_matches_thresh1;
    num_matches_thresh2_ = num_matches_thresh2;
}

     以及MatchesInfo的组织体定义:

Structure containing information about matches between two images. It’s assumed that there is a homography between those images.
    struct CV_EXPORTS MatchesInfo
    {
        MatchesInfo();
        MatchesInfo(const MatchesInfo &other);
        const MatchesInfo& operator =(const MatchesInfo &other);
        int src_img_idx, dst_img_idx; // Images indices (optional)
        std::vector<DMatch> matches;
        std::vector<uchar> inliers_mask; // Geometrically consistent matches mask
        int num_inliers; // Number of geometrically consistent matches
        Mat H; // Estimated homography
        double confidence; // Confidence two images are from the same panorama
    };

     
求出图像匹配的结果如下(具体配合参见sift特征点匹配),任意两帧图都开展匹配(3*3=9),其中1-》2跟2-》1就算了一致糟,以1-》2为以,:

     710官方网站 58

       3.2
根据随便两轴图匹配的买入信度,将有着置信度高于conf_thresh(默认是1.0)的图像合并及一个全都集中。

       我们通过函数的参数
save_graph打印匹配结果如下(我聊修改了一下出口):

"outimage101.jpg" -- "outimage102.jpg"[label="Nm=866, Ni=637, C=2.37864"];
"outimage101.jpg" -- "outimage103.jpg"[label="Nm=165, Ni=59, C=1.02609"];
"outimage102.jpg" -- "outimage103.jpg"[label="Nm=460, Ni=260, C=1.78082"];

      Nm代表匹配的数,NI代表是匹配的数目,C表示请信度

      通过查并集方法,查并集介绍请参见博文:

      http://blog.csdn.net/skeeee/article/details/20471949

    vector<int> indices = leaveBiggestComponent(features, pairwise_matches, conf_thresh);//将置信度高于门限的所有匹配合并到一个集合中
    vector<Mat> img_subset;
    vector<string> img_names_subset;
    vector<Size> full_img_sizes_subset;
    for (size_t i = 0; i < indices.size(); ++i)
    {
        img_names_subset.push_back(img_names[indices[i]]);
        img_subset.push_back(images[indices[i]]);
        full_img_sizes_subset.push_back(full_img_sizes[indices[i]]);
    }

    images = img_subset;
    img_names = img_names_subset;
    full_img_sizes = full_img_sizes_subset;

       4.基于单应性矩阵粗略估算有相机参数(焦距)

        4.1 焦距参数的量

       
根据前求来的随意两幅图的配合,我们根据两帧图的单应性矩阵H,求出符合条件的f,(4副图,16独相当,求来8独符合条件的f),然后要出f的均值或者中值当成所有图形的大概估算的f。

estimateFocal(features, pairwise_matches, focals);

       函数的严重性源码如下:

    for (int i = 0; i < num_images; ++i)
    {
        for (int j = 0; j < num_images; ++j)
        {
            const MatchesInfo &m = pairwise_matches[i*num_images + j];
            if (m.H.empty())
                continue;
            double f0, f1;
            bool f0ok, f1ok;
            focalsFromHomography(m.H, f0, f1, f0ok, f1ok);//Tries to estimate focal lengths from the given homography  79
            //under the assumption that the camera undergoes rotations around its centre only.
            if (f0ok && f1ok)
                all_focals.push_back(sqrt(f0 * f1));
        }
    }

      其中函数focalsFromHomography的概念如下:

Tries to estimate focal lengths from the given homography
    under the assumption that the camera undergoes rotations around its centre only.  
    Parameters
    H – Homography.
    f0 – Estimated focal length along X axis.
    f1 – Estimated focal length along Y axis.
    f0_ok – True, if f0 was estimated successfully, false otherwise.
    f1_ok – True, if f1 was estimated successfully, false otherwise.

     函数的源码:

void focalsFromHomography(const Mat& H, double &f0, double &f1, bool &f0_ok, bool &f1_ok)
{
    CV_Assert(H.type() == CV_64F && H.size() == Size(3, 3));//Checks a condition at runtime and throws exception if it fails

    const double* h = reinterpret_cast<const double*>(H.data);//强制类型转换

    double d1, d2; // Denominators
    double v1, v2; // Focal squares value candidates

    //具体的计算过程有点看不懂啊
    f1_ok = true;
    d1 = h[6] * h[7];
    d2 = (h[7] - h[6]) * (h[7] + h[6]);
    v1 = -(h[0] * h[1] + h[3] * h[4]) / d1;
    v2 = (h[0] * h[0] + h[3] * h[3] - h[1] * h[1] - h[4] * h[4]) / d2;
    if (v1 < v2) std::swap(v1, v2);
    if (v1 > 0 && v2 > 0) f1 = sqrt(std::abs(d1) > std::abs(d2) ? v1 : v2);
    else if (v1 > 0) f1 = sqrt(v1);
    else f1_ok = false;

    f0_ok = true;
    d1 = h[0] * h[3] + h[1] * h[4];
    d2 = h[0] * h[0] + h[1] * h[1] - h[3] * h[3] - h[4] * h[4];
    v1 = -h[2] * h[5] / d1;
    v2 = (h[5] * h[5] - h[2] * h[2]) / d2;
    if (v1 < v2) std::swap(v1, v2);
    if (v1 > 0 && v2 > 0) f0 = sqrt(std::abs(d1) > std::abs(d2) ? v1 : v2);
    else if (v1 > 0) f0 = sqrt(v1);
    else f0_ok = false;
}

      求出之焦距有8只

     710官方网站 59

      求出底焦距取中值或者平均值,然后就是有着图片的焦距。

      并构建camera参数,将矩阵写副camera:

        cameras.assign(num_images, CameraParams());
        for (int i = 0; i < num_images; ++i)
            cameras[i].focal = focals[i];

     4.2
根据匹配的内点构建最充分生成树,然后广度优先找求出干净节点,并告来camera的R矩阵,K矩阵以及光轴中心

      camera其他参数:

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

      K矩阵的价值:

710官方网站 60

Mat CameraParams::K() const
{
    Mat_<double> k = Mat::eye(3, 3, CV_64F);
    k(0,0) = focal; k(0,2) = ppx;
    k(1,1) = focal * aspect; k(1,2) = ppy;
    return k;
}

      R矩阵的价:

     710官方网站 61

    void operator ()(const GraphEdge &edge)
    {
        int pair_idx = edge.from * num_images + edge.to;

        Mat_<double> K_from = Mat::eye(3, 3, CV_64F);
        K_from(0,0) = cameras[edge.from].focal;
        K_from(1,1) = cameras[edge.from].focal * cameras[edge.from].aspect;
        K_from(0,2) = cameras[edge.from].ppx;
        K_from(0,2) = cameras[edge.from].ppy;

        Mat_<double> K_to = Mat::eye(3, 3, CV_64F);
        K_to(0,0) = cameras[edge.to].focal;
        K_to(1,1) = cameras[edge.to].focal * cameras[edge.to].aspect;
        K_to(0,2) = cameras[edge.to].ppx;
        K_to(0,2) = cameras[edge.to].ppy;

        Mat R = K_from.inv() * pairwise_matches[pair_idx].H.inv() * K_to;
        cameras[edge.to].R = cameras[edge.from].R * R;
    }

         光轴中心的值

    for (int i = 0; i < num_images; ++i)
    {
        cameras[i].ppx += 0.5 * features[i].img_size.width;
        cameras[i].ppy += 0.5 * features[i].img_size.height;
    }

      5.采用Bundle Adjustment方法对持有图片展开相机参数校正

      Bundle
Adjustment(光束法平差)算法主要是解决有相机参数的一路。这是全景拼接必须的如出一辙步,因为差不多独成对的单应性矩阵合成全景图时,会忽略全局的限制,造成累积误差。因此各一个图像都要添加光束法平差值,使图像为初始化成相同的转动和焦距长度。

     
光束法平差的靶子函数是一个装有鲁棒性的照射误差的平方和函数。即各级一个特征点都设投到外的图像被,计算产生而误差的平方和极端小的照相机参数。具体的演绎过程可参见Automatic
Panoramic Image Stitching using Invariant
Features.pdf的第五节,本文仅介绍一下opencv实现之历程(完整的驳斥以及公式
暂时还从未看明白,希望有人能一起交流)

   
 opencv中误差描述函数有有限种植如下:(opencv默认是BundleAdjusterRay ):

   BundleAdjusterReproj // BundleAdjusterBase(7, 2)//最小投影误差平方和
    Implementation of the camera parameters refinement algorithm which minimizes sum of the reprojection error squares

    BundleAdjusterRay //  BundleAdjusterBase(4, 3)//最小特征点与相机中心点的距离和
    Implementation of the camera parameters refinement algorithm which minimizes sum of the distances between the
    rays passing through the camera center and a feature.

      5.1 首先计算cam_params_的值:

    setUpInitialCameraParams(cameras);

      函数要源码:

    cam_params_.create(num_images_ * 4, 1, CV_64F);
    SVD svd;//奇异值分解
    for (int i = 0; i < num_images_; ++i)
    {
        cam_params_.at<double>(i * 4, 0) = cameras[i].focal;

        svd(cameras[i].R, SVD::FULL_UV);
        Mat R = svd.u * svd.vt;
        if (determinant(R) < 0)
            R *= -1;

        Mat rvec;
        Rodrigues(R, rvec);
        CV_Assert(rvec.type() == CV_32F);
        cam_params_.at<double>(i * 4 + 1, 0) = rvec.at<float>(0, 0);
        cam_params_.at<double>(i * 4 + 2, 0) = rvec.at<float>(1, 0);
        cam_params_.at<double>(i * 4 + 3, 0) = rvec.at<float>(2, 0);
    }

     
计算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 删除置信度小于门限的相当对

    // Leave only consistent image pairs
    edges_.clear();
    for (int i = 0; i < num_images_ - 1; ++i)
    {
        for (int j = i + 1; j < num_images_; ++j)
        {
            const MatchesInfo& matches_info = pairwise_matches_[i * num_images_ + j];
            if (matches_info.confidence > conf_thresh_)
                edges_.push_back(make_pair(i, j));
        }
    }

       5.3 使用LM算法计算camera参数。

       首先初始化LM的参数(具体理论还尚无看懂)

//计算所有内点之和
    for (size_t i = 0; i < edges_.size(); ++i)
        total_num_matches_ += static_cast<int>(pairwise_matches[edges_[i].first * num_images_ +
                                                                edges_[i].second].num_inliers);

    CvLevMarq solver(num_images_ * num_params_per_cam_,
                     total_num_matches_ * num_errs_per_measurement_,
                     term_criteria_);

    Mat err, jac;
    CvMat matParams = cam_params_;
    cvCopy(&matParams, solver.param);

    int iter = 0;
    for(;;)//类似于while(1),但是比while(1)效率高
    {
        const CvMat* _param = 0;
        CvMat* _jac = 0;
        CvMat* _err = 0;

        bool proceed = solver.update(_param, _jac, _err);

        cvCopy(_param, &matParams);

        if (!proceed || !_err)
            break;

        if (_jac)
        {
            calcJacobian(jac); //构造雅阁比行列式
            CvMat tmp = jac;
            cvCopy(&tmp, _jac);
        }

        if (_err)
        {
            calcError(err);//计算err
            LOG_CHAT(".");
            iter++;
            CvMat tmp = err;
            cvCopy(&tmp, _err);
        }
    }

       计算camera

 obtainRefinedCameraParams(cameras);//Gets the refined camera parameters.

       函数源代码:

void BundleAdjusterRay::obtainRefinedCameraParams(vector<CameraParams> &cameras) const
{
    for (int i = 0; i < num_images_; ++i)
    {
        cameras[i].focal = cam_params_.at<double>(i * 4, 0);

        Mat rvec(3, 1, CV_64F);
        rvec.at<double>(0, 0) = cam_params_.at<double>(i * 4 + 1, 0);
        rvec.at<double>(1, 0) = cam_params_.at<double>(i * 4 + 2, 0);
        rvec.at<double>(2, 0) = cam_params_.at<double>(i * 4 + 3, 0);
        Rodrigues(rvec, cameras[i].R);

        Mat tmp;
        cameras[i].R.convertTo(tmp, CV_32F);
        cameras[i].R = tmp;
    }
}

       求出清节点,然后由一化为旋转矩阵R

    // Normalize motion to center image
    Graph span_tree;
    vector<int> span_tree_centers;
    findMaxSpanningTree(num_images_, pairwise_matches, span_tree, span_tree_centers);
    Mat R_inv = cameras[span_tree_centers[0]].R.inv();
    for (int i = 0; i < num_images_; ++i)
        cameras[i].R = R_inv * cameras[i].R;

     6 波形校正

   
 前面几节约将相机旋转矩阵计算出来,但是还有一个元素需要考虑,就是由于拍摄者拍摄图片的上不自然是水平的,轻微的斜会造成全景图像出现飞机曲线,因此我们如果针对图像进行波形校正,主要是找每幅图的“上升向量”(up_vector),使他校正成710官方网站 62

710官方网站 63

浪校正的效果图

         opencv实现的源码如下(也是临时无看明白,囧)

   vector<Mat> rmats;
        for (size_t i = 0; i < cameras.size(); ++i)
            rmats.push_back(cameras[i].R);
        waveCorrect(rmats, wave_correct);//Tries to make panorama more horizontal (or vertical).
        for (size_t i = 0; i < cameras.size(); ++i)
            cameras[i].R = rmats[i];

       其中waveCorrect(rmats, wave_correct)源码如下:

void waveCorrect(vector<Mat> &rmats, WaveCorrectKind kind)
{
    LOGLN("Wave correcting...");
#if ENABLE_LOG
    int64 t = getTickCount();
#endif

    Mat moment = Mat::zeros(3, 3, CV_32F);
    for (size_t i = 0; i < rmats.size(); ++i)
    {
        Mat col = rmats[i].col(0);
        moment += col * col.t();//相机R矩阵第一列转置相乘然后相加
    }
    Mat eigen_vals, eigen_vecs;
    eigen(moment, eigen_vals, eigen_vecs);//Calculates eigenvalues and eigenvectors of a symmetric matrix.

    Mat rg1;
    if (kind == WAVE_CORRECT_HORIZ)
        rg1 = eigen_vecs.row(2).t();//如果是水平校正,去特征向量的第三行
    else if (kind == WAVE_CORRECT_VERT)
        rg1 = eigen_vecs.row(0).t();//如果是垂直校正,特征向量第一行
    else
        CV_Error(CV_StsBadArg, "unsupported kind of wave correction");

    Mat img_k = Mat::zeros(3, 1, CV_32F);
    for (size_t i = 0; i < rmats.size(); ++i)
        img_k += rmats[i].col(2);//R函数第3列相加
    Mat rg0 = rg1.cross(img_k);//rg1与img_k向量积
    rg0 /= norm(rg0);//归一化?

    Mat rg2 = rg0.cross(rg1);

    double conf = 0;
    if (kind == WAVE_CORRECT_HORIZ)
    {
        for (size_t i = 0; i < rmats.size(); ++i)
            conf += rg0.dot(rmats[i].col(0));//Computes a dot-product of two vectors.数量积
        if (conf < 0)
        {
            rg0 *= -1;
            rg1 *= -1;
        }
    }
    else if (kind == WAVE_CORRECT_VERT)
    {
        for (size_t i = 0; i < rmats.size(); ++i)
            conf -= rg1.dot(rmats[i].col(0));
        if (conf < 0)
        {
            rg0 *= -1;
            rg1 *= -1;
        }
    }

    Mat R = Mat::zeros(3, 3, CV_32F);
    Mat tmp = R.row(0);
    Mat(rg0.t()).copyTo(tmp);
    tmp = R.row(1);
    Mat(rg1.t()).copyTo(tmp);
    tmp = R.row(2);
    Mat(rg2.t()).copyTo(tmp);

    for (size_t i = 0; i < rmats.size(); ++i)
        rmats[i] = R * rmats[i];

    LOGLN("Wave correcting, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");
}

     7.单应性矩阵变换

      由图像匹配,Bundle
Adjustment算法以及波形校验,求来了图像的相机参数与旋转矩阵,接下去便针对图片进行单应性矩阵变换,亮度的增量上及多波段融合(图像金字塔)。首先介绍的就是单应性矩阵变换:

     
源图像的点(x,y,z=1),图像的旋矩阵R,图像的相机参数矩阵K,经过变换后底相同坐标(x_,y_,z_),然后照到球形坐标(u,v,w),他们中的关联如下:

void SphericalProjector::mapForward(float x, float y, float &u, float &v)
{
    float x_ = r_kinv[0] * x + r_kinv[1] * y + r_kinv[2];
    float y_ = r_kinv[3] * x + r_kinv[4] * y + r_kinv[5];
    float z_ = r_kinv[6] * x + r_kinv[7] * y + r_kinv[8];

    u = scale * atan2f(x_, z_);
    float w = y_ / sqrtf(x_ * x_ + y_ * y_ + z_ * z_);
    v = scale * (static_cast<float>(CV_PI) - acosf(w == w ? w : 0));
}

    710官方网站 64

     
 根据映射公式,对图像的光景左右季独境界求映射后底坐标,然后确定变换后图像的左上角与右上较量的坐标,

       如果是球形拼接,则用再次增长判断(暂时还不曾研究透):

    float tl_uf = static_cast<float>(dst_tl.x);
    float tl_vf = static_cast<float>(dst_tl.y);
    float br_uf = static_cast<float>(dst_br.x);
    float br_vf = static_cast<float>(dst_br.y);

    float x = projector_.rinv[1];
    float y = projector_.rinv[4];
    float z = projector_.rinv[7];
    if (y > 0.f)
    {
        float x_ = (projector_.k[0] * x + projector_.k[1] * y) / z + projector_.k[2];
        float y_ = projector_.k[4] * y / z + projector_.k[5];
        if (x_ > 0.f && x_ < src_size.width && y_ > 0.f && y_ < src_size.height)
        {
            tl_uf = min(tl_uf, 0.f); tl_vf = min(tl_vf, static_cast<float>(CV_PI * projector_.scale));
            br_uf = max(br_uf, 0.f); br_vf = max(br_vf, static_cast<float>(CV_PI * projector_.scale));
        }
    }

    x = projector_.rinv[1];
    y = -projector_.rinv[4];
    z = projector_.rinv[7];
    if (y > 0.f)
    {
        float x_ = (projector_.k[0] * x + projector_.k[1] * y) / z + projector_.k[2];
        float y_ = projector_.k[4] * y / z + projector_.k[5];
        if (x_ > 0.f && x_ < src_size.width && y_ > 0.f && y_ < src_size.height)
        {
            tl_uf = min(tl_uf, 0.f); tl_vf = min(tl_vf, static_cast<float>(0));
            br_uf = max(br_uf, 0.f); br_vf = max(br_vf, static_cast<float>(0));
        }
    }

     
然后使用反影将图纸反影至变换的图像上,像从计算默认是二维线性插值。

     反影的公式:

void SphericalProjector::mapBackward(float u, float v, float &x, float &y)
{
    u /= scale;
    v /= scale;

    float sinv = sinf(static_cast<float>(CV_PI) - v);
    float x_ = sinv * sinf(u);
    float y_ = cosf(static_cast<float>(CV_PI) - v);
    float z_ = sinv * cosf(u);

    float z;
    x = k_rinv[0] * x_ + k_rinv[1] * y_ + k_rinv[2] * z_;
    y = k_rinv[3] * x_ + k_rinv[4] * y_ + k_rinv[5] * z_;
    z = k_rinv[6] * x_ + k_rinv[7] * y_ + k_rinv[8] * z_;

    if (z > 0) { x /= z; y /= z; }
    else x = y = -1;
}

然后将反投影求出的x,y值写副Mat矩阵xmap和ymap中

 xmap.create(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, CV_32F);
    ymap.create(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, CV_32F);

    float x, y;
    for (int v = dst_tl.y; v <= dst_br.y; ++v)
    {
        for (int u = dst_tl.x; u <= dst_br.x; ++u)
        {
            projector_.mapBackward(static_cast<float>(u), static_cast<float>(v), x, y);
            xmap.at<float>(v - dst_tl.y, u - dst_tl.x) = x;
            ymap.at<float>(v - dst_tl.y, u - dst_tl.x) = y;
        }
    }

末了使opencv自带的remap函数将图像重新绘制:

remap(src, dst, xmap, ymap, interp_mode, border_mode);//重映射,xmap,yamp分别是u,v反投影对应的x,y值,默认是双线性插值

       
      8.光照补偿

     
图像拼接中,由于拍摄之图纸发或坐光圈或者光的问题,导致附近图片重叠区域出现亮度差,所以在拼接时便用针对图像进行亮度上,(opencv只对重叠区域进行了亮度上,这样见面导致图像融合处虽然光照渐变,但是图像整体的光强没有和平的衔接)。

      首先,将享有图像,mask矩阵进行分块(大小在32*32附近)。

  for (int img_idx = 0; img_idx < num_images; ++img_idx)
    {
        Size bl_per_img((images[img_idx].cols + bl_width_ - 1) / bl_width_,
                        (images[img_idx].rows + bl_height_ - 1) / bl_height_);
        int bl_width = (images[img_idx].cols + bl_per_img.width - 1) / bl_per_img.width;
        int bl_height = (images[img_idx].rows + bl_per_img.height - 1) / bl_per_img.height;
        bl_per_imgs[img_idx] = bl_per_img;
        for (int by = 0; by < bl_per_img.height; ++by)
        {
            for (int bx = 0; bx < bl_per_img.width; ++bx)
            {
                Point bl_tl(bx * bl_width, by * bl_height);
                Point bl_br(min(bl_tl.x + bl_width, images[img_idx].cols),
                            min(bl_tl.y + bl_height, images[img_idx].rows));

                block_corners.push_back(corners[img_idx] + bl_tl);
                block_images.push_back(images[img_idx](Rect(bl_tl, bl_br)));
                block_masks.push_back(make_pair(masks[img_idx].first(Rect(bl_tl, bl_br)),
                                                masks[img_idx].second));
            }
        }
    }

      然后,求出自由两块图像的层区域的平分光强,

//计算每一块区域的光照均值sqrt(sqrt(R)+sqrt(G)+sqrt(B))
    //光照均值是对称矩阵,所以一次循环计算两个光照值,(i,j),与(j,i)
    for (int i = 0; i < num_images; ++i)
    {
        for (int j = i; j < num_images; ++j)
        {
            Rect roi;
            //判断image[i]与image[j]是否有重叠部分
            if (overlapRoi(corners[i], corners[j], images[i].size(), images[j].size(), roi))
            {
                subimg1 = images[i](Rect(roi.tl() - corners[i], roi.br() - corners[i]));
                subimg2 = images[j](Rect(roi.tl() - corners[j], roi.br() - corners[j]));

                submask1 = masks[i].first(Rect(roi.tl() - corners[i], roi.br() - corners[i]));
                submask2 = masks[j].first(Rect(roi.tl() - corners[j], roi.br() - corners[j]));
                intersect = (submask1 == masks[i].second) & (submask2 == masks[j].second);

                N(i, j) = N(j, i) = max(1, countNonZero(intersect));

                double Isum1 = 0, Isum2 = 0;
                for (int y = 0; y < roi.height; ++y)
                {
                    const Point3_<uchar>* r1 = subimg1.ptr<Point3_<uchar> >(y);
                    const Point3_<uchar>* r2 = subimg2.ptr<Point3_<uchar> >(y);
                    for (int x = 0; x < roi.width; ++x)
                    {
                        if (intersect(y, x))
                        {
                            Isum1 += sqrt(static_cast<double>(sqr(r1[x].x) + sqr(r1[x].y) + sqr(r1[x].z)));
                            Isum2 += sqrt(static_cast<double>(sqr(r2[x].x) + sqr(r2[x].y) + sqr(r2[x].z)));
                        }
                    }
                }
                I(i, j) = Isum1 / N(i, j);
                I(j, i) = Isum2 / N(i, j);
            }
        }
    }

     建立方程,求来每个光强的调系数

    Mat_<double> A(num_images, num_images); A.setTo(0);
    Mat_<double> b(num_images, 1); b.setTo(0);//beta*N(i,j)
    for (int i = 0; i < num_images; ++i)
    {
        for (int j = 0; j < num_images; ++j)
        {
            b(i, 0) += beta * N(i, j);
            A(i, i) += beta * N(i, j);
            if (j == i) continue;
            A(i, i) += 2 * alpha * I(i, j) * I(i, j) * N(i, j);
            A(i, j) -= 2 * alpha * I(i, j) * I(j, i) * N(i, j);
        }
    }

    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快的平分亮度值,

710官方网站 65

参数alpha和beta,默认值是alpha=0.01,beta=100.

A矩阵,大小n*n,公式图片未净

710官方网站 66

b矩阵,大小n*1,

710官方网站 67

然后根据求解矩阵

gains_矩阵,大小1*n,A*gains = B

        然后将gains_进展线性滤波

    Mat_<float> ker(1, 3);
    ker(0,0) = 0.25; ker(0,1) = 0.5; ker(0,2) = 0.25;

    int bl_idx = 0;
    for (int img_idx = 0; img_idx < num_images; ++img_idx)
    {
        Size bl_per_img = bl_per_imgs[img_idx];
        gain_maps_[img_idx].create(bl_per_img);

        for (int by = 0; by < bl_per_img.height; ++by)
            for (int bx = 0; bx < bl_per_img.width; ++bx, ++bl_idx)
                gain_maps_[img_idx](by, bx) = static_cast<float>(gains[bl_idx]);
        //用分解的核函数对图像做卷积。首先,图像的每一行与一维的核kernelX做卷积;然后,运算结果的每一列与一维的核kernelY做卷积
        sepFilter2D(gain_maps_[img_idx], gain_maps_[img_idx], CV_32F, ker, ker);
        sepFilter2D(gain_maps_[img_idx], gain_maps_[img_idx], CV_32F, ker, ker);
    }

     
然后构建一个gain_maps的老三维矩阵,gain_main[图像的个数][图像分块的行数][图像分块的列数],然后针对无一可图像的gain进行滤波。

   

   9.Seam Estimation

     缝隙估计起6种方式,默认就是第三栽方式,seam_find_type ==
“gc_color”,该办法是利用最酷流动方法检测。

 if (seam_find_type == "no")
        seam_finder = new detail::NoSeamFinder();//Stub seam estimator which does nothing.
    else if (seam_find_type == "voronoi")
        seam_finder = new detail::VoronoiSeamFinder();//Voronoi diagram-based seam estimator. 泰森多边形缝隙估计
    else if (seam_find_type == "gc_color")
    {
#ifdef HAVE_OPENCV_GPU
        if (try_gpu && gpu::getCudaEnabledDeviceCount() > 0)
            seam_finder = new detail::GraphCutSeamFinderGpu(GraphCutSeamFinderBase::COST_COLOR);
        else
#endif
            seam_finder = new detail::GraphCutSeamFinder(GraphCutSeamFinderBase::COST_COLOR);//Minimum graph cut-based seam estimator
    }
    else if (seam_find_type == "gc_colorgrad")
    {
#ifdef HAVE_OPENCV_GPU
        if (try_gpu && gpu::getCudaEnabledDeviceCount() > 0)
            seam_finder = new detail::GraphCutSeamFinderGpu(GraphCutSeamFinderBase::COST_COLOR_GRAD);
        else
#endif
            seam_finder = new detail::GraphCutSeamFinder(GraphCutSeamFinderBase::COST_COLOR_GRAD);
    }
    else if (seam_find_type == "dp_color")
        seam_finder = new detail::DpSeamFinder(DpSeamFinder::COLOR);
    else if (seam_find_type == "dp_colorgrad")
        seam_finder = new detail::DpSeamFinder(DpSeamFinder::COLOR_GRAD);
    if (seam_finder.empty())
    {
        cout << "Can't create the following seam finder '" << seam_find_type << "'\n";
        return 1;
    }

      程序解读而参见:

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也使计算上。

  // Read image and resize it if necessary
        full_img = imread(img_names[img_idx]);
        if (!is_compose_scale_set)
        {
            if (compose_megapix > 0)
                compose_scale = min(1.0, sqrt(compose_megapix * 1e6 / full_img.size().area()));
            is_compose_scale_set = true;

            // Compute relative scales
            //compose_seam_aspect = compose_scale / seam_scale;
            compose_work_aspect = compose_scale / work_scale;

            // Update warped image scale
            warped_image_scale *= static_cast<float>(compose_work_aspect);
            warper = warper_creator->create(warped_image_scale);

            // Update corners and sizes
            for (int i = 0; i < num_images; ++i)
            {
                // Update intrinsics
                cameras[i].focal *= compose_work_aspect;
                cameras[i].ppx *= compose_work_aspect;
                cameras[i].ppy *= compose_work_aspect;

                // Update corner and size
                Size sz = full_img_sizes[i];
                if (std::abs(compose_scale - 1) > 1e-1)
                {
                    sz.width = cvRound(full_img_sizes[i].width * compose_scale);//取整
                    sz.height = cvRound(full_img_sizes[i].height * compose_scale);
                }

                Mat K;
                cameras[i].K().convertTo(K, CV_32F);
                Rect roi = warper->warpRoi(sz, K, cameras[i].R);//Returns Projected image minimum bounding box
                corners[i] = roi.tl();//! the top-left corner
                sizes[i] = roi.size();//! size of the real buffer
            }
        }
        if (abs(compose_scale - 1) > 1e-1)
            resize(full_img, img, Size(), compose_scale, compose_scale);
        else
            img = full_img;
        full_img.release();
        Size img_size = img.size();

        Mat K;
        cameras[img_idx].K().convertTo(K, CV_32F);

        // Warp the current image
        warper->warp(img, K, cameras[img_idx].R, INTER_LINEAR, BORDER_REFLECT, img_warped);
        // Warp the current image mask
        mask.create(img_size, CV_8U);
        mask.setTo(Scalar::all(255));
        warper->warp(mask, K, cameras[img_idx].R, INTER_NEAREST, BORDER_CONSTANT, mask_warped);
        // Compensate exposure
        compensator->apply(img_idx, corners[img_idx], img_warped, mask_warped);//光照补偿
        img_warped.convertTo(img_warped_s, CV_16S);
        img_warped.release();
        img.release();
        mask.release();

        dilate(masks_warped[img_idx], dilated_mask, Mat());
        resize(dilated_mask, seam_mask, mask_warped.size());
        mask_warped = seam_mask & mask_warped;

     对图像进行光照补偿,就是以针对诺区域就以gain:

//块光照补偿
void BlocksGainCompensator::apply(int index, Point /*corner*/, Mat &image, const Mat &/*mask*/)
{
    CV_Assert(image.type() == CV_8UC3);

    Mat_<float> gain_map;
    if (gain_maps_[index].size() == image.size())
        gain_map = gain_maps_[index];
    else
        resize(gain_maps_[index], gain_map, image.size(), 0, 0, INTER_LINEAR);

    for (int y = 0; y < image.rows; ++y)
    {
        const float* gain_row = gain_map.ptr<float>(y);
        Point3_<uchar>* row = image.ptr<Point3_<uchar> >(y);
        for (int x = 0; x < image.cols; ++x)
        {
            row[x].x = saturate_cast<uchar>(row[x].x * gain_row[x]);
            row[x].y = saturate_cast<uchar>(row[x].y * gain_row[x]);
            row[x].z = saturate_cast<uchar>(row[x].z * gain_row[x]);
        }
    }
}

   
 进行多波段融合,首先初始化blend,确定blender的齐心协力之法门,默认是多波段融合MULTI_BAND,以及基于corners顶点和图像的轻重缓急确定最后全景图的尺寸。

   //初始化blender
        if (blender.empty())
        {
            blender = Blender::createDefault(blend_type, try_gpu);
            Size dst_sz = resultRoi(corners, sizes).size();//计算最后图像的大小
            float blend_width = sqrt(static_cast<float>(dst_sz.area())) * blend_strength / 100.f;
            if (blend_width < 1.f)
                blender = Blender::createDefault(Blender::NO, try_gpu);
            else if (blend_type == Blender::MULTI_BAND)
            {
                MultiBandBlender* mb = dynamic_cast<MultiBandBlender*>(static_cast<Blender*>(blender));
                mb->setNumBands(static_cast<int>(ceil(log(blend_width)/log(2.)) - 1.));
                LOGLN("Multi-band blender, number of bands: " << mb->numBands());
            }
            else if (blend_type == Blender::FEATHER)
            {
                FeatherBlender* fb = dynamic_cast<FeatherBlender*>(static_cast<Blender*>(blender));
                fb->setSharpness(1.f/blend_width);
                LOGLN("Feather blender, sharpness: " << fb->sharpness());
            }
            blender->prepare(corners, sizes);//根据corners顶点和图像的大小确定最终全景图的尺寸
        }

      然后针对各级幅图图形构建金字塔,层数可以由输入的系数确定,默认是5叠。

     
先对极以及图像的有余和强开处理,使其能够被2^num_bands除尽,这样才会将进行向下采样num_bands次,首先由源图像pyrDown向下采样,在由最底部的图像pyrUp向上采样,把针对应层得及采样和生采样的相减,就赢得了图像的累累分量,存储到各个一个金字塔中。然后因mask,将诸幅图像的各个层金字塔分别写副尾声的金字塔层src_pyr_laplace中。

      最后用各层的金字塔叠加,就拿走了最终整的全景图。

 blender->feed(img_warped_s, mask_warped, corners[img_idx]);//将图像写入金字塔中

      源码:

void MultiBandBlender::feed(const Mat &img, const Mat &mask, Point tl)
{
    CV_Assert(img.type() == CV_16SC3 || img.type() == CV_8UC3);
    CV_Assert(mask.type() == CV_8U);

    // Keep source image in memory with small border
    int gap = 3 * (1 << num_bands_);
    Point tl_new(max(dst_roi_.x, tl.x - gap),
                 max(dst_roi_.y, tl.y - gap));
    Point br_new(min(dst_roi_.br().x, tl.x + img.cols + gap),
                 min(dst_roi_.br().y, tl.y + img.rows + gap));

    // Ensure coordinates of top-left, bottom-right corners are divided by (1 << num_bands_).
    // After that scale between layers is exactly 2.
    //
    // We do it to avoid interpolation problems when keeping sub-images only. There is no such problem when
    // image is bordered to have size equal to the final image size, but this is too memory hungry approach.
    //将顶点调整成能被2^num_bank次方除尽的值
    tl_new.x = dst_roi_.x + (((tl_new.x - dst_roi_.x) >> num_bands_) << num_bands_);
    tl_new.y = dst_roi_.y + (((tl_new.y - dst_roi_.y) >> num_bands_) << num_bands_);
    int width = br_new.x - tl_new.x;
    int height = br_new.y - tl_new.y;
    width += ((1 << num_bands_) - width % (1 << num_bands_)) % (1 << num_bands_);
    height += ((1 << num_bands_) - height % (1 << num_bands_)) % (1 << num_bands_);
    br_new.x = tl_new.x + width;
    br_new.y = tl_new.y + height;
    int dy = max(br_new.y - dst_roi_.br().y, 0);
    int dx = max(br_new.x - dst_roi_.br().x, 0);
    tl_new.x -= dx; br_new.x -= dx;
    tl_new.y -= dy; br_new.y -= dy;

    int top = tl.y - tl_new.y;
    int left = tl.x - tl_new.x;
    int bottom = br_new.y - tl.y - img.rows;
    int right = br_new.x - tl.x - img.cols;

    // Create the source image Laplacian pyramid
    Mat img_with_border;
    copyMakeBorder(img, img_with_border, top, bottom, left, right,
                   BORDER_REFLECT);//给图像设置一个边界,BORDER_REFLECT边界颜色任意
    vector<Mat> src_pyr_laplace;
    if (can_use_gpu_ && img_with_border.depth() == CV_16S)
        createLaplacePyrGpu(img_with_border, num_bands_, src_pyr_laplace);
    else
        createLaplacePyr(img_with_border, num_bands_, src_pyr_laplace);//创建高斯金字塔,每一层保存的全是高频信息

    // Create the weight map Gaussian pyramid
    Mat weight_map;
    vector<Mat> weight_pyr_gauss(num_bands_ + 1);

    if(weight_type_ == CV_32F)
    {
        mask.convertTo(weight_map, CV_32F, 1./255.);//将mask的0,255归一化成0,1
    }
    else// weight_type_ == CV_16S
    {
        mask.convertTo(weight_map, CV_16S);
        add(weight_map, 1, weight_map, mask != 0);
    }

    copyMakeBorder(weight_map, weight_pyr_gauss[0], top, bottom, left, right, BORDER_CONSTANT);

    for (int i = 0; i < num_bands_; ++i)
        pyrDown(weight_pyr_gauss[i], weight_pyr_gauss[i + 1]);

    int y_tl = tl_new.y - dst_roi_.y;
    int y_br = br_new.y - dst_roi_.y;
    int x_tl = tl_new.x - dst_roi_.x;
    int x_br = br_new.x - dst_roi_.x;

    // Add weighted layer of the source image to the final Laplacian pyramid layer
    if(weight_type_ == CV_32F)
    {
        for (int i = 0; i <= num_bands_; ++i)
        {
            for (int y = y_tl; y < y_br; ++y)
            {
                int y_ = y - y_tl;
                const Point3_<short>* src_row = src_pyr_laplace[i].ptr<Point3_<short> >(y_);
                Point3_<short>* dst_row = dst_pyr_laplace_[i].ptr<Point3_<short> >(y);
                const float* weight_row = weight_pyr_gauss[i].ptr<float>(y_);
                float* dst_weight_row = dst_band_weights_[i].ptr<float>(y);

                for (int x = x_tl; x < x_br; ++x)
                {
                    int x_ = x - x_tl;
                    dst_row[x].x += static_cast<short>(src_row[x_].x * weight_row[x_]);
                    dst_row[x].y += static_cast<short>(src_row[x_].y * weight_row[x_]);
                    dst_row[x].z += static_cast<short>(src_row[x_].z * weight_row[x_]);
                    dst_weight_row[x] += weight_row[x_];
                }
            }
            x_tl /= 2; y_tl /= 2;
            x_br /= 2; y_br /= 2;
        }
    }
    else// weight_type_ == CV_16S
    {
        for (int i = 0; i <= num_bands_; ++i)
        {
            for (int y = y_tl; y < y_br; ++y)
            {
                int y_ = y - y_tl;
                const Point3_<short>* src_row = src_pyr_laplace[i].ptr<Point3_<short> >(y_);
                Point3_<short>* dst_row = dst_pyr_laplace_[i].ptr<Point3_<short> >(y);
                const short* weight_row = weight_pyr_gauss[i].ptr<short>(y_);
                short* dst_weight_row = dst_band_weights_[i].ptr<short>(y);

                for (int x = x_tl; x < x_br; ++x)
                {
                    int x_ = x - x_tl;
                    dst_row[x].x += short((src_row[x_].x * weight_row[x_]) >> 8);
                    dst_row[x].y += short((src_row[x_].y * weight_row[x_]) >> 8);
                    dst_row[x].z += short((src_row[x_].z * weight_row[x_]) >> 8);
                    dst_weight_row[x] += weight_row[x_];
                }
            }
            x_tl /= 2; y_tl /= 2;
            x_br /= 2; y_br /= 2;
        }
    }
}

        其中,金字塔构建的源码:

void createLaplacePyr(const Mat &img, int num_levels, vector<Mat> &pyr)
{
#ifdef HAVE_TEGRA_OPTIMIZATION
    if(tegra::createLaplacePyr(img, num_levels, pyr))
        return;
#endif

    pyr.resize(num_levels + 1);

    if(img.depth() == CV_8U)
    {
        if(num_levels == 0)
        {
            img.convertTo(pyr[0], CV_16S);
            return;
        }

        Mat downNext;
        Mat current = img;
        pyrDown(img, downNext);

        for(int i = 1; i < num_levels; ++i)
        {
            Mat lvl_up;
            Mat lvl_down;

            pyrDown(downNext, lvl_down);
            pyrUp(downNext, lvl_up, current.size());
            subtract(current, lvl_up, pyr[i-1], noArray(), CV_16S);

            current = downNext;
            downNext = lvl_down;
        }

        {
            Mat lvl_up;
            pyrUp(downNext, lvl_up, current.size());
            subtract(current, lvl_up, pyr[num_levels-1], noArray(), CV_16S);

            downNext.convertTo(pyr[num_levels], CV_16S);
        }
    }
    else
    {
        pyr[0] = img;
        //构建高斯金字塔
        for (int i = 0; i < num_levels; ++i)
            pyrDown(pyr[i], pyr[i + 1]);//先高斯滤波,在亚采样,得到比pyr【i】缩小一半的图像
        Mat tmp;
        for (int i = 0; i < num_levels; ++i)
        {
            pyrUp(pyr[i + 1], tmp, pyr[i].size());//插值(偶数行,偶数列赋值为0),然后高斯滤波,核是5*5。
            subtract(pyr[i], tmp, pyr[i]);//pyr[i] = pyr[i]-tmp,得到的全是高频信息
        }
    }
}

      最终把所有层得金字塔叠加的次序:

    Mat result, result_mask;
    blender->blend(result, result_mask);//将多层金字塔图形叠加

     源码:

void MultiBandBlender::blend(Mat &dst, Mat &dst_mask)
{
    for (int i = 0; i <= num_bands_; ++i)
        normalizeUsingWeightMap(dst_band_weights_[i], dst_pyr_laplace_[i]);

    if (can_use_gpu_)
        restoreImageFromLaplacePyrGpu(dst_pyr_laplace_);
    else
        restoreImageFromLaplacePyr(dst_pyr_laplace_);

    dst_ = dst_pyr_laplace_[0];
    dst_ = dst_(Range(0, dst_roi_final_.height), Range(0, dst_roi_final_.width));
    dst_mask_ = dst_band_weights_[0] > WEIGHT_EPS;
    dst_mask_ = dst_mask_(Range(0, dst_roi_final_.height), Range(0, dst_roi_final_.width));
    dst_pyr_laplace_.clear();
    dst_band_weights_.clear();

    Blender::blend(dst, dst_mask);
}

 

4.3测试效果….10

5.小结..11

1.概述

全景视图是赖当一个稳的观察点,能够提供水平方向上位角360度过,垂直方向达成180度的擅自浏览
(简化的全景只能提供水平方向360过的浏览)。

眼下市场被的全景摄像机主要分为两栽:鱼眼全景摄像机和多镜头全景摄像机。鱼眼全景摄像机是由于单传感器配套特殊之超广角鱼眼镜头,并借助图像校正技术还原图像的鱼眼全景摄像机。鱼眼全景摄像机最终生成的全景图像即使通过校正也还是在一定程度之失真和不自然。多镜头全景摄像机可以避鱼眼镜头图像失真的缺陷,但是还是多还是遗失吗会存在融合边缘效果不真正、角度发生过错或分融合后发”附加”感的缺撼。

正文档中根据当下所查看找到的素材,对大多镜头全景视图拼接算法原理进行简易的介绍。

2.第一步骤

2.1.图像获取

通过相机获得图像。通常要基于失真较充分之鱼眼镜头和失真较小之窄视角镜头决定算法处理方式。单镜头及多镜头相机在算法处理达成吗会见有必然差距。

2.2鲜鱼眼图像矫正

假若相机镜头吧鱼类眼镜头,则图像需要开展一定的失真展开处理。

2.3图纸匹配

根据资料图片被互重叠的有的估算图片里匹配关系。主要匹配方式分点儿种植:

A.和特征无关的配合方式。最广大的就算为相关性匹配。

B.根据特征进行匹配的不二法门。最广泛的饶为依据SIFT,SURF等资料图片被有些特征点,匹配相邻图片中之特征点,估算图像中投影变换矩阵。

2.4图拼接

依据步骤2.3所得图片相互关系,将附近图片拼接至一头。

2.5图像融合

针对拼接得到的全景图进行融合处理。

2.6全景图像投射

拿合成后底全景图投射到球面、柱面或立方体上连树立适宜的视点,实现普的视图浏览。

贪图1:opencv stitching模块进行图像拼接的处理流程

(部分步骤可选)

3.算法技术点介绍

3.1图像获取

出于鱼眼镜头与例行镜头在扭转全景图方面在好酷差别,其校正算法完全不同,因此用分开讨论。但是校正后底图像进行拼接步骤时的拍卖方法肯定水平达但通用。

A.单常规镜头拍多摆设图纸方式(手握紧)

该措施特别普遍,在当下又无线电话上均有连带全景功能。

B.差不多单正常镜头组合的相机(或单镜头旋转扫描方式)

3.2鱼眼图像矫正

只要为鱼类眼镜头采集的届的图像,必须对图像进行矫正。鱼眼镜头图像校正算法通常发生些许种:一种是球面坐标定位法,一栽是听映射法。具体推导过程呈现参考资料【1】《鱼眼照片生成全景图算法的研究及落实》,矫正效果使下图所示:

3.3图片匹配

3.3.1与特点无关的相当方式

及风味无关之相当方式大的也罢相关性匹配,一般还用于没有复杂变换的图像拼接情况下。该措施计算简单,仅为常见的灰度模板匹配。具体细节表现参考文档【2】《全景图生成技术研讨》。

710官方网站 68

3.3.2据悉特征进行匹配的艺术

基于特征的匹配首先由图像及选特征信息,然后识别出些许轴图像对应之风味信息。常用的表征信息来特点轮廓,特征曲线,特征点,多应用特征点匹配法。

展开特征点匹配的第一步是提取所有资料图片的有特征点。普遍来讲,一布置图纸所涵盖的特征点通常就是四周蕴藏较生信息量的触及,而只有经过这些有着特色的一部分,基本就可推测出整布置图。常见的特征点包括SIFT,FAST,SURF等。

710官方网站 69

由特征点由特征向量表示,所以图被每个特征点显示为一个箭头。

形成特征向量之后下一个问题不怕是什么匹配了。最核心的计可叫做“最贴近搜索”(Nearest
Neighbour),实际上为就是是寻找在128维空间上直线距离最近底底特征向量,这个求直线距离的方法跟2维无异,最近的特征向量也就吃看是互为配合。SIFT原笔者用的办法是增多了k-d
tree算法来高效率地成功高维度上之极度贴近搜索。特征点匹配效果使下图所示

710官方网站 70

3.4

图形拼接

以上述步骤中拿走了图像里的配合关系,就可以依据这些关乎展开图像的拼凑了。按照图像匹配的不比方法,拼接处理为分割点儿良接近:

A.因模板匹配的方法,可得到图片见的位移(或者包括缩放)参数,继而根据参数进行图像拼接操作;

B.因特征点匹配的法门,则动用这些匹配的点来估算只应矩阵“(Homography
Estimation),也便是拿中同样摆放通过个关联性和其它一样布置匹配的计。单应矩阵H效果如下:

由此就应矩阵H,可以用原先图像被随机像素点坐标转换为新为标点,转换后底图像即为符合拼接的结果图像。下图虽为寻来适合几何约束的特征点之后,通过就应矩阵来对联合两摆图片的内容。

710官方网站 71

710官方网站 72

3.5图像融合

图像拼接后,需要针对图像重叠部分进行融合处理。图像融合技术控制了最后图像合成质量,常用的发出平均叠加法,线性法,加权法,多段融合法等。具体表现参考文档【2】《全景图生成技术研讨》。

3.5.1等分叠加法

平均叠加法是一直对图像进行平均叠加。这是极端简便的融合方法,会面世明显的拼接缝隙。

3.5.2线性法

柱面图像的拼凑多利用精炼的丝性法。图像投射到柱面坐标下,图像中就是简简单单的纯平面平移变换,局部对准后,对重叠区域用线性法融合。该办法可柱面全景图生成,或者就具有平移变换的一定量幅图像融合。

3.5.3加以权函数法

加权函数法与线性法类似,也是广泛应用的同甘共苦方法之一。该方法会使得去除边界缝隙,但当拼合区往往出现叠影模糊的观。

3.5.4几近段落融合法(多分辨率样条)

差不多段落融合法是眼下比好之同甘共苦方法,拼接成的图像既清晰而光无缝,能幸免缝隙问题跟叠影现象。另外,如果选好的最佳缝隙线,还能处理发生一线运动物体的图像拼接。但欠法运算量大是其显著缺陷。

3.6全景图像投射

3.6.1柱面全景图

永恒视点,使相机在档次对内转悠一圆拍摄场景,得到同组有重叠区域之连年环视图像序列,将及时组图像序列无缝拼合,并投影至柱面空间坐标,就获了衣服柱面全景图。柱面投影就是提图像投影至柱面上,它是同样种透视投影而休平行投影,通俗的称就是是如果活的从影子中心立即或多或少达成观测图像于柱面上的成像。下图表示以三维空间达到之点(X,Y,Z)映射到柱面模型上得相应为柱面模型上之点(x,y,z)的进程。其中θ为观察视域中心与X轴夹角,h为柱面模型高度,(x,y,z)为(X,Y,Z)在柱面模型上的影子。

710官方网站 73

其中,

710官方网站 74

710官方网站 75

3.6.2球面全景图

球面全景图是经过求取图像投射到球面的参数,将图像投射到球面模型上,然后抱的平面反展开图就是球面全景图或者有球面全景图。

710官方网站 76

710官方网站 77

4.起源图像算法库OpenCV拼接模块

开源图像算法库OpenCV在2.4.0本子后并了一个全景图拼接模块stitch,其中一个于详细的样例代码stitching_detail.cpp简要介绍如下:

4.1
stitching_detail程序运行流程

1.命令行调用程序,输入源图像及程序的参数

2.特征点检测,判断是使用surf还是orb,默认是surf。

3.对图像的特征点进行匹配,使用以来邻和次近邻方法,将有限单极度良好的相当的采购信度保存下来。

4.针对图像进行排序和将置信度高之图像保存到同一个集中,删除置信度比较低的图像里的相当,得到能科学匹配的图像序列。这样以采购信度高于门限的兼具匹配合并至一个集合中。

5.针对性具有图像进行相机参数粗略估算,然后要来旋转矩阵

6.动光束平均法更加精准的量起旋转矩阵。

7.波形校正,水平要垂直

8.拼接

9.齐心协力,多频段融合,光照补偿,

4.2
stitching_detail程序接口介绍

img1 img2 img3输入图像

–preview因预览模式运行程序,比正规模式使及早,但输出图像分辨率低,拼接的分辨率compose_megapix设置为0.6

–try_gpu(yes|no)是否使用GPU(图形处理器),默认为no

/*挪动估计参数*/

–work_megapix <–work_megapix
>图像匹配的分辨率大小,图像的面积尺寸变为work_megapix*100000,默认为0.6

–features (surf|orb)选择surf或者orb算法进行特征点计算,默认为surf

–match_conf
特征点检测置信等级,最近邻匹配距离与坏接近邻匹配距离的比值,surf默看0.65,orb默看0.3

–conf_thresh 两轴图来自同一全景图的购买信度,默认为1.0

–ba (reproj|ray)光束平均法的误差函数选择,默认是ray方法

–ba_refine_mask (mask)—————

–wave_correct (no|horiz|vert)波形校验 水平,垂直或者没 默认是horiz

–save_graph
将匹配的图以点的样式保留至文件中,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 拼接缝像从的分寸 默认是0.1 ————

–seam (no|voronoi|gc_color|gc_colorgrad)拼接缝隙估计方法
默认是gc_color

–compose_megapix 拼接分辨率,默认为-1

–expos_comp (no|gain|gain_blocks)光照补偿方式,默认是gain_blocks

–blend (no|feather|multiband)融合方法,默认是大抵频段融合

–blend_strength 融合强度,0-100.默认是5.

–output 输出图像的文件称,默认是result,jpg

命令下实例,以及程序运行时的唤起:

710官方网站 78

相关文章