二.根本步骤..贰,输入源图像以及程序的参数

并从未都读完,可是感觉依旧有要求做八个笔记的,终归那只是小说不是小说,所以能够有稍许写多少,也总算工作总计了,最根本的是那一个幸而可以,达成具有有意义文书档案的检索,比起自个儿的word来说高级很多~~~。

1.概述..2

业已不担负图像拼接相关工作,有技巧难点请自个儿消除,谢谢。

以下壹些剧情

二.首要步骤..二

一、stitching_detail程序运营流程

以下

二.一.图像获取….二

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

全文翻译。可是本人实在是少数星星在读。。。

2.二鱼眼图像校勘….2

      贰.特征点检验,判断是应用surf照旧orb,暗中同意是surf。

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

二.3图片相配….3

     
三.对图像的特征点实行相配,使用以来邻和次近邻方法,将多个最优的相称的置信度保存下去。

 

二.四图纸拼接….三

     
四.对图像进行排序以及将置信度高的图像保存到同二个成团中,删除置信度比较低的图像间的特出,获得能科学相配的图像连串。那样将置信度高于门限的享有匹合作并到四个汇合中。

一、简介:

段1:大家的拼接不须要手动校直。

段2:依据文献,图像自动对其和拼接有两套形式:一、直接的 二、基于特征的

间接的法子:一、提供十三分规范的一直,2、不过急需发轫化处理。

根据特征的配准:一、不须要开端化,二、不过可信性上欠缺。(有不小可能是哈Rees)

段三:

大家的相称具有:

一、可信相称,在转悠,缩放和光照都不安宁的情形下。

2、给出无关图,也能找到有关图之间的涉及。【

有东西须求达成一下了,刚刚的拼接,并从未控制变量。

再试2次。

真正发现了autostitch具有越来越高的可相信性,stitch会跑出卓殊,不平稳。不知晓 stitching_detail具有多高的安居。

为此大家能够先stitch 然后再 detail 最后再 autostitch?

3、多波段融合,结果更是可相信。

四、增益补偿

五、自动校直

陆、捆绑调整, 呈现了对于几个波段的交汇图像举行多波段融合。

 

段四:

别的部分结构如下:

其次片段认证,所斟酌难点的几何学,以及选取不变特征的缘由。

其3局部介绍了:图像相称方法(RANSAC)和表达图像相称的可能率模型。

第四有的:描述图像对准算法(捆绑调整),共同优化录制头的参数。

5到7有个别:处理进度 自动校直,增益补偿,多波段融合

【八片段去哪儿了?】

九有些:结论和前程展望

 

2.5图像融合….三

     5.对具备图像实行相机参数粗略估量,然后求出旋转矩阵

贰、特征相配:

段一:

全景识别算法的率先步是在装有图像之间提取和相配SITF特征点。

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

Sift特征点使得特征尺度和样子被明确。为衡量提供了一个相似不变的结构。对于仿射变化有所一定的健壮性。【什么是放射变换?】

空中累积计算达成移动不变性。

亮度不变性由梯度(解决偏差)

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

段2:sift在转悠和标准变化时是不变的,因而得以处理缩放图像,旋转图像。

哈里斯角点图像修补的相关性特征相称技术就不能够实现。

由此笔者是否也要试一下旋转图像方向。

价值观技艺完毕的相关性在图像旋转时候是浮动的。

哈Rees不仅在打转时候变化,在尺度变换时,也是浮动的。

段三:

假诺相机绕光学中央旋转,图像的变化群是四个应和矩阵的与众不一样群。由叁个旋转矢量710官方网站 1和焦距将每种录像头参数化,就付出了成对的对应矩阵,710官方网站 2,其中

710官方网站 3(1)

并且710官方网站 4710官方网站 5是均匀的图像右侧(710官方网站 6,其中710官方网站 7是2维的图像坐标)。四参数相机模型定义为:

710官方网站 8(2)

 

对旋转使用指数表示:

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

 

段四:

 

在那些变换群中,理想的尺度下将会动用不变的图像特点。然而,在图像左边中对此小的转换表示如下:

 

710官方网站 11(4)

只怕等价于710官方网站 12其中,

710官方网站 13(5)

是经过3个关于710官方网站 14的对应线性化获得的仿射变化。那象征每种小的图像修补经过逐壹仿射变换,并且创设运用了在仿射变换下1些不变的SIFT特征。

段五:

假使从拥有n个图像中领到特征点后(线性时间内),需对特征点进行匹配。由于七个图像恐怕重叠在一个纯粹的光柱上,在特点空间内种种特征点需和它如今的k个临域点同盟(k=4),通过动用k-d树算法找到类似近期的临域点,时间复杂度为O(nlogn)。K-d树是壹种周对其的二进制空间划分,它在平均最高方差维递归划分特征空间。

 

2.陆全景图像投射….三

     6.使用光束平均法进一步精准的估摸出旋转矩阵。

3、图像相称

段一:图像相配的靶子是找到全部相称(例如重叠)图像,然后图像相称联通集会成为全景图。由于每一种图大概和任意其余三个相称,这么些难点一初叶就显示出是图像数的1遍方。【所以那一个数字应该是几?对于我们早就输入的n张图片来讲。】为了博取1个好的拼凑结果,对于图像而言,每种图像只要求和个别重叠的图像来协作。【借使重叠那拼接结果自然最佳】

段1计算:难题难度巨大,倘诺有重叠就好办。

段贰:从特征相配步骤中,已经找出大量同盟的图像间全数大批量相配点的图像。对于当前图像,咱们将m幅图形作为可能的万分图像(如果m=陆),那m幅图像与方今图像有最大数量的特点匹配点。首先,使用RANSAC算法选择一体系和图像间对应矩阵包容的内点,然后选用可能率模型做特别的印证。

段二总括:用RANSAC算法。

【备选图像有m张,那么到底要不要入选,要求RANSAC算法来决定】

三.算法技艺点介绍..三

     七.波形校对,水平依旧垂直

三.1选取RANSAC算法的硬朗对应矩阵估算

段三:RANSAC(随机取样1致性算法)【RANSAC算法是使用最少的壹组自由采集样品相称点的一种健康揣度过程,那么些看下原来的小说那那句话是想说,那种算法选用的相称点少,并且也要命结实?】用来打量图像变换参数,并找到与数据颇具最佳一致性的化解方案。【使用RANSAC算法来估算图像变换参数,来找到最棒化解方案。】在全景图的景况下,选取r=四对匹配特征点,使用直接线性别变化换(DLT)方法总结图像间的呼应矩阵H。【r是哪些?君越代表相配特征点,此处数量为四.】使用直接线性变换方式总结图像间的照应矩阵H。【什么事直接线性别变化换格局?图像间的应和矩阵又是指什么?】重复试验500次,选取内点数量最大的缓解方案(在像素相对误差范围内,其预测和H是同等的)。倘使1对男才女貌图像间的性格相配点的正确率(内点可能率)为,n次试验后找到正确变换可能率为:

710官方网站 15 (6)

 

透过大批量试行后,找到科学对应矩阵的概率非常的大。例如,对于内点可能率710官方网站 16=0.五,在500次试验后,未找到正确对应矩阵的票房价值为710官方网站 17。【这这几个地方大家是或不是足以稍微弄得小一些。以增强运营的频率,抑或然因为壹些用途,可以三番五次让那个数字变小。不过肯定已经远非要求了,又或然我们得以调动这些比重参数,在少数重合区域的图像内点可能率就会十分的大,那么实验次数就能够对应的小,来减小计算压力,对于有个别内点可能率小,地点大家得以供给尝试次数尽量的大,所以怎么是内点呢?】

 

段三计算:当内点选用某1数值的情事下,我们得以因而扩充试行次数来进步相应矩阵的相配程度。

 

段四:RANSAC算法本质上是一种估量H的采集样品方法,假诺用对数似然和的最大化代替内点数量的最大化,结果是最大似然测度(MLE)。其它,尽管变换参数的先验值是实惠的,能够计算出最大后验可能率(MAP)。这个算法被分级成为MLESAC和MAPSAC。

 

段4总括:RANSAC本质上是一种估摸图像间对应矩阵的采集样品方法。假如用对数似然和的最大化代替内点数量的最大化,结果是最大似然估计(MLE)。假如参数先验值有效,可以总结出最大后验概率MAP。对应的算法分别成为:MLESAC 和MAPSAC。但是为啥对数似然和的最大化为何能够代替内点数量的最大化。内点指得是何等?

 

三.一图像获取….三

     8.拼接

三.贰图像匹配关系验证的概率模型

 

段伍:对两两图像间是或不是留存非常关系,大家运用一密密麻麻几何一样的特点相称点(RANSAC内点)和一雨后春笋在重叠区域内,但不同的特征点(RANSAC外点)来证实。【小编就好像想起来了,RANSAC貌似是一种算法,界外的正是外点,属于不沾边的点,不相匹配的点,而RANSAC内点,便是特点相称点。】验证模型通过相比这几个科学相称发生的一层层内点和错误相称发生的1种种外点的可能率来拓展认证。

 

段5总括:是还是不是相称 通过
特征1致点和特点不一致点来判定。那下边就要说判定方法了:

 

段陆:对于壹副给定的图像,重叠区域内总的相称特征点数量为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,则足以应用贝叶斯规则(式拾、1一)总结科学图像相称的先验可能率。

710官方网站 33(10)

       =710官方网站 34(11)

万1满足710官方网站 35

 

 

710官方网站 36(12)

我们能够达成图像拼接。假定710官方网站 37710官方网站 38,进一步得出正确图像相配的判定条件:

710官方网站 39 (13)

其中,710官方网站 40=8.0,710官方网站 41=0.三。纵然在那大家挑选了710官方网站 42710官方网站 43710官方网站 44710官方网站 45710官方网站 46是大家只要的值,这一个多少大家得以转移的,只怕是可以让她更为纯粹的。比如,大家得以经过在大的数额集中总括壹部分相配点和正确的对应矩阵相平等,来猜测710官方网站 47。最终只是为着总计那样二个公式710官方网站 48,要想正确相配,必须知足内点数大于这么二个与总相称特征点的一个函数。如果低于等于,就无法达成科学相配。当中710官方网站 49。(内点数必须低于总点数。)

段七:

倘使图像间的匹配点对
显明,小编么你能够找到全景种类作为三番五次相称图像集,它能够辨别图像集中的八个全景,拒绝
不相称的噪声 图像。(见图贰)

710官方网站 50

 

710官方网站 51

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

图一,从具有图像中提取SIFT特征点。使用k-d数格外全数特征点后,对于叁个加以图像,用最多特点相配点的m幅图像举行图像相配。首先执行RANSAC算法总计出对应矩阵,然后利用概率模型验证基于内点数的图像相称,在这么些事例中,输入图像是5壹七*37四,有二4三个不错特征相称点。

710官方网站 52

710官方网站 53

图二,可识别全景图。思量叁个特色匹配点的噪声集,大家运用RANSAC算法好几率验证进度找到同样的图像相配(a),种种图像对 间的 箭头表示在图像对
间找到同样的性状相称点集,图像相称连接分量被找到(b),拼接成全景图(c);注意到该算法对 不属于全景图的噪声图像
不灵活。

 

 

捆绑调整下三个等级在处理。。。

 

3.贰鱼眼图像改进….四

     九.计出万全,多频段融合,光照补偿,

叁.三图片相配….五

二、stitching_detail程序接口介绍

三.三.1与特点非亲非故的十分情势….五

    img1 img二 img三 输入图像

三.三.二基于特征进行相配的诀要….5

     –preview
 以预览格局运作程序,比常规格局要快,但输出图像分辨率低,拼接的分辨率compose_megapix 设置为0.6

叁.4图纸拼接….陆

     –try_gpu  (yes|no)  是或不是利用GPU(图形处理器),默许为no

三.5图像融合….柒

/* 运动预计参数 */

三.伍.一平均叠加法….7

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

3.5.2线性法….7

    –features (surf|orb) 选取surf大概orb算法实行特征点总括,暗中同意为surf

三.伍.叁加权函数法….柒

    –match_conf <float>
特征点检测置信等级,近来邻相称距离与次近邻相称距离的比率,surf默许为0.65,orb私下认可为0.3

三.伍.4多段融合法(多分辨率样条)….7

    –conf_thresh <float> 两幅图来自同壹全景图的置信度,暗中同意为壹.0

三.陆全景图像投射….7

    –ba (reproj|ray) 光束平均法的截断误差函数选取,暗中同意是ray方法

三.6.一柱面全景图….七

    –ba_refine_mask (mask)  —————

叁.陆.二球面全景图….八

    –wave_correct (no|horiz|vert) 波形校验 水平,垂直可能尚未
暗许是horiz

叁.6.三多面体全景图….捌

     –save_graph <file_name>
将极度的图样以点的款型保留到文件中,
Nm代表相配的多寡,NI代表正确相称的数据,C表示置信度

四.开源图像算法库OpenCV拼接模块..玖

/*图像融合参数:*/

4.1

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

stitching_detail程序运行流程….玖

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

4.2

    采用融合的平面,暗中同意是球形

stitching_detail程序接口介绍….九

    –seam_megapix <float> 拼接缝像素的大大小小 默许是0.壹

    –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.暗中同意是伍.

   –output <result_img> 输出图像的文本名,默许是result,jpg

    命令使用实例,以及程序运营时的唤起:

710官方网站 54

3、程序代码分析

    1.参数读入

   
 程序参数读入分析,将程序运维是输入的参数以及须要拼接的图像读入内部存款和储蓄器,检查图像是或不是多于二张。

    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;
    }

      二.特征点检查测试

     
判断选择是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
    贰.叁 计算图像特征点,以及总结特征点描述子,并将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
     二.4 将源图像resize到seam_megapix*10^6,并存入image[]中

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

    三.图像相配

     
 对随意两副图形实行特征点相配,然后利用查并集法,将图纸的合营关系找出,并删除那二个不属于同一全景图的图纸。

    三.一 使用以来邻和次近邻相称,对轻易两幅图举办特征点相称。

    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

    介绍一下BestOf二NearestMatcher 函数:

      //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特征点相配),任意两幅图都进展相配(三*三=9),在那之中一-》二和二-》3头计算了3回,以一-》贰为准,:

     710官方网站 58

       三.2依据随便两幅图相配的置信度,将持有置信度高于conf_thresh(暗许是一.0)的图像合并到3个全集中。

       大家经过函数的参数
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;

       四.基于单应性矩阵粗略臆想出相机参数(焦距)

        四.一 焦距参数的测度

       
依照前边求出的随机两幅图的十三分,大家依照两幅图的单应性矩阵H,求出符合条件的f,(四副图,17个门户相当,求出七个符合条件的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;
}

      求出的焦距有7个

     710官方网站 59

      求出的焦距取中值可能平均值,然后正是全部图片的焦距。

      并创设camera参数,将矩阵写入camera:

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

     四.2依据相配的内点塑造最大生成树,然后广度优先搜索求出根节点,并求出camera的Wrangler矩阵,K矩阵以及光轴中央

      camera别的参数:

     aspect = 一.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;
}

      Wrangler矩阵的值:

     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(光束法平差)算法主即使涸泽而渔全数相机参数的同步。那是全景拼接必须的一步,因为多少个成对的单应性矩阵合成全景图时,会忽视全局的界定,造成累积固有误差。由此每3个图像都要拉长光束法平差值,使图像被伊始化成相同的旋转和焦距长度。

     
光束法平差的目的函数是一个全体鲁棒性的映射误差的平方和函数。即每3个特征点都要映射到任何的图像中,总括出使相对误差的平方和微小的照相机参数。具体的演绎进程能够参见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.

      伍.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].福特Explorer先经过奇异值分解,然后对u*vt进行罗德里gues运算,获得的rvec第叁行一个值赋给cam_params_。

     奇异值分解的定义:

在矩阵M的奇异值分解中 M = UΣV*

·U的列(columns)组成一套对M的正交”输入”或”分析”的基向量。那么些向量是MM*的特征向量。

·V的列(columns)组成一套对M的正交”输出”的基向量。这个向量是M*M的特征向量。

·Σ对角线上的因素是奇异值,可说是是在输入与出口间开展的标量的”膨胀控制”。那一个是M*M及MM*的奇异值,并与U和V的行向量绝对应。

     伍.二 删除置信度小于门限的合作对

    // 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));
        }
    }

       五.叁 使用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;
    }
}

       求出根节点,然后归1化旋转矩阵讴歌RDX

    // 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;

     陆 波形订正

   
 前边几节把相机旋转矩阵计算出来,可是还有三个要素必要思考,正是由于拍片者拍片图片的时候不肯定是程度的,轻微的倾斜会招致全景图像出现飞机曲线,因而大家要对图像举行波形改进,首借使摸索每幅图形的“上升向量”(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),图像的旋转矩阵奥德赛,图像的照相机参数矩阵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值,默认是双线性插值

       
      捌.光照补偿

     
图像拼接中,由于拍片的图片有希望因为光圈可能光线的难点,导致邻座图片重叠区域出现亮度差,所以在拼接时就须求对图像进行亮度补偿,(opencv只对重叠区域展开了亮度补偿,那样会促成图像融合处尽管光照渐变,不过图像全部的光强未有和平的接入)。

      首先,将持有图像,mask矩阵实行分块(大小在3二*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

参数阿尔法和beta,暗中同意值是阿尔法=0.0壹,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[图像的个数][图像分块的行数][图像分块的列数],然后对没1副图像的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

     

     拾.多波段融合

     
由于原先举办处理的图片都以以work_scale(3.一节有介绍)举行缩放的,所以图像的内情,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顶点和图像的大小确定最终全景图的尺寸
        }

      然后对每幅图图形塑造金字塔,层数能够由输入的周全显明,暗中同意是伍层。

     
先对终端以及图像的宽和高做处理,使其能被二^num_bands除尽,那样才能将实行向下采集样品num_bands次,首先从源图像pyrDown向下采集样品,在由最底部的图像pyrUp向上采集样品,把对应层得上采集样品和下采集样品的相减,就收获了图像的往往分量,存款和储蓄到每2个金字塔中。然后根据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);
}

 

肆.3测试效果….10

5.小结..11

1.概述

全景视图是指在三个定位的观望点,能够提供水平方向上方位角360度,垂直方向上180度的人身自由浏览
(简化的全景只好提供水平方向360度的浏览)。

此时此刻市镇中的全景摄像机首要分为二种:鱼眼全景录像机和多镜头全景录像机。鱼眼全景录像机是由单传感器配套特殊的超广角鱼近视镜头,并借助图像改良技术还原图像的鱼眼全景摄像机。鱼眼全景录像机最后生成的全景图像尽管通过校勘也照例留存一定水准的失真和不自然。多镜头全景录制机能够幸免鱼老花镜头图像失真的老毛病,但是或多或少也会存在融合边缘效果不诚实、角度有过错或瓜分融合后有”附加”感的缺撼。

正文书档案中依照当下所查找到的素材,对多镜头全景视图拼接算法原理进行简短的牵线。

2.最主要步骤

二.一.图像获取

因而相机得到图像。平时必要基于失真较大的鱼近视镜头和失真较小的窄视角镜头决定算法处理格局。单镜头和多镜头相机在算法处理上也会有早晚距离。

二.二鱼眼图像改进

若相机镜头为鱼近视镜头,则图像须求展开一定的失真展开处理。

二.叁图片相配

基于资料图片中互相重叠的有的揣摸图片间相配关系。主要相称格局分两种:

A.与风味非亲非故的十一分形式。最常见的即为相关性相配。

B.依据特征举行相称的措施。最广泛的即为依据SIFT,SUHavalF等质感图片中有些特征点,相配相邻图片中的特征点,预计图像间投影变换矩阵。

2.4图形拼接

基于步骤2.3所得图片相互关系,将紧邻图片拼接至三头。

贰.5图像融合

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

二.陆全景图像投射

将合成后的全景图投射至球面、柱面或立方体上并建立适合的视点,完成全部的视图浏览。

图1:opencv stitching模块举行图像拼接的拍卖流程

(部分步骤可选)

三.算法技术点介绍

三.一图像获取

由于鱼老花镜头和例行镜头在转变全景图方面存在一点都不小差异,其改进算法完全区别,由此需分开切磋。然而改正后的图像实行拼接步骤时的拍卖方法自然程度上可通用。

A.单常规镜头拍录多张图片情势(手持)

该格局很常见,在当下有余手机上均有相关全景成效。

B.多个常规镜头组合的相机(或单镜头旋转扫描情势)

3.2鱼眼图像修正

若为鱼近视镜头采集的到的图像,必须对图像实行校对。鱼老花镜头图像校对算法平常有二种:1种是球面坐标定位法,一种是治理映射法。具体推导进程见参考资料【一】《鱼眼照片生成全景图算法的研讨与达成》,校勘效果如下图所示:

3.3图形相配

三.3.壹与风味非亲非故的极度形式

710官方网站,与风味无关的相配情势广泛的为相关性相称,1般都用来没有复杂变换的图像拼接情状下。该措施总结不难,仅为常见的灰度模板相称。具体细节见参考文书档案【二】《全景图生成技术钻探》。

710官方网站 68

3.3.2据他们说特征进行相称的主意

据说特征的非凡首先从图像上摘取特征音讯,然后识别出两幅图像对应的风味音讯。常用的表征消息有特色轮廓,特征曲线,特征点,多选取特征点相配法。

展开特征点相称的率先步是领取全体素材图片的部分特征点。普遍来讲,一张图纸所涵盖的特征点日常就是周边蕴藏较大音信量的点,而仅透过这么些具有特色的局地,基本就能够推论出整张图片。常见的特征点包罗SIFT,FAST,SURF等。

710官方网站 69

由于特征点由特征向量表示,所以图中各个特征点展现为三个箭头。

形成特征向量之后下3个难点正是怎么样匹配了。最大旨的办法能够称作“最贴近搜索”(Nearest
Neighbour),实际上也正是找在128维空间上直线距离近来的的特征向量,那个求直线距离的诀窍和二维无差异,近来的特征向量也就被认为是互相合作。SIFT原来的小说者利用的章程是扩大了k-d
tree算法来高效能地做到高维度上的最贴近搜索。特征点相称效果如下图所示

710官方网站 70

3.4

图表拼接

在以上步骤中得到了图像间的相配关系,就足以依据这几个涉及进展图像的拼凑了。根据图像相称的比不上格局,拼接处理也分两大类:

A.依据模板相配的格局,可得到图片见的活动(或然包涵缩放)参数,继而依据参数进行图像拼接操作;

B.依照特征点相配的法门,则采纳这一个相配的点来估算单应矩阵“(Homography
Estimation),也便是把内部一张通过个关联性和另一张相配的方法。单应矩阵H效果如下:

透过单应矩阵H,能够将原图像中专断像素点坐标转换为新坐标点,转换后的图像即为适合拼接的结果图像。下图即为找出适合几何约束的特征点之后,通过单应矩阵来对齐两张图片的剧情。

710官方网站 71

710官方网站 72

三.5图像融合

图像拼接后,要求对图像重叠部分开始展览融合处理。图像融合技术控制了最终图像合成品质,常用的有平均叠加法,线性法,加权法,多段融合法等。具体见参考文书档案【2】《全景图生成技术商量》。

叁.五.1平均叠加法

平均叠加法是直接对图像举行平均叠加。那是最简单易行的众志成城方法,会产出显著的拼凑缝隙。

3.5.2线性法

柱面图像的拼凑多使用简单的线性法。图像投射到柱面坐标下,图像间正是简单的纯平面平移变换,局地对准后,对重叠区域用线性法融合。该方法适合柱面全景图生成,或许仅具有平移变换的两幅图像融合。

3.五.3加权函数法

加权函数法与线性法类似,也是广泛应用的丹舟共济方法之1。该方法能使得去除边界缝隙,但在拼合区往往出现叠影模糊的景色。

三.5.四多段融合法(多分辨率样条)

多段融合法是日前可比好的融合方法,拼接成的图像既清晰又光滑无缝,能防止缝隙难题和叠影现象。此外,要是选拔好的最棒缝隙线,还是能处理有一线运动物体的图像拼接。但该措施运算量大是其通晓缺欠。

三.陆全景图像投射

3.六.1柱面全景图

永恒视点,使相机在档次面内转悠7日拍戏场景,获得壹组具有重叠区域的总是环视图像系列,将那组图像种类无缝拼合,并投影到柱面空间坐标,就获得了衣服柱面全景图。柱面投影便是讲图像投影到柱面上,它是壹种透视投影而非平行投影,通俗的讲正是要活的从事电影工作子核心那或多或少上观看比赛图像在柱面上的成像。下图表示将三维空间上的点(X,Y,Z)映射到柱面模型上赢得相应于柱面模型上的点(x,y,z)的经过。当中θ为观看视域大旨与X轴夹角,h为柱面模型中度,(x,y,z)为(X,Y,Z)在柱面模型上的阴影。

710官方网站 73

其中,

710官方网站 74

710官方网站 75

叁.6.贰球面全景图

球面全景图是透过求取图像投射到球面包车型地铁参数,将图像投射到球面模型上,然后拿走的平面反展开图便是球面全景图只怕有个别球面全景图。

710官方网站 76

710官方网站 77

四.开源图像算法库OpenCV拼接模块

开源图像算法库OpenCV在贰.四.0版本后集成了三个全景图拼接模块stitch,个中一个较详细的样例代码stitching_detail.cpp简要介绍如下:

4.1
stitching_detail程序运营流程

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

2.特征点检查评定,判断是使用surf依旧orb,默许是surf。

三.对图像的特征点举行相称,使用以来邻和次近邻方法,将八个最优的同盟的置信度保存下去。

四.对图像进行排序以及将置信度高的图像保存到同三个凑合中,删除置信度相比低的图像间的很是,获得能正确相称的图像系列。那样将置信度高于门限的有着匹协作并到贰个汇聚中。

伍.对富有图像实行相机参数粗略估摸,然后求出旋转矩阵

陆.利用光束平均法进一步精准的预计出旋转矩阵。

7.波形改正,水平依然垂直

8.拼接

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

4.2
stitching_detail程序接口介绍

img一 img二 img三输入图像

–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.陆5,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.壹 ————

–seam (no|voronoi|gc_color|gc_colorgrad)拼接缝隙估量方法
暗中认可是gc_color

–compose_megapix 拼接分辨率,暗中同意为-一

–expos_comp (no|gain|gain_blocks)光照补偿格局,默许是gain_blocks

–blend (no|feather|multiband)融合方法,暗中认可是多频段融合

–blend_strength 融合强度,0-100.私下认可是5.

–output 输出图像的文书名,暗许是result,jpg

一声令下使用实例,以及程序运转时的升迁:

710官方网站 78

相关文章