精华内容
下载资源
问答
  • OpenCV 常用算法

    千次阅读 2008-09-24 15:02:00
    OpenCV 常用算法function StorePage(){d=document;t=d.selection?(d.selection.type!=None?d.selection.createRange().text:):(d.getSelection?d.getSelection():);void(keyit=window.open(http://www.365key

    OpenCV 常用算法

    <script>function StorePage(){d=document;t=d.selection?(d.selection.type!='None'?d.selection.createRange().text:''):(d.getSelection?d.getSelection():'');void(keyit=window.open('http://www.365key.com/storeit.aspx?t='+escape(d.title)+'&u='+escape(d.location.href)+'&c='+escape(t),'keyit','scrollbars=no,width=475,height=575,left=75,top=20,status=no,resizable=yes'));keyit.focus();}</script>  
    6月8日

    关于直方图

        图像直方图是图像处理中一种十分重要的图像分析工具,它描述了一幅图像的灰度级内 容。
       
        从数学上来说图像直方图是图像各灰度值统计特性与图像灰度值的函数,它统计一幅图像中各个灰度级出现的次数或概率;
       
        从图形上来说,它是一个二维图, 横坐标表示图像中各个像素点的灰度级,纵坐标为各个灰度级上图像各个像素点出现的 次数或概率。任何一幅图像的直方图都包含了丰富的信息,它主要用在图象分割,图像灰度变换等处理过程中。
       
        最合适直方图的典型特征是从左端到右端的每一个辉度级别都有相对应的图形存在,也就是说,值是平均分布的,而且图像应该服从正态分布,即中间大,两头小的情况。
     

    1.正常曝光的照片

    2.
    曝光不足的照片


    3.
    曝光过渡的照片

    4.
    反差过低的照片

    5.
    反差过高的照片

     
     

        
        没有任何两幅图像的直方图是完全一样的。 因此,利用直方图可以在目标跟踪时有所作为。
    6月1日

    你够鲁棒吗

    鲁棒是Robust的音译,也就是健壮和强壮的意思。
    鲁棒性(robustness)就是系统的健壮性。它是在异常和危险情况下系统生存的关键。比如说,计算机软件在输入错误、磁盘故障、网络过载或有意攻击情况下,能否不死机、不崩溃,就是该软件的鲁棒性。所谓“鲁棒性”,是指控制系统在一定(结构,大小)的参数摄动下,维持某些性能的特性。根据对性能的不同定义,可分为稳定鲁棒性和性能鲁棒性。以闭环系统的鲁棒性作为目标设计得到的固定控制器称为鲁棒控制器。

    4月12日

    OpenCV学习笔记(五)卡尔曼滤波器

    卡尔曼滤波器 Kalman Filter

    1    什么是卡尔曼滤波器

    What is the Kalman Filter?

    在学习卡尔曼滤波器之前,首先看看为什么叫“卡尔曼”。跟其他著名的理论(例如傅立叶变换,泰勒级数等等)一样,卡尔曼也是一个人的名字,而跟他们不同的是,他是个现代人!

    卡尔曼全名Rudolf Emil Kalman,匈牙利数学家,1930年出生于匈牙利首都布达佩斯。19531954年于麻省理工学院分别获得电机工程学士及硕士学位。1957年于哥伦比亚大学获得博士学位。我们现在要学习的卡尔曼滤波器,正是源于他的博士论文和1960年发表的论文《A New Approach to Linear Filtering and Prediction Problems》(线性滤波与预测问题的新方法)。如果对这编论文有兴趣,可以到这里的地址下载: http://www.cs.unc.edu/~welch/media/pdf/Kalman1960.pdf

    简单来说,卡尔曼滤波器是一个“optimal recursive data processing algorithm(最优化自回归数据处理算法)”。对于解决很大部分的问题,他是最优,效率最高甚至是最有用的。他的广泛应用已经超过30年,包括机器人导航,控制,传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。近年来更被应用于计算机图像处理,例如头脸识别,图像分割,图像边缘检测等等。

    2.卡尔曼滤波器的介绍

    Introduction to the Kalman Filter

    为了可以更加容易的理解卡尔曼滤波器,这里会应用形象的描述方法来讲解,而不是像大多数参考书那样罗列一大堆的数学公式和数学符号。但是,他的5条公式是其核心内容。结合现代的计算机,其实卡尔曼的程序相当的简单,只要你理解了他的那5条公式。

    在介绍他的5条公式之前,先让我们来根据下面的例子一步一步的探索。

    假设我们要研究的对象是一个房间的温度。根据你的经验判断,这个房间的温度是恒定的,也就是下一分钟的温度等于现在这一分钟的温度(假设我们用一分钟来做时间单位)。假设你对你的经验不是100%的相信,可能会有上下偏差几度。我们把这些偏差看成是高斯白噪声(White Gaussian Noise),也就是这些偏差跟前后时间是没有关系的而且符合高斯分配(Gaussian Distribution)。另外,我们在房间里放一个温度计,但是这个温度计也不准确的,测量值会比实际值偏差。我们也把这些偏差看成是高斯白噪声。

    好了,现在对于某一分钟我们有两个有关于该房间的温度值:你根据经验的预测值(系统的预测值)和温度计的值(测量值)。下面我们要用这两个值结合他们各自的噪声来估算出房间的实际温度值。

    假如我们要估算k时刻的是实际温度值。首先你要根据k-1时刻的温度值,来预测k时刻的温度。因为你相信温度是恒定的,所以你会得到k时刻的温度预测值是跟k-1时刻一样的,假设是23度,同时该值的高斯噪声的偏差是5度(5是这样得到的:如果k-1时刻估算出的最优温度值的偏差是3,你对自己预测的不确定度是4度,他们平方相加再开方,就是5)。然后,你从温度计那里得到了k时刻的温度值,假设是25度,同时该值的偏差是4度。

    由于我们用于估算k时刻的实际温度有两个温度值,分别是23度和25度。究竟实际温度是多少呢?相信自己还是相信温度计呢?究竟相信谁多一点,我们可以用他们的covariance来判断。因为Kg^2=5^2/(5^2+4^2),所以Kg=0.78,我们可以估算出k时刻的实际温度值是:23+0.78*(25-23)=24.56度。可以看出,因为温度计的covariance比较小(比较相信温度计),所以估算出的最优温度值偏向温度计的值。

    现在我们已经得到k时刻的最优温度值了,下一步就是要进入k+1时刻,进行新的最优估算。到现在为止,好像还没看到什么自回归的东西出现。对了,在进入k+1时刻之前,我们还要算出k时刻那个最优值(24.56度)的偏差。算法如下:((1-Kg)*5^2)^0.5=2.35。这里的5就是上面的k时刻你预测的那个23度温度值的偏差,得出的2.35就是进入k+1时刻以后k时刻估算出的最优温度值的偏差(对应于上面的3)。

    就是这样,卡尔曼滤波器就不断的把covariance递归,从而估算出最优的温度值。他运行的很快,而且它只保留了上一时刻的covariance。上面的Kg,就是卡尔曼增益(Kalman Gain)。他可以随不同的时刻而改变他自己的值,是不是很神奇!

    下面就要言归正传,讨论真正工程系统上的卡尔曼。

    3    卡尔曼滤波器算法

    The Kalman Filter Algorithm

    在这一部分,我们就来描述源于Dr Kalman 的卡尔曼滤波器。下面的描述,会涉及一些基本的概念知识,包括概率(Probability),随即变量(Random Variable),高斯或正态分配(Gaussian Distribution)还有State-space Model等等。但对于卡尔曼滤波器的详细证明,这里不能一一描述。

    首先,我们先要引入一个离散控制过程的系统。该系统可用一个线性随机微分方程(Linear Stochastic Difference equation)来描述:

    X(k)=A X(k-1)+B U(k)+W(k)

    再加上系统的测量值:

    Z(k)=H X(k)+V(k)

    上两式子中,X(k)k时刻的系统状态,U(k)k时刻对系统的控制量。AB是系统参数,对于多模型系统,他们为矩阵。Z(k)k时刻的测量值,H是测量系统的参数,对于多测量系统,H为矩阵。W(k)V(k)分别表示过程和测量的噪声。他们被假设成高斯白噪声(White Gaussian Noise),他们的covariance 分别是QR(这里我们假设他们不随系统状态变化而变化)。

    对于满足上面的条件(线性随机微分系统,过程和测量都是高斯白噪声),卡尔曼滤波器是最优的信息处理器。下面我们来用他们结合他们的covariances 来估算系统的最优化输出(类似上一节那个温度的例子)。

    首先我们要利用系统的过程模型,来预测下一状态的系统。假设现在的系统状态是k,根据系统的模型,可以基于系统的上一状态而预测出现在状态:

    X(k|k-1)=A X(k-1|k-1)+B U(k) ……….. (1)

    (1)中,X(k|k-1)是利用上一状态预测的结果,X(k-1|k-1)是上一状态最优的结果,U(k)为现在状态的控制量,如果没有控制量,它可以为0

    到现在为止,我们的系统结果已经更新了,可是,对应于X(k|k-1)covariance还没更新。我们用P表示covariance

    P(k|k-1)=A P(k-1|k-1) A+Q ……… (2)

    (2)中,P(k|k-1)X(k|k-1)对应的covarianceP(k-1|k-1)X(k-1|k-1)对应的covarianceA’表示A的转置矩阵,Q是系统过程的covariance。式子12就是卡尔曼滤波器5个公式当中的前两个,也就是对系统的预测。

    现在我们有了现在状态的预测结果,然后我们再收集现在状态的测量值。结合预测值和测量值,我们可以得到现在状态(k)的最优化估算值X(k|k)

    X(k|k)= X(k|k-1)+Kg(k) (Z(k)-H X(k|k-1)) ……… (3)

    其中Kg为卡尔曼增益(Kalman Gain)

    Kg(k)= P(k|k-1) H / (H P(k|k-1) H + R) ……… (4)

    到现在为止,我们已经得到了k状态下最优的估算值X(k|k)。但是为了要另卡尔曼滤波器不断的运行下去直到系统过程结束,我们还要更新k状态下X(k|k)covariance

    P(k|k)=I-Kg(k) HP(k|k-1) ……… (5)

    其中I 1的矩阵,对于单模型单测量,I=1。当系统进入k+1状态时,P(k|k)就是式子(2)P(k-1|k-1)。这样,算法就可以自回归的运算下去。

    卡尔曼滤波器的原理基本描述了,式子12345就是他的5 个基本公式。根据这5个公式,可以很容易的实现计算机的程序

     

    4月4日

    OpenCV学习笔记(四)运动物体跟踪的camshift算法

    CamShift算法
    简介
    CamShift算法,即"Continuously Apative Mean-Shift"算法,是一种运动跟踪算法。它主要通过视频图像中运动物体的颜色信息来达到跟踪的目的。我把这个算法分解成三个部分,便于理解:

    Back Projection计算。

    Mean Shift算法

    CamShift算法

    1 Back Projection计算
    计算Back Projection的步骤是这样的:

    1. 计算被跟踪目标的色彩直方图。在各种色彩空间中,只有HSI空间(或与HSI类似的色彩空间)中的H分量可以表示颜色信息。所以在具体的计算过程中,首先将其他的色彩空间的值转化到HSI空间,然后会其中的H分量做1D直方图计算。

    2. 根据获得的色彩直方图将原始图像转化成色彩概率分布图像,这个过程就被称作"Back Projection"。

    在OpenCV中的直方图函数中,包含Back Projection的函数,函数原型是:

       void cvCalcBackProject(IplImage** img, CvArr** backproject, const CvHistogram* hist);

    传递给这个函数的参数有三个:

    1. IplImage** img:存放原始图像,输入。

    2. CvArr** backproject:存放Back Projection结果,输出。

    3. CvHistogram* hist:存放直方图,输入

    下面就给出计算Back Projection的OpenCV代码。

    1.准备一张只包含被跟踪目标的图片,将色彩空间转化到HSI空间,获得其中的H分量:

      IplImage* target=cvLoadImage("target.bmp",-1);  //装载图片

      IplImage* target_hsv=cvCreateImage( cvGetSize(target), IPL_DEPTH_8U, 3 );

      IplImage* target_hue=cvCreateImage( cvGetSize(target), IPL_DEPTH_8U, 3 );

      cvCvtColor(target,target_hsv,CV_BGR2HSV);       //转化到HSV空间

      cvSplit( target_hsv, target_hue, NULL, NULL, NULL );    //获得H分量

    2.计算H分量的直方图,即1D直方图:

      IplImage* h_plane=cvCreateImage( cvGetSize(target_hsv),IPL_DEPTH_8U,1 );

      int hist_size[]={255};          //将H分量的值量化到[0,255]

      float* ranges[]={ {0,360} };    //H分量的取值范围是[0,360)

      CvHistogram* hist=cvCreateHist(1, hist_size, ranges, 1);

      cvCalcHist(&target_hue, hist, 0, NULL);

    在这里需要考虑H分量的取值范围的问题,H分量的取值范围是[0,360),这个取值范围的值不能用一个byte来表示,为了能用一个byte表示,需要将H值做适当的量化处理,在这里我们将H分量的范围量化到[0,255].

    4.计算Back Projection:

      IplImage* rawImage;

      //----------------------------------------------

      //get from video frame,unsigned byte,one channel

      //----------------------------------------------

      IplImage* result=cvCreateImage(cvGetSize(rawImage),IPL_DEPTH_8U,1);

      cvCalcBackProject(&rawImage,result,hist);

    5.结果:result即为我们需要的.

    2) Mean Shift算法
     

    这里来到了CamShift算法,OpenCV实现的第二部分,这一次重点讨论Mean Shift算法。

    在讨论Mean Shift算法之前,首先讨论在2D概率分布图像中,如何计算某个区域的重心(Mass Center)的问题,重心可以通过以下公式来计算:

    1.计算区域内0阶矩

    for(int i=0;i<height;i++)

      for(int j=0;j<width;j++)

         M00+=I(i,j)

    2.区域内1阶矩:

    for(int i=0;i<height;i++)

      for(int j=0;j<width;j++)

      {

        M10+=i*I(i,j);

        M01+=j*I(i,j);

      }

    3.则Mass Center为:

    Xc=M10/M00; Yc=M01/M00

    接下来,讨论Mean Shift算法的具体步骤,Mean Shift算法可以分为以下4步:

    1.选择窗的大小和初始位置.

    2.计算此时窗口内的Mass Center.

    3.调整窗口的中心到Mass Center.

    4.重复2和3,直到窗口中心"会聚",即每次窗口移动的距离小于一定的阈值。

    在OpenCV中,提供Mean Shift算法的函数,函数的原型是:

    int cvMeanShift(IplImage* imgprob,CvRect windowIn,

                        CvTermCriteria criteria,CvConnectedComp* out);

    需要的参数为:

    1.IplImage* imgprob:2D概率分布图像,传入;

    2.CvRect windowIn:初始的窗口,传入;

    3.CvTermCriteria criteria:停止迭代的标准,传入;

    4.CvConnectedComp* out:查询结果,传出。

    (注:构造CvTermCriteria变量需要三个参数,一个是类型,另一个是迭代的最大次数,最后一个表示特定的阈值。例如可以这样构造criteria:criteria=cvTermCriteria(CV_TERMCRIT_ITER|CV_TERMCRIT_EPS,10,0.1)。)

    返回的参数:

    1.int:迭代的次数。

    实现代码:暂时缺

    3) CamShift算法
    1.原理

    在了解了MeanShift算法以后,我们将MeanShift算法扩展到连续图像序列(一般都是指视频图像序列),这样就形成了CamShift算法。CamShift算法的全称是"Continuously Apaptive Mean-SHIFT",它的基本思想是视频图像的所有帧作MeanShift运算,并将上一帧的结果(即Search Window的中心和大小)作为下一帧MeanShift算法的Search Window的初始值,如此迭代下去,就可以实现对目标的跟踪。整个算法的具体步骤分5步:

    Step 1:将整个图像设为搜寻区域。

    Step 2:初始话Search Window的大小和位置。

    Step 3:计算Search Window内的彩色概率分布,此区域的大小比Search Window要稍微大一点。

    Step 4:运行MeanShift。获得Search Window新的位置和大小。

    Step 5:在下一帧视频图像中,用Step 3获得的值初始化Search Window的位置和大小。跳转到Step 3继续运行。

    2.实现

    在OpenCV中,有实现CamShift算法的函数,此函数的原型是:

      cvCamShift(IplImage* imgprob, CvRect windowIn,

                    CvTermCriteria criteria,

                    CvConnectedComp* out, CvBox2D* box=0);

    其中:

       imgprob:色彩概率分布图像。

       windowIn:Search Window的初始值。

       Criteria:用来判断搜寻是否停止的一个标准。

       out:保存运算结果,包括新的Search Window的位置和面积。

       box:包含被跟踪物体的最小矩形。

    说明:

    1.在OpenCV 4.0 beta的目录中,有CamShift的例子。遗憾的是这个例子目标的跟踪是半自动的,即需要人手工选定一个目标。我正在努力尝试全自动的目标跟踪,希望可以和大家能在这方面与大家交流。


     

    3月31日

    OpenCV学习笔记(三)人脸检测的代码分析

     
    OpenCV学习笔记(三)人脸检测的代码分析
    一、预备知识:
    1、动态内存存储及操作函数
    CvMemStorage
    typedef struct CvMemStorage
    {
        struct CvMemBlock* bottom;/* first allocated block */
        struct CvMemBlock* top; /* the current memory block - top of the stack */
        struct CvMemStorage* parent; /* borrows new blocks from */
        int block_size; /* block size */
        int free_space; /* free space in the top block (in bytes) */
    } CvMemStorage;
    内存存储器是一个可用来存储诸如序列,轮廓,图形,子划分等动态增长数据结构的底层结构。它是由一系列以同等大小的内存块构成,呈列表型 ---bottom 域指的是列首,top 域指的是当前指向的块但未必是列尾.在bottom和top之间所有的块(包括bottom, 不包括top)被完全占据了空间;在 top和列尾之间所有的块(包括块尾,不包括top)则是空的;而top块本身则被占据了部分空间 -- free_space 指的是top块剩余的空字节数。新分配的内存缓冲区(或显示的通过 cvMemStorageAlloc 函数分配,或隐示的通过 cvSeqPush, cvGraphAddEdge等高级函数分配)总是起始于当前块(即top块)的剩余那部分,如果剩余那部分能满足要求(够分配的大小)。分配后,free_space 就减少了新分配的那部分内存大小,外加一些用来保存适当列型的附加大小。当top块的剩余空间无法满足被分配的块(缓冲区)大小时,top块的下一个存储块被置为当前块(新的top块) --  free_space 被置为先前分配的整个块的大小。如果已经不存在空的存储块(即:top块已是列尾),则必须再分配一个新的块(或从parent那继承,见 cvCreateChildMemStorage)并将该块加到列尾上去。于是,存储器(memory storage)就如同栈(Stack)那样, bottom指向栈底,(top, free_space)对指向栈顶。栈顶可通过 cvSaveMemStoragePos保存,通过 cvRestoreMemStoragePos 恢复指向, 通过 cvClearStorage 重置。
    CvMemBlock
    内存存储块结构
    typedef struct CvMemBlock
    {
        struct CvMemBlock* prev;
        struct CvMemBlock* next;
    } CvMemBlock;
    CvMemBlock 代表一个单独的内存存储块结构。 内存存储块中的实际数据存储在 header块 之后(即:存在一个头指针 head 指向的块 header ,该块不存储数据),于是,内存块的第 i 个字节可以通过表达式 ((char*)(mem_block_ptr+1))[i] 获得。然而,通常没必要直接去获得存储结构的域。
    CvMemStoragePos
    内存存储块地址
    typedef struct CvMemStoragePos
    {
        CvMemBlock* top;
        int free_space;
    } CvMemStoragePos;
    该结构(如以下所说)保存栈顶的地址,栈顶可以通过 cvSaveMemStoragePos 保存,也可以通过 cvRestoreMemStoragePos 恢复。
    ________________________________________
    cvCreateMemStorage
    创建内存块
    CvMemStorage* cvCreateMemStorage( int block_size=0 );
     block_size:存储块的大小以字节表示。如果大小是 0 byte, 则将该块设置成默认值  当前默认大小为64k.
    函数 cvCreateMemStorage 创建一内存块并返回指向块首的指针。起初,存储块是空的。头部(即:header)的所有域值都为 0,除了 block_size 外.
    ________________________________________
    cvCreateChildMemStorage
    创建子内存块
    CvMemStorage* cvCreateChildMemStorage( CvMemStorage* parent );
    parent    父内存块
    函数 cvCreateChildMemStorage 创建一类似于普通内存块的子内存块,除了内存分配/释放机制不同外。当一个子存储块需要一个新的块加入时,它就试图从parent 那得到这样一个块。如果 parent 中 还未被占据空间的那些块中的第一个块是可获得的,就获取第一个块(依此类推),再将该块从 parent  那里去除。如果不存在这样的块,则 parent 要么分配一个,要么从它自己 parent (即:parent 的 parent) 那借个过来。换句话说,完全有可能形成一个链或更为复杂的结构,其中的内存存储块互为 child/ parent 关系(父子关系)。当子存储结构被释放或清除,它就把所有的块还给各自的 parent. 在其他方面,子存储结构同普通存储结构一样。
    子存储结构在下列情况中是非常有用的。想象一下,如果用户需要处理存储在某个块中的动态数据,再将处理的结果存放在该块中。在使用了最简单的方法处理后,临时数据作为输入和输出数据被存放在了同一个存储块中,于是该存储块看上去就类似下面处理后的样子: Dynamic data processing without using child storage. 结果,在存储块中,出现了垃圾(临时数据)。然而,如果在开始处理数据前就先建立一个子存储块,将临时数据写入子存储块中并在最后释放子存储块,那么最终在 源/目的存储块 (source / destination storage) 中就不会出现垃圾, 于是该存储块看上去应该是如下形式:Dynamic data processing using a child storage.
    cvReleaseMemStorage
    释放内存块
    void cvReleaseMemStorage( CvMemStorage** storage );
    storage: 指向被释放了的存储块的指针
    函数 cvReleaseMemStorage 释放所有的存储(内存)块 或者 将它们返回给各自的 parent(如果需要的话)。 接下来再释放 header块(即:释放头指针 head 指向的块 = free(head))并清除指向该块的指针(即:head = NULL)。在释放作为 parent 的块之前,先清除各自的 child 块。
    cvClearMemStorage
    清空内存存储块
    void cvClearMemStorage( CvMemStorage* storage );
    storage:存储存储块
    函数 cvClearMemStorage 将存储块的 top 置到存储块的头部(注:清空存储块中的存储内容)。该函数并不释放内存(仅清空内存)。假使该内存块有一个父内存块(即:存在一内存块与其有父子关系),则函数就将所有的块返回给其 parent.
    cvMemStorageAlloc
    在存储块中分配以内存缓冲区
    void* cvMemStorageAlloc( CvMemStorage* storage, size_t size );
    storage:内存块.
    size:缓冲区的大小.
    函数 cvMemStorageAlloc 在存储块中分配一内存缓冲区。该缓冲区的大小不能超过内存块的大小,否则就会导致运行时错误。缓冲区的地址被调整为CV_STRUCT_ALIGN 字节 (当前为 sizeof(double)).
    cvMemStorageAllocString
    在存储块中分配一文本字符串
    typedef struct CvString
    {
        int len;
        char* ptr;
    }
    CvString;
    CvString cvMemStorageAllocString( CvMemStorage* storage, const char* ptr, int len=-1 );
    storage:存储块
    ptr:字符串
    len:字符串的长度(不计算'/0')。如果参数为负数,函数就计算该字符串的长度。
    函数 cvMemStorageAlloString 在存储块中创建了一字符串的拷贝。它返回一结构,该结构包含字符串的长度(该长度或通过用户传递,或通过计算得到)和指向被拷贝了的字符串的指针。
    cvSaveMemStoragePos
    保存内存块的位置(地址)
    void cvSaveMemStoragePos( const CvMemStorage* storage, CvMemStoragePos* pos );
    storage:内存块.
    pos:内存块顶部位置。
    函数 cvSaveMemStoragePos 将存储块的当前位置保存到参数 pos 中。 函数 cvRestoreMemStoragePos 可进一步获取该位置(地址)。
    cvRestoreMemStoragePos
    恢复内存存储块的位置
    void cvRestoreMemStoragePos( CvMemStorage* storage, CvMemStoragePos* pos );
    storage:内存块.
    pos:新的存储块的位置
    函数 cvRestoreMemStoragePos 通过参数 pos 恢复内存块的位置。该函数和函数 cvClearMemStorage 是释放被占用内存块的唯一方法。注意:没有什么方法可去释放存储块中被占用的部分内存。
    2、分类器结构及操作函数:
    CvHaarFeature
    #define CV_HAAR_FEATURE_MAX  3
    typedef struct CvHaarFeature
    {
        int  tilted; 
        struct
        {
            CvRect r;
            float weight;
    } rect[CV_HAAR_FEATURE_MAX];
     
    }
    CvHaarFeature;
    一个 harr 特征由 2-3 个具有相应权重的矩形组成
    titled :/* 0 means up-right feature, 1 means 45--rotated feature */
    rect[CV_HAAR_FEATURE_MAX];  /* 2-3 rectangles with weights of opposite signs and       with absolute values inversely proportional to the areas of the rectangles. if rect[2].weight !=0, then  the feature consists of 3 rectangles, otherwise it consists of 2 */
    CvHaarClassifier
    typedef struct CvHaarClassifier
    {
        int count; 
        CvHaarFeature* haar_feature;
        float* threshold;
        int* left;
        int* right;
        float* alpha;
    }
    CvHaarClassifier;
    /* a single tree classifier (stump in the simplest case) that returns the response for the feature   at the particular image location (i.e. pixel sum over subrectangles of the window) and gives out a value depending on the responce */
    int count;  /* number of nodes in the decision tree */
     /* these are "parallel" arrays. Every index i corresponds to a node of the decision tree (root has 0-th index).
    left[i] - index of the left child (or negated index if the left child is a leaf)
    right[i] - index of the right child (or negated index if the right child is a leaf)
    threshold[i] - branch threshold. if feature responce is <= threshold, left branch                      is chosen, otherwise right branch is chosed.
    alpha[i] - output value correponding to the leaf. */
    CvHaarStageClassifier
    typedef struct CvHaarStageClassifier
    {
        int  count;  /* number of classifiers in the battery */
        float threshold; /* threshold for the boosted classifier */
        CvHaarClassifier* classifier; /* array of classifiers */
        /* these fields are used for organizing trees of stage classifiers,
           rather than just stright cascades */
        int next;
        int child;
        int parent;
    }
    CvHaarStageClassifier;

    /* a boosted battery of classifiers(=stage classifier): the stage classifier returns 1 if the sum of the classifiers' responces is greater than threshold and 0 otherwise */
    int  count;  /* number of classifiers in the battery */
    float threshold; /* threshold for the boosted classifier */
    CvHaarClassifier* classifier; /* array of classifiers */
    /* these fields are used for organizing trees of stage classifiers, rather than just stright cascades */
    CvHaarClassifierCascade
    typedef struct CvHidHaarClassifierCascade CvHidHaarClassifierCascade;
    typedef struct CvHaarClassifierCascade
    {
        int  flags;
        int  count;
        CvSize orig_window_size;
        CvSize real_window_size;
        double scale;
        CvHaarStageClassifier* stage_classifier;
        CvHidHaarClassifierCascade* hid_cascade;
    }
    CvHaarClassifierCascade;

    /* cascade or tree of stage classifiers */
    int  flags; /* signature */
    int  count; /* number of stages */
    CvSize orig_window_size; /* original object size (the cascade is trained for) */
    /* these two parameters are set by cvSetImagesForHaarClassifierCascade */
    CvSize real_window_size; /* current object size */
    double scale; /* current scale */
    CvHaarStageClassifier* stage_classifier; /* array of stage classifiers */
    CvHidHaarClassifierCascade* hid_cascade; /* hidden optimized representation of the cascade, created by cvSetImagesForHaarClassifierCascade */
    所有的结构都代表一个级联boosted Haar分类器。级联有下面的等级结构:
        Cascade:
            Stage1:
                Classifier11:
                    Feature11
                Classifier12:
                    Feature12
                ...
            Stage2:
                Classifier21:
                    Feature21
                ...
            ...
    整个等级可以手工构建,也可以利用函数cvLoadHaarClassifierCascade从已有的磁盘文件或嵌入式基中导入。
    特征检测用到的函数:
    cvLoadHaarClassifierCascade
    从文件中装载训练好的级联分类器或者从OpenCV中嵌入的分类器数据库中导入
    CvHaarClassifierCascade* cvLoadHaarClassifierCascade(
                             const char* directory,
                             CvSize orig_window_size );
    directory :训练好的级联分类器的路径
    orig_window_size:级联分类器训练中采用的检测目标的尺寸。因为这个信息没有在级联分类器中存储,所有要单独指出。
    函数 cvLoadHaarClassifierCascade 用于从文件中装载训练好的利用海尔特征的级联分类器,或者从OpenCV中嵌入的分类器数据库中导入。分类器的训练可以应用函数haartraining(详细察看opencv/apps/haartraining)
    函数 已经过时了。现在的目标检测分类器通常存储在  XML 或 YAML 文件中,而不是通过路径导入。从文件中导入分类器,可以使用函数 cvLoad 。
    cvReleaseHaarClassifierCascade
    释放haar classifier cascade。
    void cvReleaseHaarClassifierCascade( CvHaarClassifierCascade** cascade );
    cascade :双指针类型指针指向要释放的cascade. 指针由函数声明。
    函数 cvReleaseHaarClassifierCascade 释放cascade的动态内存,其中cascade的动态内存或者是手工创建,或者通过函数 cvLoadHaarClassifierCascade 或 cvLoad分配。
    cvHaarDetectObjects
    检测图像中的目标
    typedef struct CvAvgComp
    {
    CvRect rect; /* bounding rectangle for the object (average rectangle of a group) */
    int neighbors; /* number of neighbor rectangles in the group */
    }
    CvAvgComp;
    CvSeq* cvHaarDetectObjects( const CvArr* image,
    CvHaarClassifierCascade* cascade,
                                CvMemStorage* storage,
                                double scale_factor=1.1,
                                int min_neighbors=3, int flags=0,
                                CvSize min_size=cvSize(0,0) );
    image 被检图像
    cascade harr 分类器级联的内部标识形式
    storage 用来存储检测到的一序列候选目标矩形框的内存区域。
    scale_factor 在前后两次相继的扫描中,搜索窗口的比例系数。例如1.1指将搜索窗口依次扩大10%。
    min_neighbors 构成检测目标的相邻矩形的最小个数(缺省-1)。如果组成检测目标的小矩形的个数和小于min_neighbors-1 都会被排除。如果min_neighbors 为 0, 则函数不做任何操作就返回所有的被检候选矩形框,这种设定值一般用在用户自定义对检测结果的组合程序上。
    flags 操作方式。当前唯一可以定义的操作方式是 CV_HAAR_DO_CANNY_PRUNING。如果被设定,函数利用Canny边缘检测器来排除一些边缘很少或者很多的图像区域,因为这样的区域一般不含被检目标。人脸检测中通过设定阈值使用了这种方法,并因此提高了检测速度。 
    min_size 检测窗口的最小尺寸。缺省的情况下被设为分类器训练时采用的样本尺寸(人脸检测中缺省大小是~20×20)。
    函数 cvHaarDetectObjects 使用针对某目标物体训练的级联分类器在图像中找到包含目标物体的矩形区域,并且将这些区域作为一序列的矩形框返回。函数以不同比例大小的扫描窗口对图像进行几次搜索(察看cvSetImagesForHaarClassifierCascade)。 每次都要对图像中的这些重叠区域利用cvRunHaarClassifierCascade进行检测。 有时候也会利用某些继承(heuristics)技术以减少分析的候选区域,例如利用 Canny 裁减 (prunning)方法。 函数在处理和收集到候选的方框(全部通过级联分类器各层的区域)之后,接着对这些区域进行组合并且返回一系列各个足够大的组合中的平均矩形。调节程序中的缺省参数(scale_factor=1.1, min_neighbors=3, flags=0)用于对目标进行更精确同时也是耗时较长的进一步检测。为了能对视频图像进行更快的实时检测,参数设置通常是:scale_factor=1.2, min_neighbors=2, flags=CV_HAAR_DO_CANNY_PRUNING, min_size=<minimum possible face size> (例如, 对于视频会议的图像区域).
    cvSetImagesForHaarClassifierCascade
    为隐藏的cascade(hidden cascade)指定图像
    void cvSetImagesForHaarClassifierCascade( CvHaarClassifierCascade* cascade,
                                              const CvArr* sum, const CvArr* sqsum,
                                              const CvArr* tilted_sum, double scale );
    cascade 隐藏 Harr 分类器级联 (Hidden Haar classifier cascade), 由函数 cvCreateHidHaarClassifierCascade生成
    sum 32-比特,单通道图像的积分图像(Integral (sum) 单通道 image of 32-比特 integer format). 这幅图像以及随后的两幅用于对快速特征的评价和亮度/对比度的归一化。 它们都可以利用函数 cvIntegral从8-比特或浮点数 单通道的输入图像中得到。
    sqsum 单通道64比特图像的平方和图像
    tilted_sum 单通道32比特整数格式的图像的倾斜和(Tilted sum)
    scale cascade的窗口比例. 如果 scale=1, 就只用原始窗口尺寸检测 (只检测同样尺寸大小的目标物体) - 原始窗口尺寸在函数cvLoadHaarClassifierCascade中定义 (在 "<default_face_cascade>"中缺省为24x24), 如果scale=2, 使用的窗口是上面的两倍 (在face cascade中缺省值是48x48 )。 这样尽管可以将检测速度提高四倍,但同时尺寸小于48x48的人脸将不能被检测到。
    函数 cvSetImagesForHaarClassifierCascade 为hidden classifier cascade 指定图像 and/or 窗口比例系数。 如果图像指针为空,会继续使用原来的图像(i.e. NULLs 意味这"不改变图像")。比例系数没有 "protection" 值,但是原来的值可以通过函数 cvGetHaarClassifierCascadeScale 重新得到并使用。这个函数用于对特定图像中检测特定目标尺寸的cascade分类器的设定。函数通过cvHaarDetectObjects进行内部调用,但当需要在更低一层的函数cvRunHaarClassifierCascade中使用的时候,用户也可以自行调用。
    cvRunHaarClassifierCascade
    在给定位置的图像中运行 cascade of boosted classifier
    int cvRunHaarClassifierCascade( CvHaarClassifierCascade* cascade,
                                    CvPoint pt, int start_stage=0 );
    cascade Haar 级联分类器
    pt 待检测区域的左上角坐标。待检测区域大小为原始窗口尺寸乘以当前设定的比例系数。当前窗口尺寸可以通过cvGetHaarClassifierCascadeWindowSize重新得到。
    start_stage 级联层的初始下标值(从0开始计数)。函数假定前面所有每层的分类器都已通过。这个特征通过函数cvHaarDetectObjects内部调用,用于更好的处理器高速缓冲存储器。
    函数 cvRunHaarHaarClassifierCascade 用于对单幅图片的检测。在函数调用前首先利用 cvSetImagesForHaarClassifierCascade设定积分图和合适的比例系数 (=> 窗口尺寸)。当分析的矩形框全部通过级联分类器每一层的时返回正值(这是一个候选目标),否则返回0或负值。
    二、例程分析:
    例子:利用级联的Haar classifiers寻找检测目标(e.g. faces).
    #include "cv.h"
    #include "highgui.h"
    //读取训练好的分类器。
    CvHaarClassifierCascade* load_object_detector( const char* cascade_path )
    {
        return (CvHaarClassifierCascade*)cvLoad( cascade_path );
    }

    void detect_and_draw_objects( IplImage* image,
                                  CvHaarClassifierCascade* cascade,
                                  int do_pyramids )
    {
        IplImage* small_image = image;
        CvMemStorage* storage = cvCreateMemStorage(0); //创建动态内存
        CvSeq* faces;
        int i, scale = 1;
        /* if the flag is specified, down-scale the 输入图像 to get a
           performance boost w/o loosing quality (perhaps) */
        if( do_pyramids )
        {
            small_image = cvCreateImage( cvSize(image->width/2,image->height/2), IPL_DEPTH_8U, 3 );
            cvPyrDown( image, small_image, CV_GAUSSIAN_5x5 );//函数 cvPyrDown 使用 Gaussian 金字塔分解对输入图像向下采样。首先它对输入图像用指定滤波器进行卷积,然后通过拒绝偶数的行与列来下采样图像。
            scale = 2;
        }
        /* use the fastest variant */
        faces = cvHaarDetectObjects( small_image, cascade, storage, 1.2, 2, CV_HAAR_DO_CANNY_PRUNING );
        /* draw all the rectangles */
        for( i = 0; i < faces->total; i++ )
        {
            /* extract the rectanlges only */
            CvRect face_rect = *(CvRect*)cvGetSeqElem( faces, i, 0 );
            cvRectangle( image, cvPoint(face_rect.x*scale,face_rect.y*scale),
                         cvPoint((face_rect.x+face_rect.width)*scale,
                                 (face_rect.y+face_rect.height)*scale),
                         CV_RGB(255,0,0), 3 );
        }
        if( small_image != image )
            cvReleaseImage( &small_image );
        cvReleaseMemStorage( &storage );  //释放动态内存
    }
    /* takes image filename and cascade path from the command line */
    int main( int argc, char** argv )
    {
        IplImage* image;
        if( argc==3 && (image = cvLoadImage( argv[1], 1 )) != 0 )
        {
            CvHaarClassifierCascade* cascade = load_object_detector(argv[2]);
            detect_and_draw_objects( image, cascade, 1 );
            cvNamedWindow( "test", 0 );
            cvShowImage( "test", image );
            cvWaitKey(0);
            cvReleaseHaarClassifierCascade( &cascade );
            cvReleaseImage( &image );
        }
        return 0;
    }
    关键代码很简单,装载分类器,对输入图像进行金字塔采样,然后用cv的函数进行检测目标,最后输出检测到的目标矩形。

    OpenCV学习笔记(二)基于Haar-like特征的层叠推进分类器快速目标检测


    OpenCV学习笔记之二――基于Haar-like特征的层叠推进分类器快速目标检测
    一、简介
    目标检测方法最初由Paul Viola [Viola01]提出,并由Rainer Lienhart [Lienhart02]对这一方法进行了改善。该方法的基本步骤为: 首先,利用样本(大约几百幅样本图片)的 harr 特征进行分类器训练,得到一个级联的boosted分类器。
    分类器中的"级联"是指最终的分类器是由几个简单分类器级联组成。在图像检测中,被检窗口依次通过每一级分类器, 这样在前面几层的检测中大部分的候选区域就被排除了,全部通过每一级分类器检测的区域即为目标区域。
    分类器训练完以后,就可以应用于输入图像中的感兴趣区域(与训练样本相同的尺寸)的检测。检测到目标区域(汽车或人脸)分类器输出为1,否则输出为0。为了检测整副图像,可以在图像中移动搜索窗口,检测每一个位置来确定可能的目标。 为了搜索不同大小的目标物体,分类器被设计为可以进行尺寸改变,这样比改变待检图像的尺寸大小更为有效。所以,为了在图像中检测未知大小的目标物体,扫描程序通常需要用不同比例大小的搜索窗口对图片进行几次扫描。
    目前支持这种分类器的boosting技术有四种: Discrete Adaboost, Real Adaboost, Gentle Adaboost and Logitboost。
    "boosted" 即指级联分类器的每一层都可以从中选取一个boosting算法(权重投票),并利用基础分类器的自我训练得到。
    根据上面的分析,目标检测分为三个步骤:
    1、 样本的创建
    2、 训练分类器
    3、 利用训练好的分类器进行目标检测。
    二、样本创建
    训练样本分为正例样本和反例样本,其中正例样本是指待检目标样本(例如人脸或汽车等),反例样本指其它任意图片,所有的样本图片都被归一化为同样的尺寸大小(例如,20x20)。
    负样本
    负样本可以来自于任意的图片,但这些图片不能包含目标特征。负样本由背景描述文件来描述。背景描述文件是一个文本文件,每一行包含了一个负样本图片的文件名(基于描述文件的相对路径)。该文件必须手工创建。
    e.g: 负样本描述文件的一个例子:
    假定目录结构如下:
    /img
    img1.jpg
    img2.jpg
    bg.txt
    则背景描述文件bg.txt的内容为:
    img/img1.jpg
    img/img2.jpg
    正样本
    正样本由程序craatesample程序来创建。该程序的源代码由OpenCV给出,并且在bin目录下包含了这个可执行的程序。
    正样本可以由单个的目标图片或者一系列的事先标记好的图片来创建。
    Createsamples程序的命令行参数:
    命令行参数:
    -vec <vec_file_name>
    训练好的正样本的输出文件名。
    -img<image_file_name>
    源目标图片(例如:一个公司图标)
    -bg<background_file_name>
    背景描述文件。
    -num<number_of_samples>
    要产生的正样本的数量,和正样本图片数目相同。
    -bgcolor<background_color>
    背景色(假定当前图片为灰度图)。背景色制定了透明色。对于压缩图片,颜色方差量由bgthresh参数来指定。则在bgcolor-bgthresh和bgcolor+bgthresh中间的像素被认为是透明的。
    -bgthresh<background_color_threshold>
    -inv
    如果指定,颜色会反色
    -randinv
    如果指定,颜色会任意反色
    -maxidev<max_intensity_deviation>
    背景色最大的偏离度。
    -maxangel<max_x_rotation_angle>
    -maxangle<max_y_rotation_angle>,
    -maxzangle<max_x_rotation_angle>
    最大旋转角度,以弧度为单位。
    -show
    如果指定,每个样本会被显示出来,按下"esc"会关闭这一开关,即不显示样本图片,而创建过程继续。这是个有用的debug选项。
    -w<sample_width>
    输出样本的宽度(以像素为单位)
    -h《sample_height》
    输出样本的高度,以像素为单位。
    注:正样本也可以从一个预先标记好的图像集合中获取。这个集合由一个文本文件来描述,类似于背景描述文件。每一个文本行对应一个图片。每行的第一个元素是图片文件名,第二个元素是对象实体的个数。后面紧跟着的是与之匹配的矩形框(x, y, 宽度,高度)。
    下面是一个创建样本的例子:
    假定我们要进行人脸的检测,有5个正样本图片文件img1.bmp,…img5.bmp;有2个背景图片文件:bg1.bmp,bg2.bmp,文件目录结构如下:
    positive
      img1.bmp
      ……
      Img5.bmp
    negative
      bg1.bmp
      bg2.bmp
    info.dat
    bg.txt
    正样本描述文件info.dat的内容如下:
    Positive/imag1.bmp 1 0 0 24 28
    ……
    Positive/imag5.bmp 1 0 0 24 28
    图片img1.bmp包含了单个目标对象实体,矩形为(0,0,24,28)。
    注意:要从图片集中创建正样本,要用-info参数而不是用-img参数。
    -info <collect_file_name>
    标记特征的图片集合的描述文件。
    背景(负样本)描述文件的内容如下:
    nagative/bg1.bmp
    nagative/bg2.bmp
    我们用一个批处理文件run.bat来进行正样本的创建:该文件的内容如下:
    cd  e:/face/bin
    CreateSamples   -vec e:/face/a.vec
     -info e:/face/info.dat
    -bg e:/face/bg.txt
    -num 5
    -show
    -w 24
     -h 28
    其中e:/face/bin目录包含了createsamples可执行程序,生成的正样本文件a.vec在e:/face目录下。
    三、训练分类器
    样本创建之后,接下来要训练分类器,这个过程是由haartraining程序来实现的。该程序源码由OpenCV自带,且可执行程序在OpenCV安装目录的bin目录下。
    Haartraining的命令行参数如下:
    -data<dir_name>
    存放训练好的分类器的路径名。
    -vec<vec_file_name>
    正样本文件名(由trainingssamples程序或者由其他的方法创建的)
    -bg<background_file_name>
    背景描述文件。
    -npos<number_of_positive_samples>,
    -nneg<number_of_negative_samples>
    用来训练每一个分类器阶段的正/负样本。合理的值是:nPos = 7000;nNeg = 3000
    -nstages<number_of_stages>
    训练的阶段数。
    -nsplits<number_of_splits>
    决定用于阶段分类器的弱分类器。如果1,则一个简单的stump classifier被使用。如果是2或者更多,则带有number_of_splits个内部节点的CART分类器被使用。
    -mem<memory_in_MB>
    预先计算的以MB为单位的可用内存。内存越大则训练的速度越快。
    -sym(default)
    -nonsym
    指定训练的目标对象是否垂直对称。垂直对称提高目标的训练速度。例如,正面部是垂直对称的。
    -minhitrate《min_hit_rate》
    每个阶段分类器需要的最小的命中率。总的命中率为min_hit_rate的number_of_stages次方。
    -maxfalsealarm<max_false_alarm_rate>
    没有阶段分类器的最大错误报警率。总的错误警告率为max_false_alarm_rate的number_of_stages次方。
    -weighttrimming<weight_trimming>
    指定是否使用权修正和使用多大的权修正。一个基本的选择是0.9
    -eqw
    -mode<basic(default)|core|all>
    选择用来训练的haar特征集的种类。basic仅仅使用垂直特征。all使用垂直和45度角旋转特征。
    -w《sample_width》
    -h《sample_height》
    训练样本的尺寸,(以像素为单位)。必须和训练样本创建的尺寸相同。
    一个训练分类器的例子:
    同上例,分类器训练的过程用一个批处理文件run2.bat来完成:
    cd e:/face/bin
    haartraining -data e:/face/data
    -vec e:/face/a.vec
    -bg e:/face/bg.txt
    -npos 5
    -nneg 2
     -w 24
     -h 28
    训练结束后,会在目录data下生成一些子目录,即为训练好的分类器。
    注:OpenCv的某些版本可以将这些目录中的分类器直接转换成xml文件。但在实际的操作中,haartraining程序却好像永远不会停止,而且没有生成xml文件,后来在OpenCV的yahoo论坛上找到一个haarconv的程序,才将分类器转换为xml文件,其中的原因尚待研究。
    四、目标检测
    OpenCV的cvHaarDetectObjects()函数(在haarFaceDetect演示程序中示例)被用来做侦测。关于该检测的详细分析,将在下面的笔记中详细描述。

     
    3月29日

    一个合格的程序员该做的事情


    2006.03.15  来自:CSDN 选自Mailbomb 的blog  
     

      程序员每天该做的事

    1、总结自己一天任务的完成情况

    最好的方式是写工作日志,把自己今天完成了什么事情,遇见了什么问题都记录下来,日后翻看好处多多

    2、考虑自己明天应该做的主要工作

    把明天要做的事情列出来,并按照优先级排列,第二天应该把自己效率最高的时间分配给最重要的工作

    3、考虑自己一天工作中失误的地方,并想出避免下一次再犯的方法

    出错不要紧,最重要的是不要重复犯相同的错误,那是愚蠢

    4、考虑自己一天工作完成的质量和效率能否还能提高

    一天只提高1%,365天你的效率就能提高多少倍你知道吗? (1+0.01)^365 = 37 倍

    5、看一个有用的新闻网站或读一张有用的报纸,了解业界动态

    闭门造车是不行的,了解一下别人都在做什么,对自己能带来很多启示

    6、记住一位同事的名字及其特点

    你认识公司的所有同事吗?你了解他们吗?

    7、清理自己的代码

    今天完成的代码,把中间的调试信息,测试代码清理掉,按照编码风格整理好,注释都写好了吗?

    8、清理自己的桌面

    当日事当日毕,保持清洁干劲的桌面才能让你工作时不分心,程序员特别要把电脑的桌面清理干净

    程序员每月该做的事

    1、至少和一个同事一起吃饭或喝茶
    不光了解自己工作伙伴的工作,还要了解他们的生活

    2、自我考核一次
    相对正式地考核自己一下,你对得起这个月的工资吗?

    3、对你的同事考核一次
    你的同事表现怎么样?哪些人值得学习,哪些人需要帮助?

    3、制定下月的计划,确定下月的工作重点

    4、总结自己工作质量改进状况
    自己的质量提高了多少?

    5、有针对性地对一项工作指标做深入地分析并得出改进的方案
    可以是对自己的,也可以是对公司的,一定要深入地分析后拿出自己的观点来。要想在老板面前说得上话,做的成事,工作上功夫要做足。

    6、与老板沟通一次
    最好是面对面地沟通,好好表现一下自己,虚心听取老板的意见,更重要的是要了解老板当前关心的重点

    程序员每年该做的事

    1、年终总结
    每个公司都会做的事情,但你真正认真地总结过自己吗?

    2、兑现给自己、给家人的承诺
    给老婆、儿子的新年礼物买了没有?给自己的呢?

    3、下年度工作规划
    好好想想自己明年的发展目标,争取升职/加薪、跳槽还是自己出来干?

    4、掌握一项新技术
    至少是一项,作为程序员一年要是一项新技术都学不到手,那就一定会被淘汰。
    掌握可不是看本书就行的,要真正懂得应用,最好你能够写一篇教程发表到你的blog

    5、推出一种新产品
    可以是一个真正的产品,也可以只是一个类库,只要是你创造的东西就行,让别人使用它,也为世界作点贡献。当然如果真的很有价值,收点注册费也是应该的

    6、与父母团聚一次
    常回家看看,常回家看看

     

    3月28日

    OpenCV学习笔记(一)概述和系统配置


    OpenCV学习笔记(一)概述和系统配置

    一、概述
    OpenCV是英特尔公司于1999年在俄罗斯设立的软件开发中心“Software Development Center”开发的。该公司一直致力于基于个人电脑的计算机视觉应用的开发,可以实时追踪的视觉用户接口技术的普及为目标。初步拟定应用于Human-Computer Interaction(HCI,人机互动)、物体确定、面孔识别、表情识别,移动物体追踪、自主运动(Ego-motion)、移动机器人等领域。因此,“OpenCV 2.1将提供给玩具制造商及机器人制造商等从事计算机视觉相关技术的各类企业/团体”
        OpenCV将以公开源码的方式提供,也就是接受方有权在修改之后另行向第三方提供。源代码(C语言)中包括有库(Library)的所有功能。详细情况刊登在OpenCV的WWW站点上。
    英特尔公司解释说,“速度更高的微处理器、廉价的数码相机以及USB 2等技术使高速视频捕获(Video Capture)成为可能,因此,基于普通个人电脑的实时计算机视觉将有望实现”。
    OpenCV的最新版本为beta5。
    二、OpenCV组成部分
    目前OpenCV包含下面几个部分:
        cvcore:一些基本函数(各种数据类型的基本运算等)
        cv:    图像处理和计算机视觉功能(图像处理,结构分析,运动分析,物体跟踪,模式识别,摄像机定标)
        cvaux:一些实验性的函数(view morphing,三维跟踪,pca,hmm)
        highgui 用户交互部分(GUI,图像视频I/O,系统调用函数)
        另外还有cvcam,不过linux版本已经抛弃。windows版本中将directX支持加入highgui后,cvcam将被彻底去掉。
    三、OpenCV特点
    OpenCV是Intel公司开发的图像处理和计算机视觉函数库,它有以下特点:
    1) 开放C源码
    2) 基于Intel处理器指令集开发的优化代码
    3) 统一的结构和功能定义
    4) 强大的图像和矩阵运算能力
    5) 方便灵活的用户接口
    6)同时支持MS-WINDOWS、LINUX平台
    四、下载OpenCV
        http://www.sourceforge.net/projects/opencvlibrary
    五、参考资料
        =》源代码及文档下载:SOURCEFORGE.NET
    http://sourceforge.net/projects/opencvlibrary/
        =》INTEL的OPENCV主页:
    http://www.intel.com/research/mrl/research/opencv/
        =》YAHOO OPENCV 的邮件列表:
    http://groups.yahoo.com/group/OpenCV/
        =》CMU(卡耐基-梅隆大学)的计算机视觉主页:
    http://www-2.cs.cmu.edu/afs/cs/project/cil/ftp/html/vision.html
        =》OPENCV 更为详细的介绍
    http://www.assuredigit.com//incoming/sourcecode/opencv/chinese_docs/index.htm
        =》OPENCV 的常用问题与解答
    http://www.assuredigit.com//incoming/sourcecode/opencv/chinese_docs/faq.htm
        =》OPENCV 的安装指南
    http://www.assuredigit.com//incoming/sourcecode/opencv/chinese_docs/install
        =》更多的最新资料,请访问
    http://blog.csdn.net/hunnishhttp://rocee.bokee.com
    六、创建一个 DeveloperStudio 项目来开始 OpenCV
    1. 在 Developer Studio 中创建新的应用程序:
    选择菜单 "File"->"New..."->"Projects" . 选择 "Win32 Application" 或 "Win32 console application" - 后者是更简单的方法。
    键入项目名称,并且选择存储位置
    可以为项目创建一个单独的 workspace ("Create new workspace") , 也可以将新的项目加入到当前的 workspace 中 ("Add to current workspace").
    单击 "next" 
    选择 "An empty project", 点击 "Finish", "OK".
    经过以上步骤,Developer Studio 会创建一个项目目录 (缺省情况下,目录名就是项目名), .dsp 文件以及.dsw,.ncb ... ,如果你创建自己的workspace.
    2添加文件到 project 中:
    选择菜单"File"->"New..."->"Files" .
    选择"C++ Source File", 键入文件名,点击"OK"
    增加 OpenCV 相关的 头文件目录 #include :
            #include "cv.h"
            /* #inlcude "cvaux.h" // experimental stuff (if need) */
            #include "highgui.h"
         
    或者你可以拷贝部分已有的文件 (如:opencv/samples/c/morphology.c) 到项目目录中,打开它,并且加入到项目中 (右键点击编辑器的视图 -> "Insert File into Project" -> ).
    3配置项目:
    选择菜单"Project"->"Settings..."以激活项目配置对话框 .
    在左边选择你的项目.
    调节设置,对 Release 和 Debug 配置都有效:
    选择 "Settings For:"->"All Configurations"
    选择 "C/C++" tab -> "Preprocessor" category -> "Additional Include Directories:". 加入用逗号分隔的相对路径 (对文件 .dsp 而言) 或绝对路径
     d:/opencv/cxcore/include,d:/opencv/cv/include,d:/opencv/otherlibs/highgui, d:/opencv/cvaux/include(optionally,)
    选择 "Link" tab -> "Input" category -> "Additional library path:".
    加入输入库所在的路径 (cxcore[d].lib cv[d].lib hihghui[d].lib cvaux[d].lib)
    d:/opencv/lib
    调节 "Debug" 配置:
    选择 "Settings For:"->"Win32 Debug".
    选择 "Link" tab -> "General" category -> "Object/library modules". 加入空格分隔的 cvd.lib,cxcored.lib highguid.lib,cvauxd.lib (optionally)
    可以改变输出文件的名称和位置。如想把产生的 .exe 文件放置于项目目录而不是Debug/ 子目录下,可在 "Link" tab -> "General" category -> "Output file name:" 中键入 ./d.exe 
    调节 "Release" 配置
    选择 "Settings For:"->"Win32 Release".
    选择 "Link" tab -> "General" category -> "Object/library modules". 加入空格分隔的cv.lib cxcore.lib highgui.lib cvaux.lib (optionally)
    4增加从属性项目到 workspace 中:
    选择菜单: "Project" -> "Insert project into workspace".
    选择 opencv/cv/make/cv.dsp.
    同样步骤对 opencv/cvaux/make/cvaux.dsp, opencv/otherlibs/highgui/highgui.dsp.
    设置从属性:
    选择菜单: "Project" -> "Dependencies..."
    对 "cv" 选择 "cxcore",
    对 "cvaux" 选择 "cv", "cxcore",
    对 "highgui" 选择 "cxcore",
    对你的项目,选择所有的: "cxcore", "cv", "cvaux", "highgui".
    从属性配置保证了在源代码被改变的情况下,自动重新编译 opencv 库.
    5就这么多。可以编译并且运行一切了。
    七、库设置:
       静态库设置:
       Opencv程序需要静态库设置,其release版本的静态库在系统的lib目录下,其debug版本的静态库需要重新全编译所有的程序。
        动态库设置:
       OPenCV启动时需要一些动态库的支持,这些动态库必须放在系统目录下或者当前目录下。
        Cv097.dll,cvaux097.dll,cvcam097.dll,cxcore097.dll highguid097.dll,libguide40.dll


     

    出租司机给"我"上的MBA课

    Posted on 2006-03-15 01:14 刘润
    我要从徐家汇赶去机场,于是匆匆结束了一个会议,在美罗大厦前搜索出租车。一辆大众发现了我,非常专业的、径直的停在我的面前。这一停,于是有了后面的这个让我深感震撼的故事,象上了一堂生动的MBA案例课。为了忠实于这名出租车司机的原意,我凭记忆尽量重复他原来的话。
    “去哪里……好的,机场。我在徐家汇就喜欢做美罗大厦的生意。这里我只做两个地方。美罗大厦,均瑶大厦。你知道吗?接到你之前,我在美罗大厦门口兜了两圈,终于被我看到你了!从写字楼里出来的,肯定去的不近~~~”
    “哦?你很有方法嘛!”我附和了一下。
    “做出租车司机,也要用科学的方法。”他说。我一愣,顿时很有些兴趣“什么科学的方法?”
    “要懂得统计。我做过精确的计算。我说给你听啊。我每天开17个小时的车,每小时成本34.5元……”
    “怎么算出来的?”我追问。
    “你算啊,我每天要交380元,油费大概210元左右。一天17小时,平均每小时固定成本22元,交给公司,平均每小时12.5元油费。这是不是就是34.5元?”,我有些惊讶。我打了10年的车,第一次听到有出租车司机这么计算成本。以前的司机都和我说,每公里成本0.3元,另外每天交多少钱之类的。
    “成本是不能按公里算的,只能按时间算。你看,计价器有一个“检查”功能。你可以看到一天的详细记录。我做过数据分析,每次载客之间的空驶时间平均为7分钟。如果上来一个起步价,10元,大概要开10分钟。也就是每一个10元的客人要花17分钟的成本,就是9.8元。不赚钱啊!如果说做浦东、杭州、青浦的客人是吃饭,做10元的客人连吃菜都算不上,只能算是撒了些味精。”
    强!这位师傅听上去真不象出租车司机,到象是一位成本核算师。“那你怎么办呢?”我更感兴趣了,继续问。看来去机场的路上还能学到新东西。
    “千万不能被客户拉了满街跑。而是通过选择停车的地点,时间,和客户,主动地决定你要去的地方。”我非常惊讶,这听上去很有意思。“有人说做出租车司机是靠运气吃饭的职业。我以为不是。你要站在客户的位置上,从客户的角度去思考。”这句话听上去很专业,有点象很多商业管理培训老师说的“put yourself into others' shoes.”
    “给你举个例子,医院门口,一个拿着药的,一个拿着脸盆的,你带哪一个。”我想了想,说不知道。
    “你要带那个拿脸盆的。一般人病小痛的到医院看一看,拿点药,不一定会去很远的医院。拿着脸盆打车的,那是出院的。住院哪有不死人的?今天二楼的谁死了,明天三楼又死了一个。从医院出来的人通常会有一种重获新生的感觉,重新认识生命的意义,健康才最重要。那天这个说:走,去青浦。眼睛都不眨一下。你说他会打车到人民广场,再去做青浦线吗?绝对不会!”
    我不由得开始佩服。
    “再给你举个例子。那天人民广场,三个人在前面招手。一个年轻女子,拿着小包,刚买完东西。还有一对青年男女,一看就是逛街的。第三个是个里面穿绒衬衫的,外面羽绒服的男子,拿着笔记本包。我看一个人只要3秒钟。我毫不犹豫地停在这个男子面前。这个男的上车后说:延安高架、南北高架~~~还没说后面就忍不住问,为什么你毫不犹豫地开到我面前?前面还有两个人,他们要是想上车,我也不好意思和他们抢。我回答说,中午的时候,还有十几分钟就1点了。那个女孩子是中午溜出来买东西的,估计公司很近;那对男女是游客,没拿什么东西,不会去很远;你是出去办事的,拿着笔记本包,一看就是公务。而且这个时候出去,估计应该不会近。那个男的就说,你说对了,去宝山。”
    “那些在超市门口,地铁口打车,穿着睡衣的人可能去很远吗?可能去机场吗?机场也不会让她进啊。”
    有道理!我越听越有意思。
    “很多司机都抱怨,生意不好做啊,油价又涨了啊,都从别人身上找原因。我说,你永远从别人身上找原因,你永远不能提高。从自己身上找找看,问题出在哪里。”这话听起来好熟,好像是“如果你不能改变世界,就改变你自己”,或者Steven Corvey的“影响圈和关注圈”的翻版。“有一次,在南丹路一个人拦车,去田林。后来又有一次,一个人在南丹路拦车,还是去田林。我就问了,怎么你们从南丹路出来的人,很多都是去田林呢?人家说,在南丹路有一个公共汽车总站,我们都是坐公共汽车从浦东到这里,然后搭车去田林的。我恍然大悟。比如你看我们开过的这条路,没有写字楼,没有酒店,什么都没有,只有公共汽车站,站在这里拦车的多半都是刚下公共汽车的,再选择一条最短路经打车。在这里拦车的客户通常不会高于15元。”
    “所以我说,态度决定一切!”我听十几个总裁讲过这句话,第一次听出租车司机这么说。
    “要用科学的方法,统计学来做生意。天天等在地铁站口排队,怎么能赚到钱?每个月就赚500块钱怎么养活老婆孩子?这就是在谋杀啊!慢性谋杀你的全家。要用知识武装自己。学习知识可以把一个人变成聪明的人,一个聪明的人学习知识可以变成很聪明的人。一个很聪明的人学习知识,可以变成天才。”
    “有一次一个人打车去火车站,问怎么走。他说这么这么走。我说慢,上高架,再这么这么走。他说,这就绕远了。我说,没关系,你经常走你有经验,你那么走50块,你按我的走法,等里程表50块了,我就翻表。你只给50快就好了,多的算我的。按你说的那么走要50分钟,我带你这么走只要25分钟。最后,按我的路走,多走了4公里,快了25分钟,我只收了50块。乘客很高兴,省了10元钱左右。这4公里对我来说就是1块多钱的油钱。我相当于用1元多钱买了25分钟。我刚才说了,我一小时的成本34.5块,我多合算啊!”
    “在大众公司,一般一个司机3、4千,拿回家。做的好的大概5千左右。顶级的司机大概每月能有7000。全大众2万个司机,大概只有2-3个司机,万里挑一,每月能拿到8000以上。我就是这2-3个人中间的一个。而且很稳定,基本不会大的波动。”
    太强了!到此为止,我越来越佩服这个出租车司机。
    “我常常说我是一个快乐的车夫。有人说,你是因为赚的钱多,所以当然快乐。我对他们说,你们正好错了。是因为我有快乐、积极的心态,所以赚的钱多。”
    说的多好啊!
    “要懂得体味工作带给你的美。堵在人民广场的时候,很多司机抱怨,又堵车了!真是倒霉。千万不要这样,用心体会一下这个城市的美,外面有很多漂亮的女孩子经过,非常现代的高楼大厦,虽然买不起,但是却可以用欣赏的眼光去享受。开车去机场,看着两边的绿色,冬天是白色的,多美啊。再看看里程表,100多了,就更美了!每一样工作都有她美丽的地方,我们要懂得从工作中体会这种美丽。”
    “我10年前是强生公司的总教练。8年前在公司作过三个不同部门的部门经理。后来我不干了,一个月就3、5千块,没意思。就主动来做司机。我愿意做一个快乐的车夫。哈哈哈哈。”
    到了机场,我给他留了一张名片,说:“你有没有兴趣这个星期五,到我办公室,给微软的员工讲一讲你怎么开出租车的?你就当打着表,60公里一小时,你讲多久,我就付你多少钱。给我电话。”
    我迫不及待的在飞机上记录下他这堂生动的MBA课。
    展开全文
  • opencv stitching算法分析

    千次阅读 2017-04-17 17:10:33
    #ifdef HAVE_OPENCV_GPU     if  (try_use_gpu && getCudaEnabledDeviceCount() > 0)   impl_ =  new  GpuMatcher(match_conf);    else   #else    ( void )try_use_gpu;  #...

    原文地址:http://blog.csdn.net/skeeee/article/details/19480693?utm_source=tuicool

    一、stitching_detail程序运行流程

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

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

          3.对图像的特征点进行匹配,使用最近邻和次近邻方法,将两个最优的匹配的置信度保存下来。

          4.对图像进行排序以及将置信度高的图像保存到同一个集合中,删除置信度比较低的图像间的匹配,得到能正确匹配的图像序列。这样将置信度高于门限的所有匹配合并到一个集合中。

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

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

         7.波形校正,水平或者垂直

         8.拼接

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


    二、stitching_detail程序接口介绍

        img1 img2 img3 输入图像

         --preview  以预览模式运行程序,比正常模式要快,但输出图像分辨率低,拼接的分辨率compose_megapix 设置为0.6

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

    /* 运动估计参数 */

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

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

        --match_conf <float> 特征点检测置信等级,最近邻匹配距离与次近邻匹配距离的比值,surf默认为0.65,orb默认为0.3

        --conf_thresh <float> 两幅图来自同一全景图的置信度,默认为1.0

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

        --ba_refine_mask (mask)  ---------------

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

         --save_graph <file_name> 将匹配的图形以点的形式保存到文件中, Nm代表匹配的数量,NI代表正确匹配的数量,C表示置信度

    /*图像融合参数:*/

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

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

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

        --seam_megapix <float> 拼接缝像素的大小 默认是0.1 ------------

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

        --compose_megapix <float> 拼接分辨率,默认为-1

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

        --blend (no|feather|multiband) 融合方法,默认是多频段融合

        --blend_strength <float> 融合强度,0-100.默认是5.

       --output <result_img> 输出图像的文件名,默认是result,jpg

        命令使用实例,以及程序运行时的提示:


    三、程序代码分析

        1.参数读入

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

    [cpp]  view plain copy
    1. int retval = parseCmdArgs(argc, argv);  
    2. if (retval)  
    3.     return retval;  
    4.   
    5. // Check if have enough images  
    6. int num_images = static_cast<int>(img_names.size());  
    7. if (num_images < 2)  
    8. {  
    9.     LOGLN("Need more images");  
    10.     return -1;  
    11. }  

          2.特征点检测

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

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

    [cpp]  view plain copy
    1. work_scale = min(1.0, sqrt(work_megapix * 1e6 / full_img.size().area()));  
    [cpp]  view plain copy
    1. resize(full_img, img, Size(), work_scale, work_scale);  
          图像大小是740*830,面积大于6*10^5,所以计算出work_scale = 0.98,然后对图像resize。 

         

         2.2 计算seam_scale,也是根据图像的面积小于seam_megapix*10^6,(seam_megapix 默认是0.1),seam_work_aspect目前还没用到

    [cpp]  view plain copy
    1. seam_scale = min(1.0, sqrt(seam_megapix * 1e6 / full_img.size().area()));  
    2. seam_work_aspect = seam_scale / work_scale; //seam_megapix = 0.1 seam_work_aspect = 0.69  
         
        2.3 计算图像特征点,以及计算特征点描述子,并将img_idx设置为i。

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

    [cpp]  view plain copy
    1. struct detail::ImageFeatures  
    2. Structure containing image keypoints and descriptors.  
    3. struct CV_EXPORTS ImageFeatures  
    4. {  
    5.     int img_idx;//   
    6.     Size img_size;  
    7.     std::vector<KeyPoint> keypoints;  
    8.     Mat descriptors;  
    9. };  

        
         2.4 将源图像resize到seam_megapix*10^6,并存入image[]中

    [cpp]  view plain copy
    1. resize(full_img, img, Size(), seam_scale, seam_scale);  
    2. images[i] = img.clone();  
        3.图像匹配

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

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

    [cpp]  view plain copy
    1. vector<MatchesInfo> pairwise_matches;//Structure containing information about matches between two images.   
    2. BestOf2NearestMatcher matcher(try_gpu, match_conf);//最近邻和次近邻法  
    3. matcher(features, pairwise_matches);//对每两个图片进行matcher,20-》400 matchers.cpp 502  
        介绍一下BestOf2NearestMatcher 函数:

    [cpp]  view plain copy
    1.    //Features matcher which finds two best matches for each feature and leaves the best one only if the ratio between descriptor distances is greater than the threshold match_conf.  
    2.   detail::BestOf2NearestMatcher::BestOf2NearestMatcher(bool try_use_gpu=false,float match_conf=0.3f,  
    3.                                                    intnum_matches_thresh1=6, int num_matches_thresh2=6)  
    4.   Parameters:   try_use_gpu – Should try to use GPU or not  
    5. match_conf – Match distances ration threshold  
    6. num_matches_thresh1 – Minimum number of matches required for the 2D projective  
    7. transform estimation used in the inliers classification step  
    8. num_matches_thresh2 – Minimum number of matches required for the 2D projective  
    9. transform re-estimation on inliers  
         函数的定义(只是设置一下参数,属于构造函数):

    [cpp]  view plain copy
    1. BestOf2NearestMatcher::BestOf2NearestMatcher(bool try_use_gpu, float match_conf, int num_matches_thresh1, int num_matches_thresh2)  
    2. {  
    3. #ifdef HAVE_OPENCV_GPU  
    4.     if (try_use_gpu && getCudaEnabledDeviceCount() > 0)  
    5.         impl_ = new GpuMatcher(match_conf);  
    6.     else  
    7. #else  
    8.     (void)try_use_gpu;  
    9. #endif  
    10.         impl_ = new CpuMatcher(match_conf);  
    11.   
    12.     is_thread_safe_ = impl_->isThreadSafe();  
    13.     num_matches_thresh1_ = num_matches_thresh1;  
    14.     num_matches_thresh2_ = num_matches_thresh2;  
    15. }  

         以及MatchesInfo的结构体定义:

    [cpp]  view plain copy
    1. Structure containing information about matches between two images. It’s assumed that there is a homography between those images.  
    2.     struct CV_EXPORTS MatchesInfo  
    3.     {  
    4.         MatchesInfo();  
    5.         MatchesInfo(const MatchesInfo &other);  
    6.         const MatchesInfo& operator =(const MatchesInfo &other);  
    7.         int src_img_idx, dst_img_idx; // Images indices (optional)  
    8.         std::vector<DMatch> matches;  
    9.         std::vector<uchar> inliers_mask; // Geometrically consistent matches mask  
    10.         int num_inliers; // Number of geometrically consistent matches  
    11.         Mat H; // Estimated homography  
    12.         double confidence; // Confidence two images are from the same panorama  
    13.     };  
          求出图像匹配的结果如下(具体匹配参见sift特征点匹配),任意两幅图都进行匹配(3*3=9),其中1-》2和2-》1只计算了一次,以1-》2为准,:

         

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

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

    [cpp]  view plain copy
    1. "outimage101.jpg" -- "outimage102.jpg"[label="Nm=866, Ni=637, C=2.37864"];  
    2. "outimage101.jpg" -- "outimage103.jpg"[label="Nm=165, Ni=59, C=1.02609"];  
    3. "outimage102.jpg" -- "outimage103.jpg"[label="Nm=460, Ni=260, C=1.78082"];  
          Nm代表匹配的数量,NI代表正确匹配的数量,C表示置信度

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

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

    [cpp]  view plain copy
    1. vector<int> indices = leaveBiggestComponent(features, pairwise_matches, conf_thresh);//将置信度高于门限的所有匹配合并到一个集合中  
    2. vector<Mat> img_subset;  
    3. vector<string> img_names_subset;  
    4. vector<Size> full_img_sizes_subset;  
    5. for (size_t i = 0; i < indices.size(); ++i)  
    6. {  
    7.     img_names_subset.push_back(img_names[indices[i]]);  
    8.     img_subset.push_back(images[indices[i]]);  
    9.     full_img_sizes_subset.push_back(full_img_sizes[indices[i]]);  
    10. }  
    11.   
    12. images = img_subset;  
    13. img_names = img_names_subset;  
    14. full_img_sizes = full_img_sizes_subset;  
           4.根据单应性矩阵粗略估计出相机参数(焦距)

            4.1 焦距参数的估计

            根据前面求出的任意两幅图的匹配,我们根据两幅图的单应性矩阵H,求出符合条件的f,(4副图,16个匹配,求出8个符合条件的f),然后求出f的均值或者中值当成所有图形的粗略估计的f。

    [cpp]  view plain copy
    1. estimateFocal(features, pairwise_matches, focals);  
           函数的主要源码如下:

    [cpp]  view plain copy
    1.  for (int i = 0; i < num_images; ++i)  
    2.  {  
    3.      for (int j = 0; j < num_images; ++j)  
    4.      {  
    5.          const MatchesInfo &m = pairwise_matches[i*num_images + j];  
    6.          if (m.H.empty())  
    7.              continue;  
    8.          double f0, f1;  
    9.          bool f0ok, f1ok;  
    10. focalsFromHomography(m.H, f0, f1, f0ok, f1ok);//Tries to estimate focal lengths from the given homography  79  
    11. //under the assumption that the camera undergoes rotations around its centre only.  
    12.          if (f0ok && f1ok)  
    13.              all_focals.push_back(sqrt(f0 * f1));  
    14.      }  
    15.  }  

          其中函数focalsFromHomography的定义如下:

    [cpp]  view plain copy
    1. Tries to estimate focal lengths from the given homography  
    2.     under the assumption that the camera undergoes rotations around its centre only.    
    3.     Parameters  
    4.     H – Homography.  
    5.     f0 – Estimated focal length along X axis.  
    6.     f1 – Estimated focal length along Y axis.  
    7.     f0_ok – True, if f0 was estimated successfully, false otherwise.  
    8.     f1_ok – True, if f1 was estimated successfully, false otherwise.  
         函数的源码:

    [cpp]  view plain copy
    1. void focalsFromHomography(const Mat& H, double &f0, double &f1, bool &f0_ok, bool &f1_ok)  
    2. {  
    3.     CV_Assert(H.type() == CV_64F && H.size() == Size(3, 3));//Checks a condition at runtime and throws exception if it fails  
    4.   
    5.     const double* h = reinterpret_cast<const double*>(H.data);//强制类型转换  
    6.   
    7.     double d1, d2; // Denominators  
    8.     double v1, v2; // Focal squares value candidates  
    9.   
    10.     //具体的计算过程有点看不懂啊  
    11.     f1_ok = true;  
    12.     d1 = h[6] * h[7];  
    13.     d2 = (h[7] - h[6]) * (h[7] + h[6]);  
    14.     v1 = -(h[0] * h[1] + h[3] * h[4]) / d1;  
    15.     v2 = (h[0] * h[0] + h[3] * h[3] - h[1] * h[1] - h[4] * h[4]) / d2;  
    16.     if (v1 < v2) std::swap(v1, v2);  
    17.     if (v1 > 0 && v2 > 0) f1 = sqrt(std::abs(d1) > std::abs(d2) ? v1 : v2);  
    18.     else if (v1 > 0) f1 = sqrt(v1);  
    19.     else f1_ok = false;  
    20.   
    21.     f0_ok = true;  
    22.     d1 = h[0] * h[3] + h[1] * h[4];  
    23.     d2 = h[0] * h[0] + h[1] * h[1] - h[3] * h[3] - h[4] * h[4];  
    24.     v1 = -h[2] * h[5] / d1;  
    25.     v2 = (h[5] * h[5] - h[2] * h[2]) / d2;  
    26.     if (v1 < v2) std::swap(v1, v2);  
    27.     if (v1 > 0 && v2 > 0) f0 = sqrt(std::abs(d1) > std::abs(d2) ? v1 : v2);  
    28.     else if (v1 > 0) f0 = sqrt(v1);  
    29.     else f0_ok = false;  
    30. }  

          求出的焦距有8个

         

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

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

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

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

          camera其他参数:

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

          K矩阵的值:


    [cpp]  view plain copy
    1. Mat CameraParams::K() const  
    2. {  
    3.     Mat_<double> k = Mat::eye(3, 3, CV_64F);  
    4.     k(0,0) = focal; k(0,2) = ppx;  
    5.     k(1,1) = focal * aspect; k(1,2) = ppy;  
    6.     return k;  
    7. }  

          R矩阵的值:

         

    [cpp]  view plain copy
    1. void operator ()(const GraphEdge &edge)  
    2. {  
    3.     int pair_idx = edge.from * num_images + edge.to;  
    4.   
    5.     Mat_<double> K_from = Mat::eye(3, 3, CV_64F);  
    6.     K_from(0,0) = cameras[edge.from].focal;  
    7.     K_from(1,1) = cameras[edge.from].focal * cameras[edge.from].aspect;  
    8.     K_from(0,2) = cameras[edge.from].ppx;  
    9.     K_from(0,2) = cameras[edge.from].ppy;  
    10.   
    11.     Mat_<double> K_to = Mat::eye(3, 3, CV_64F);  
    12.     K_to(0,0) = cameras[edge.to].focal;  
    13.     K_to(1,1) = cameras[edge.to].focal * cameras[edge.to].aspect;  
    14.     K_to(0,2) = cameras[edge.to].ppx;  
    15.     K_to(0,2) = cameras[edge.to].ppy;  
    16.   
    17.     Mat R = K_from.inv() * pairwise_matches[pair_idx].H.inv() * K_to;  
    18.     cameras[edge.to].R = cameras[edge.from].R * R;  
    19. }  
             光轴中心的值

    [cpp]  view plain copy
    1. for (int i = 0; i < num_images; ++i)  
    2. {  
    3.     cameras[i].ppx += 0.5 * features[i].img_size.width;  
    4.     cameras[i].ppy += 0.5 * features[i].img_size.height;  
    5. }  

          5.使用Bundle Adjustment方法对所有图片进行相机参数校正

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

          光束法平差的目标函数是一个具有鲁棒性的映射误差的平方和函数。即每一个特征点都要映射到其他的图像中,计算出使误差的平方和最小的相机参数。具体的推导过程可以参见Automatic Panoramic Image Stitching using Invariant Features.pdf的第五章,本文只介绍一下OpenCV实现的过程(完整的理论和公式 暂时还没看懂,希望有人能一起交流)

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

    [cpp]  view plain copy
    1. BundleAdjusterReproj // BundleAdjusterBase(7, 2)//最小投影误差平方和  
    2. Implementation of the camera parameters refinement algorithm which minimizes sum of the reprojection error squares  
    3.   
    4. BundleAdjusterRay //  BundleAdjusterBase(4, 3)//最小特征点与相机中心点的距离和  
    5. Implementation of the camera parameters refinement algorithm which minimizes sum of the distances between the  
    6. rays passing through the camera center and a feature.  

          5.1 首先计算cam_params_的值:

    [cpp]  view plain copy
    1. setUpInitialCameraParams(cameras);  

          函数主要源码:

    [cpp]  view plain copy
    1. cam_params_.create(num_images_ * 4, 1, CV_64F);  
    2. SVD svd;//奇异值分解  
    3. for (int i = 0; i < num_images_; ++i)  
    4. {  
    5.     cam_params_.at<double>(i * 4, 0) = cameras[i].focal;  
    6.   
    7.     svd(cameras[i].R, SVD::FULL_UV);  
    8.     Mat R = svd.u * svd.vt;  
    9.     if (determinant(R) < 0)  
    10.         R *= -1;  
    11.   
    12.     Mat rvec;  
    13.     Rodrigues(R, rvec);  
    14.     CV_Assert(rvec.type() == CV_32F);  
    15.     cam_params_.at<double>(i * 4 + 1, 0) = rvec.at<float>(0, 0);  
    16.     cam_params_.at<double>(i * 4 + 2, 0) = rvec.at<float>(1, 0);  
    17.     cam_params_.at<double>(i * 4 + 3, 0) = rvec.at<float>(2, 0);  
    18. }  
          计算cam_params_的值,先初始化cam_params(num_images_*4,1,CV_64F);

          cam_params_[i*4+0] =  cameras[i].focal;

          cam_params_后面3个值,是cameras[i].R先经过奇异值分解,然后对u*vt进行Rodrigues运算,得到的rvec第一行3个值赋给cam_params_。

         奇异值分解的定义:

    在矩阵M的奇异值分解中 M = UΣV*
    ·U的列(columns)组成一套对M的正交"输入"或"分析"的基向量。这些向量是MM*的特征向量。
    ·V的列(columns)组成一套对M的正交"输出"的基向量。这些向量是M*M的特征向量。
    ·Σ对角线上的元素是奇异值,可视为是在输入与输出间进行的标量的"膨胀控制"。这些是M*M及MM*的奇异值,并与U和V的行向量相对应。

         5.2 删除置信度小于门限的匹配对

    [cpp]  view plain copy
    1. // Leave only consistent image pairs  
    2. edges_.clear();  
    3. for (int i = 0; i < num_images_ - 1; ++i)  
    4. {  
    5.     for (int j = i + 1; j < num_images_; ++j)  
    6.     {  
    7.         const MatchesInfo& matches_info = pairwise_matches_[i * num_images_ + j];  
    8.         if (matches_info.confidence > conf_thresh_)  
    9.             edges_.push_back(make_pair(i, j));  
    10.     }  
    11. }  
           5.3 使用LM算法计算camera参数。

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

    [cpp]  view plain copy
    1. //计算所有内点之和  
    2.     for (size_t i = 0; i < edges_.size(); ++i)  
    3.         total_num_matches_ += static_cast<int>(pairwise_matches[edges_[i].first * num_images_ +  
    4.                                                                 edges_[i].second].num_inliers);  
    5.   
    6.     CvLevMarq solver(num_images_ * num_params_per_cam_,  
    7.                      total_num_matches_ * num_errs_per_measurement_,  
    8.                      term_criteria_);  
    9.   
    10.     Mat err, jac;  
    11.     CvMat matParams = cam_params_;  
    12.     cvCopy(&matParams, solver.param);  
    13.   
    14.     int iter = 0;  
    15.     for(;;)//类似于while(1),但是比while(1)效率高  
    16.     {  
    17.         const CvMat* _param = 0;  
    18.         CvMat* _jac = 0;  
    19.         CvMat* _err = 0;  
    20.   
    21.         bool proceed = solver.update(_param, _jac, _err);  
    22.   
    23.         cvCopy(_param, &matParams);  
    24.   
    25.         if (!proceed || !_err)  
    26.             break;  
    27.   
    28.         if (_jac)  
    29.         {  
    30.             calcJacobian(jac); //构造雅阁比行列式  
    31.             CvMat tmp = jac;  
    32.             cvCopy(&tmp, _jac);  
    33.         }  
    34.   
    35.         if (_err)  
    36.         {  
    37.             calcError(err);//计算err  
    38.             LOG_CHAT(".");  
    39.             iter++;  
    40.             CvMat tmp = err;  
    41.             cvCopy(&tmp, _err);  
    42.         }  
    43.     }  
           计算camera

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

    [cpp]  view plain copy
    1. void BundleAdjusterRay::obtainRefinedCameraParams(vector<CameraParams> &cameras) const  
    2. {  
    3.     for (int i = 0; i < num_images_; ++i)  
    4.     {  
    5.         cameras[i].focal = cam_params_.at<double>(i * 4, 0);  
    6.   
    7.         Mat rvec(3, 1, CV_64F);  
    8.         rvec.at<double>(0, 0) = cam_params_.at<double>(i * 4 + 1, 0);  
    9.         rvec.at<double>(1, 0) = cam_params_.at<double>(i * 4 + 2, 0);  
    10.         rvec.at<double>(2, 0) = cam_params_.at<double>(i * 4 + 3, 0);  
    11.         Rodrigues(rvec, cameras[i].R);  
    12.   
    13.         Mat tmp;  
    14.         cameras[i].R.convertTo(tmp, CV_32F);  
    15.         cameras[i].R = tmp;  
    16.     }  
    17. }  
           求出根节点,然后归一化旋转矩阵R
    [cpp]  view plain copy
    1. // Normalize motion to center image  
    2. Graph span_tree;  
    3. vector<int> span_tree_centers;  
    4. findMaxSpanningTree(num_images_, pairwise_matches, span_tree, span_tree_centers);  
    5. Mat R_inv = cameras[span_tree_centers[0]].R.inv();  
    6. for (int i = 0; i < num_images_; ++i)  
    7.     cameras[i].R = R_inv * cameras[i].R;  
         6 波形校正

         前面几节把相机旋转矩阵计算出来,但是还有一个因素需要考虑,就是由于拍摄者拍摄图片的时候不一定是水平的,轻微的倾斜会导致全景图像出现飞机曲线,因此我们要对图像进行波形校正,主要是寻找每幅图形的“上升向量”(up_vector),使他校正成

    波形校正的效果图

             opencv实现的源码如下(也是暂时没看懂,囧)

    [cpp]  view plain copy
    1. <span style="white-space:pre">    </span>vector<Mat> rmats;  
    2.        for (size_t i = 0; i < cameras.size(); ++i)  
    3.            rmats.push_back(cameras[i].R);  
    4.        waveCorrect(rmats, wave_correct);//Tries to make panorama more horizontal (or vertical).  
    5.        for (size_t i = 0; i < cameras.size(); ++i)  
    6.            cameras[i].R = rmats[i];  
           其中waveCorrect(rmats, wave_correct)源码如下:

    [cpp]  view plain copy
    1. void waveCorrect(vector<Mat> &rmats, WaveCorrectKind kind)  
    2. {  
    3.     LOGLN("Wave correcting...");  
    4. #if ENABLE_LOG  
    5.     int64 t = getTickCount();  
    6. #endif  
    7.   
    8.     Mat moment = Mat::zeros(3, 3, CV_32F);  
    9.     for (size_t i = 0; i < rmats.size(); ++i)  
    10.     {  
    11.         Mat col = rmats[i].col(0);  
    12.         moment += col * col.t();//相机R矩阵第一列转置相乘然后相加  
    13.     }  
    14.     Mat eigen_vals, eigen_vecs;  
    15.     eigen(moment, eigen_vals, eigen_vecs);//Calculates eigenvalues and eigenvectors of a symmetric matrix.  
    16.   
    17.     Mat rg1;  
    18.     if (kind == WAVE_CORRECT_HORIZ)  
    19.         rg1 = eigen_vecs.row(2).t();//如果是水平校正,去特征向量的第三行  
    20.     else if (kind == WAVE_CORRECT_VERT)  
    21.         rg1 = eigen_vecs.row(0).t();//如果是垂直校正,特征向量第一行  
    22.     else  
    23.         CV_Error(CV_StsBadArg, "unsupported kind of wave correction");  
    24.   
    25.     Mat img_k = Mat::zeros(3, 1, CV_32F);  
    26.     for (size_t i = 0; i < rmats.size(); ++i)  
    27.         img_k += rmats[i].col(2);//R函数第3列相加  
    28.     Mat rg0 = rg1.cross(img_k);//rg1与img_k向量积  
    29.     rg0 /= norm(rg0);//归一化?  
    30.   
    31.     Mat rg2 = rg0.cross(rg1);  
    32.   
    33.     double conf = 0;  
    34.     if (kind == WAVE_CORRECT_HORIZ)  
    35.     {  
    36.         for (size_t i = 0; i < rmats.size(); ++i)  
    37.             conf += rg0.dot(rmats[i].col(0));//Computes a dot-product of two vectors.数量积  
    38.         if (conf < 0)  
    39.         {  
    40.             rg0 *= -1;  
    41.             rg1 *= -1;  
    42.         }  
    43.     }  
    44.     else if (kind == WAVE_CORRECT_VERT)  
    45.     {  
    46.         for (size_t i = 0; i < rmats.size(); ++i)  
    47.             conf -= rg1.dot(rmats[i].col(0));  
    48.         if (conf < 0)  
    49.         {  
    50.             rg0 *= -1;  
    51.             rg1 *= -1;  
    52.         }  
    53.     }  
    54.   
    55.     Mat R = Mat::zeros(3, 3, CV_32F);  
    56.     Mat tmp = R.row(0);  
    57.     Mat(rg0.t()).copyTo(tmp);  
    58.     tmp = R.row(1);  
    59.     Mat(rg1.t()).copyTo(tmp);  
    60.     tmp = R.row(2);  
    61.     Mat(rg2.t()).copyTo(tmp);  
    62.   
    63.     for (size_t i = 0; i < rmats.size(); ++i)  
    64.         rmats[i] = R * rmats[i];  
    65.   
    66.     LOGLN("Wave correcting, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");  
    67. }  

         7.单应性矩阵变换

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

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

    [cpp]  view plain copy
    1. void SphericalProjector::mapForward(float x, float y, float &u, float &v)  
    2. {  
    3.     float x_ = r_kinv[0] * x + r_kinv[1] * y + r_kinv[2];  
    4.     float y_ = r_kinv[3] * x + r_kinv[4] * y + r_kinv[5];  
    5.     float z_ = r_kinv[6] * x + r_kinv[7] * y + r_kinv[8];  
    6.   
    7.     u = scale * atan2f(x_, z_);  
    8.     float w = y_ / sqrtf(x_ * x_ + y_ * y_ + z_ * z_);  
    9.     v = scale * (static_cast<float>(CV_PI) - acosf(w == w ? w : 0));  
    10. }  

        

           根据映射公式,对图像的上下左右四个边界求映射后的坐标,然后确定变换后图像的左上角和右上角的坐标,

           如果是球形拼接,则需要再加上判断(暂时还没研究透):

    [cpp]  view plain copy
    1. float tl_uf = static_cast<float>(dst_tl.x);  
    2. float tl_vf = static_cast<float>(dst_tl.y);  
    3. float br_uf = static_cast<float>(dst_br.x);  
    4. float br_vf = static_cast<float>(dst_br.y);  
    5.   
    6. float x = projector_.rinv[1];  
    7. float y = projector_.rinv[4];  
    8. float z = projector_.rinv[7];  
    9. if (y > 0.f)  
    10. {  
    11.     float x_ = (projector_.k[0] * x + projector_.k[1] * y) / z + projector_.k[2];  
    12.     float y_ = projector_.k[4] * y / z + projector_.k[5];  
    13.     if (x_ > 0.f && x_ < src_size.width && y_ > 0.f && y_ < src_size.height)  
    14.     {  
    15.         tl_uf = min(tl_uf, 0.f); tl_vf = min(tl_vf, static_cast<float>(CV_PI * projector_.scale));  
    16.         br_uf = max(br_uf, 0.f); br_vf = max(br_vf, static_cast<float>(CV_PI * projector_.scale));  
    17.     }  
    18. }  
    19.   
    20. x = projector_.rinv[1];  
    21. y = -projector_.rinv[4];  
    22. z = projector_.rinv[7];  
    23. if (y > 0.f)  
    24. {  
    25.     float x_ = (projector_.k[0] * x + projector_.k[1] * y) / z + projector_.k[2];  
    26.     float y_ = projector_.k[4] * y / z + projector_.k[5];  
    27.     if (x_ > 0.f && x_ < src_size.width && y_ > 0.f && y_ < src_size.height)  
    28.     {  
    29.         tl_uf = min(tl_uf, 0.f); tl_vf = min(tl_vf, static_cast<float>(0));  
    30.         br_uf = max(br_uf, 0.f); br_vf = max(br_vf, static_cast<float>(0));  
    31.     }  
    32. }  
          然后利用反投影将图形反投影到变换的图像上,像素计算默认是二维线性插值。

         反投影的公式:

    [cpp]  view plain copy
    1. void SphericalProjector::mapBackward(float u, float v, float &x, float &y)  
    2. {  
    3.     u /= scale;  
    4.     v /= scale;  
    5.   
    6.     float sinv = sinf(static_cast<float>(CV_PI) - v);  
    7.     float x_ = sinv * sinf(u);  
    8.     float y_ = cosf(static_cast<float>(CV_PI) - v);  
    9.     float z_ = sinv * cosf(u);  
    10.   
    11.     float z;  
    12.     x = k_rinv[0] * x_ + k_rinv[1] * y_ + k_rinv[2] * z_;  
    13.     y = k_rinv[3] * x_ + k_rinv[4] * y_ + k_rinv[5] * z_;  
    14.     z = k_rinv[6] * x_ + k_rinv[7] * y_ + k_rinv[8] * z_;  
    15.   
    16.     if (z > 0) { x /= z; y /= z; }  
    17.     else x = y = -1;  
    18. }  
    然后将反投影求出的x,y值写入Mat矩阵xmap和ymap中

    [cpp]  view plain copy
    1. xmap.create(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, CV_32F);  
    2.    ymap.create(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, CV_32F);  
    3.   
    4.    float x, y;  
    5.    for (int v = dst_tl.y; v <= dst_br.y; ++v)  
    6.    {  
    7.        for (int u = dst_tl.x; u <= dst_br.x; ++u)  
    8.        {  
    9.            projector_.mapBackward(static_cast<float>(u), static_cast<float>(v), x, y);  
    10.            xmap.at<float>(v - dst_tl.y, u - dst_tl.x) = x;  
    11.            ymap.at<float>(v - dst_tl.y, u - dst_tl.x) = y;  
    12.        }  
    13.    }  
    最后使用opencv自带的remap函数将图像重新绘制:

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

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

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

    [cpp]  view plain copy
    1. for (int img_idx = 0; img_idx < num_images; ++img_idx)  
    2.   {  
    3.       Size bl_per_img((images[img_idx].cols + bl_width_ - 1) / bl_width_,  
    4.                       (images[img_idx].rows + bl_height_ - 1) / bl_height_);  
    5.       int bl_width = (images[img_idx].cols + bl_per_img.width - 1) / bl_per_img.width;  
    6.       int bl_height = (images[img_idx].rows + bl_per_img.height - 1) / bl_per_img.height;  
    7.       bl_per_imgs[img_idx] = bl_per_img;  
    8.       for (int by = 0; by < bl_per_img.height; ++by)  
    9.       {  
    10.           for (int bx = 0; bx < bl_per_img.width; ++bx)  
    11.           {  
    12.               Point bl_tl(bx * bl_width, by * bl_height);  
    13.               Point bl_br(min(bl_tl.x + bl_width, images[img_idx].cols),  
    14.                           min(bl_tl.y + bl_height, images[img_idx].rows));  
    15.   
    16.               block_corners.push_back(corners[img_idx] + bl_tl);  
    17.               block_images.push_back(images[img_idx](Rect(bl_tl, bl_br)));  
    18.               block_masks.push_back(make_pair(masks[img_idx].first(Rect(bl_tl, bl_br)),  
    19.                                               masks[img_idx].second));  
    20.           }  
    21.       }  
    22.   }  

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

    [cpp]  view plain copy
    1. //计算每一块区域的光照均值sqrt(sqrt(R)+sqrt(G)+sqrt(B))  
    2.     //光照均值是对称矩阵,所以一次循环计算两个光照值,(i,j),与(j,i)  
    3.     for (int i = 0; i < num_images; ++i)  
    4.     {  
    5.         for (int j = i; j < num_images; ++j)  
    6.         {  
    7.             Rect roi;  
    8.             //判断image[i]与image[j]是否有重叠部分  
    9.             if (overlapRoi(corners[i], corners[j], images[i].size(), images[j].size(), roi))  
    10.             {  
    11.                 subimg1 = images[i](Rect(roi.tl() - corners[i], roi.br() - corners[i]));  
    12.                 subimg2 = images[j](Rect(roi.tl() - corners[j], roi.br() - corners[j]));  
    13.   
    14.                 submask1 = masks[i].first(Rect(roi.tl() - corners[i], roi.br() - corners[i]));  
    15.                 submask2 = masks[j].first(Rect(roi.tl() - corners[j], roi.br() - corners[j]));  
    16.                 intersect = (submask1 == masks[i].second) & (submask2 == masks[j].second);  
    17.   
    18.                 N(i, j) = N(j, i) = max(1, countNonZero(intersect));  
    19.   
    20.                 double Isum1 = 0, Isum2 = 0;  
    21.                 for (int y = 0; y < roi.height; ++y)  
    22.                 {  
    23.                     const Point3_<uchar>* r1 = subimg1.ptr<Point3_<uchar> >(y);  
    24.                     const Point3_<uchar>* r2 = subimg2.ptr<Point3_<uchar> >(y);  
    25.                     for (int x = 0; x < roi.width; ++x)  
    26.                     {  
    27.                         if (intersect(y, x))  
    28.                         {  
    29.                             Isum1 += sqrt(static_cast<double>(sqr(r1[x].x) + sqr(r1[x].y) + sqr(r1[x].z)));  
    30.                             Isum2 += sqrt(static_cast<double>(sqr(r2[x].x) + sqr(r2[x].y) + sqr(r2[x].z)));  
    31.                         }  
    32.                     }  
    33.                 }  
    34.                 I(i, j) = Isum1 / N(i, j);  
    35.                 I(j, i) = Isum2 / N(i, j);  
    36.             }  
    37.         }  
    38.     }  
         建立方程,求出每个光强的调整系数

    [cpp]  view plain copy
    1. Mat_<double> A(num_images, num_images); A.setTo(0);  
    2. Mat_<double> b(num_images, 1); b.setTo(0);//beta*N(i,j)  
    3. for (int i = 0; i < num_images; ++i)  
    4. {  
    5.     for (int j = 0; j < num_images; ++j)  
    6.     {  
    7.         b(i, 0) += beta * N(i, j);  
    8.         A(i, i) += beta * N(i, j);  
    9.         if (j == i) continue;  
    10.         A(i, i) += 2 * alpha * I(i, j) * I(i, j) * N(i, j);  
    11.         A(i, j) -= 2 * alpha * I(i, j) * I(j, i) * N(i, j);  
    12.     }  
    13. }  
    14.   
    15. solve(A, b, gains_);//求方程的解A*gains = B  

            gains_原理分析:

    num_images :表示图像分块的个数,零num_images = n

    N矩阵,大小n*n,N(i,j)表示第i幅图像与第j幅图像重合的像素点数,N(i,j)=N(j,i)

    I矩阵,大小n*n,I(i,j)与I(j,i)表示第i,j块区域重合部分的像素平均值,I(i,j)是重合区域中i快的平均亮度值,


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

    A矩阵,大小n*n,公式图片不全


    b矩阵,大小n*1,


    然后根据求解矩阵

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

            然后将gains_进行线性滤波

    [cpp]  view plain copy
    1.   Mat_<float> ker(1, 3);  
    2.   ker(0,0) = 0.25; ker(0,1) = 0.5; ker(0,2) = 0.25;  
    3.   
    4.   int bl_idx = 0;  
    5.   for (int img_idx = 0; img_idx < num_images; ++img_idx)  
    6.   {  
    7. Size bl_per_img = bl_per_imgs[img_idx];  
    8. gain_maps_[img_idx].create(bl_per_img);  
    9.   
    10.       for (int by = 0; by < bl_per_img.height; ++by)  
    11.           for (int bx = 0; bx < bl_per_img.width; ++bx, ++bl_idx)  
    12.               gain_maps_[img_idx](by, bx) = static_cast<float>(gains[bl_idx]);  
    13. //用分解的核函数对图像做卷积。首先,图像的每一行与一维的核kernelX做卷积;然后,运算结果的每一列与一维的核kernelY做卷积  
    14.       sepFilter2D(gain_maps_[img_idx], gain_maps_[img_idx], CV_32F, ker, ker);  
    15.       sepFilter2D(gain_maps_[img_idx], gain_maps_[img_idx], CV_32F, ker, ker);  
    16.   }  

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

       

       9.Seam Estimation

         缝隙估计有6种方法,默认就是第三种方法,seam_find_type == "gc_color",该方法是利用最大流方法检测。

    [cpp]  view plain copy
    1.  if (seam_find_type == "no")  
    2.         seam_finder = new detail::NoSeamFinder();//Stub seam estimator which does nothing.  
    3.     else if (seam_find_type == "voronoi")  
    4.         seam_finder = new detail::VoronoiSeamFinder();//Voronoi diagram-based seam estimator. 泰森多边形缝隙估计  
    5.     else if (seam_find_type == "gc_color")  
    6.     {  
    7. #ifdef HAVE_OPENCV_GPU  
    8.         if (try_gpu && gpu::getCudaEnabledDeviceCount() > 0)  
    9.             seam_finder = new detail::GraphCutSeamFinderGpu(GraphCutSeamFinderBase::COST_COLOR);  
    10.         else  
    11. #endif  
    12.             seam_finder = new detail::GraphCutSeamFinder(GraphCutSeamFinderBase::COST_COLOR);//Minimum graph cut-based seam estimator  
    13.     }  
    14.     else if (seam_find_type == "gc_colorgrad")  
    15.     {  
    16. #ifdef HAVE_OPENCV_GPU  
    17.         if (try_gpu && gpu::getCudaEnabledDeviceCount() > 0)  
    18.             seam_finder = new detail::GraphCutSeamFinderGpu(GraphCutSeamFinderBase::COST_COLOR_GRAD);  
    19.         else  
    20. #endif  
    21.             seam_finder = new detail::GraphCutSeamFinder(GraphCutSeamFinderBase::COST_COLOR_GRAD);  
    22.     }  
    23.     else if (seam_find_type == "dp_color")  
    24.         seam_finder = new detail::DpSeamFinder(DpSeamFinder::COLOR);  
    25.     else if (seam_find_type == "dp_colorgrad")  
    26.         seam_finder = new detail::DpSeamFinder(DpSeamFinder::COLOR_GRAD);  
    27.     if (seam_finder.empty())  
    28.     {  
    29.         cout << "Can't create the following seam finder '" << seam_find_type << "'\n";  
    30.         return 1;  
    31.     }  
          程序解读可参见:

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

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

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

         

         10.多波段融合

          由于以前进行处理的图片都是以work_scale(3.1节有介绍)进行缩放的,所以图像的内参,corner(同一坐标后的顶点),mask(融合的掩码)都需要重新计算。以及将之前计算的光照增强的gain也要计算进去。

    [cpp]  view plain copy
    1. // Read image and resize it if necessary  
    2.       full_img = imread(img_names[img_idx]);  
    3.       if (!is_compose_scale_set)  
    4.       {  
    5.           if (compose_megapix > 0)  
    6.               compose_scale = min(1.0, sqrt(compose_megapix * 1e6 / full_img.size().area()));  
    7.           is_compose_scale_set = true;  
    8.   
    9.           // Compute relative scales  
    10.           //compose_seam_aspect = compose_scale / seam_scale;  
    11.           compose_work_aspect = compose_scale / work_scale;  
    12.   
    13.           // Update warped image scale  
    14.           warped_image_scale *= static_cast<float>(compose_work_aspect);  
    15.           warper = warper_creator->create(warped_image_scale);  
    16.   
    17.           // Update corners and sizes  
    18.           for (int i = 0; i < num_images; ++i)  
    19.           {  
    20.               // Update intrinsics  
    21.               cameras[i].focal *= compose_work_aspect;  
    22.               cameras[i].ppx *= compose_work_aspect;  
    23.               cameras[i].ppy *= compose_work_aspect;  
    24.   
    25.               // Update corner and size  
    26.               Size sz = full_img_sizes[i];  
    27.               if (std::abs(compose_scale - 1) > 1e-1)  
    28.               {  
    29.                   sz.width = cvRound(full_img_sizes[i].width * compose_scale);//取整  
    30.                   sz.height = cvRound(full_img_sizes[i].height * compose_scale);  
    31.               }  
    32.   
    33.               Mat K;  
    34.               cameras[i].K().convertTo(K, CV_32F);  
    35.               Rect roi = warper->warpRoi(sz, K, cameras[i].R);//Returns Projected image minimum bounding box  
    36.               corners[i] = roi.tl();//! the top-left corner  
    37.               sizes[i] = roi.size();//! size of the real buffer  
    38.           }  
    39.       }  
    40.       if (abs(compose_scale - 1) > 1e-1)  
    41.           resize(full_img, img, Size(), compose_scale, compose_scale);  
    42.       else  
    43.           img = full_img;  
    44.       full_img.release();  
    45.       Size img_size = img.size();  
    46.   
    47.       Mat K;  
    48.       cameras[img_idx].K().convertTo(K, CV_32F);  
    49.   
    50.       // Warp the current image  
    51.       warper->warp(img, K, cameras[img_idx].R, INTER_LINEAR, BORDER_REFLECT, img_warped);  
    52.       // Warp the current image mask  
    53.       mask.create(img_size, CV_8U);  
    54.       mask.setTo(Scalar::all(255));  
    55.       warper->warp(mask, K, cameras[img_idx].R, INTER_NEAREST, BORDER_CONSTANT, mask_warped);  
    56.       // Compensate exposure  
    57.       compensator->apply(img_idx, corners[img_idx], img_warped, mask_warped);//光照补偿  
    58.       img_warped.convertTo(img_warped_s, CV_16S);  
    59.       img_warped.release();  
    60.       img.release();  
    61.       mask.release();  
    62.   
    63.       dilate(masks_warped[img_idx], dilated_mask, Mat());  
    64.       resize(dilated_mask, seam_mask, mask_warped.size());  
    65.       mask_warped = seam_mask & mask_warped;  
         对图像进行光照补偿,就是将对应区域乘以gain:

    [cpp]  view plain copy
    1. //块光照补偿  
    2. void BlocksGainCompensator::apply(int index, Point /*corner*/, Mat &image, const Mat &/*mask*/)  
    3. {  
    4.     CV_Assert(image.type() == CV_8UC3);  
    5.   
    6.     Mat_<float> gain_map;  
    7.     if (gain_maps_[index].size() == image.size())  
    8.         gain_map = gain_maps_[index];  
    9.     else  
    10.         resize(gain_maps_[index], gain_map, image.size(), 0, 0, INTER_LINEAR);  
    11.   
    12.     for (int y = 0; y < image.rows; ++y)  
    13.     {  
    14.         const float* gain_row = gain_map.ptr<float>(y);  
    15.         Point3_<uchar>* row = image.ptr<Point3_<uchar> >(y);  
    16.         for (int x = 0; x < image.cols; ++x)  
    17.         {  
    18.             row[x].x = saturate_cast<uchar>(row[x].x * gain_row[x]);  
    19.             row[x].y = saturate_cast<uchar>(row[x].y * gain_row[x]);  
    20.             row[x].z = saturate_cast<uchar>(row[x].z * gain_row[x]);  
    21.         }  
    22.     }  
    23. }  

         进行多波段融合,首先初始化blend,确定blender的融合的方式,默认是多波段融合MULTI_BAND,以及根据corners顶点和图像的大小确定最终全景图的尺寸。

    [cpp]  view plain copy
    1. <span>    </span>//初始化blender  
    2.         if (blender.empty())  
    3.         {  
    4.             blender = Blender::createDefault(blend_type, try_gpu);  
    5.             Size dst_sz = resultRoi(corners, sizes).size();//计算最后图像的大小  
    6.             float blend_width = sqrt(static_cast<float>(dst_sz.area())) * blend_strength / 100.f;  
    7.             if (blend_width < 1.f)  
    8.                 blender = Blender::createDefault(Blender::NO, try_gpu);  
    9.             else if (blend_type == Blender::MULTI_BAND)  
    10.             {  
    11.                 MultiBandBlender* mb = dynamic_cast<MultiBandBlender*>(static_cast<Blender*>(blender));  
    12.                 mb->setNumBands(static_cast<int>(ceil(log(blend_width)/log(2.)) - 1.));  
    13.                 LOGLN("Multi-band blender, number of bands: " << mb->numBands());  
    14.             }  
    15.             else if (blend_type == Blender::FEATHER)  
    16.             {  
    17.                 FeatherBlender* fb = dynamic_cast<FeatherBlender*>(static_cast<Blender*>(blender));  
    18.                 fb->setSharpness(1.f/blend_width);  
    19.                 LOGLN("Feather blender, sharpness: " << fb->sharpness());  
    20.             }  
    21.             blender->prepare(corners, sizes);//根据corners顶点和图像的大小确定最终全景图的尺寸  
    22.         }  
          然后对每幅图图形构建金字塔,层数可以由输入的系数确定,默认是5层。

          先对顶点以及图像的宽和高做处理,使其能被2^num_bands除尽,这样才能将进行向下采样num_bands次,首先从源图像pyrDown向下采样,在由最底部的图像pyrUp向上采样,把对应层得上采样和下采样的相减,就得到了图像的高频分量,存储到每一个金字塔中。然后根据mask,将每幅图像的各层金字塔分别写入最终的金字塔层src_pyr_laplace中。

          最后将各层的金字塔叠加,就得到了最终完整的全景图。

    [cpp]  view plain copy
    1. blender->feed(img_warped_s, mask_warped, corners[img_idx]);//将图像写入金字塔中  
          源码:

    [cpp]  view plain copy
    1. void MultiBandBlender::feed(const Mat &img, const Mat &mask, Point tl)  
    2. {  
    3.     CV_Assert(img.type() == CV_16SC3 || img.type() == CV_8UC3);  
    4.     CV_Assert(mask.type() == CV_8U);  
    5.   
    6.     // Keep source image in memory with small border  
    7.     int gap = 3 * (1 << num_bands_);  
    8.     Point tl_new(max(dst_roi_.x, tl.x - gap),  
    9.                  max(dst_roi_.y, tl.y - gap));  
    10.     Point br_new(min(dst_roi_.br().x, tl.x + img.cols + gap),  
    11.                  min(dst_roi_.br().y, tl.y + img.rows + gap));  
    12.   
    13.     // Ensure coordinates of top-left, bottom-right corners are divided by (1 << num_bands_).  
    14.     // After that scale between layers is exactly 2.  
    15.     //  
    16.     // We do it to avoid interpolation problems when keeping sub-images only. There is no such problem when  
    17.     // image is bordered to have size equal to the final image size, but this is too memory hungry approach.  
    18.     //将顶点调整成能被2^num_bank次方除尽的值  
    19.     tl_new.x = dst_roi_.x + (((tl_new.x - dst_roi_.x) >> num_bands_) << num_bands_);  
    20.     tl_new.y = dst_roi_.y + (((tl_new.y - dst_roi_.y) >> num_bands_) << num_bands_);  
    21.     int width = br_new.x - tl_new.x;  
    22.     int height = br_new.y - tl_new.y;  
    23.     width += ((1 << num_bands_) - width % (1 << num_bands_)) % (1 << num_bands_);  
    24.     height += ((1 << num_bands_) - height % (1 << num_bands_)) % (1 << num_bands_);  
    25.     br_new.x = tl_new.x + width;  
    26.     br_new.y = tl_new.y + height;  
    27.     int dy = max(br_new.y - dst_roi_.br().y, 0);  
    28.     int dx = max(br_new.x - dst_roi_.br().x, 0);  
    29.     tl_new.x -= dx; br_new.x -= dx;  
    30.     tl_new.y -= dy; br_new.y -= dy;  
    31.   
    32.     int top = tl.y - tl_new.y;  
    33.     int left = tl.x - tl_new.x;  
    34.     int bottom = br_new.y - tl.y - img.rows;  
    35.     int right = br_new.x - tl.x - img.cols;  
    36.   
    37.     // Create the source image Laplacian pyramid  
    38.     Mat img_with_border;  
    39.     copyMakeBorder(img, img_with_border, top, bottom, left, right,  
    40.                    BORDER_REFLECT);//给图像设置一个边界,BORDER_REFLECT边界颜色任意  
    41.     vector<Mat> src_pyr_laplace;  
    42.     if (can_use_gpu_ && img_with_border.depth() == CV_16S)  
    43.         createLaplacePyrGpu(img_with_border, num_bands_, src_pyr_laplace);  
    44.     else  
    45.         createLaplacePyr(img_with_border, num_bands_, src_pyr_laplace);//创建高斯金字塔,每一层保存的全是高频信息  
    46.   
    47.     // Create the weight map Gaussian pyramid  
    48.     Mat weight_map;  
    49.     vector<Mat> weight_pyr_gauss(num_bands_ + 1);  
    50.   
    51.     if(weight_type_ == CV_32F)  
    52.     {  
    53.         mask.convertTo(weight_map, CV_32F, 1./255.);//将mask的0,255归一化成0,1  
    54.     }  
    55.     else// weight_type_ == CV_16S  
    56.     {  
    57.         mask.convertTo(weight_map, CV_16S);  
    58.         add(weight_map, 1, weight_map, mask != 0);  
    59.     }  
    60.   
    61.     copyMakeBorder(weight_map, weight_pyr_gauss[0], top, bottom, left, right, BORDER_CONSTANT);  
    62.   
    63.     for (int i = 0; i < num_bands_; ++i)  
    64.         pyrDown(weight_pyr_gauss[i], weight_pyr_gauss[i + 1]);  
    65.   
    66.     int y_tl = tl_new.y - dst_roi_.y;  
    67.     int y_br = br_new.y - dst_roi_.y;  
    68.     int x_tl = tl_new.x - dst_roi_.x;  
    69.     int x_br = br_new.x - dst_roi_.x;  
    70.   
    71.     // Add weighted layer of the source image to the final Laplacian pyramid layer  
    72.     if(weight_type_ == CV_32F)  
    73.     {  
    74.         for (int i = 0; i <= num_bands_; ++i)  
    75.         {  
    76.             for (int y = y_tl; y < y_br; ++y)  
    77.             {  
    78.                 int y_ = y - y_tl;  
    79.                 const Point3_<short>* src_row = src_pyr_laplace[i].ptr<Point3_<short> >(y_);  
    80.                 Point3_<short>* dst_row = dst_pyr_laplace_[i].ptr<Point3_<short> >(y);  
    81.                 const float* weight_row = weight_pyr_gauss[i].ptr<float>(y_);  
    82.                 float* dst_weight_row = dst_band_weights_[i].ptr<float>(y);  
    83.   
    84.                 for (int x = x_tl; x < x_br; ++x)  
    85.                 {  
    86.                     int x_ = x - x_tl;  
    87.                     dst_row[x].x += static_cast<short>(src_row[x_].x * weight_row[x_]);  
    88.                     dst_row[x].y += static_cast<short>(src_row[x_].y * weight_row[x_]);  
    89.                     dst_row[x].z += static_cast<short>(src_row[x_].z * weight_row[x_]);  
    90.                     dst_weight_row[x] += weight_row[x_];  
    91.                 }  
    92.             }  
    93.             x_tl /= 2; y_tl /= 2;  
    94.             x_br /= 2; y_br /= 2;  
    95.         }  
    96.     }  
    97.     else// weight_type_ == CV_16S  
    98.     {  
    99.         for (int i = 0; i <= num_bands_; ++i)  
    100.         {  
    101.             for (int y = y_tl; y < y_br; ++y)  
    102.             {  
    103.                 int y_ = y - y_tl;  
    104.                 const Point3_<short>* src_row = src_pyr_laplace[i].ptr<Point3_<short> >(y_);  
    105.                 Point3_<short>* dst_row = dst_pyr_laplace_[i].ptr<Point3_<short> >(y);  
    106.                 const short* weight_row = weight_pyr_gauss[i].ptr<short>(y_);  
    107.                 short* dst_weight_row = dst_band_weights_[i].ptr<short>(y);  
    108.   
    109.                 for (int x = x_tl; x < x_br; ++x)  
    110.                 {  
    111.                     int x_ = x - x_tl;  
    112.                     dst_row[x].x += short((src_row[x_].x * weight_row[x_]) >> 8);  
    113.                     dst_row[x].y += short((src_row[x_].y * weight_row[x_]) >> 8);  
    114.                     dst_row[x].z += short((src_row[x_].z * weight_row[x_]) >> 8);  
    115.                     dst_weight_row[x] += weight_row[x_];  
    116.                 }  
    117.             }  
    118.             x_tl /= 2; y_tl /= 2;  
    119.             x_br /= 2; y_br /= 2;  
    120.         }  
    121.     }  
    122. }  
            其中,金字塔构建的源码:
    [cpp]  view plain copy
    1. void createLaplacePyr(const Mat &img, int num_levels, vector<Mat> &pyr)  
    2. {  
    3. #ifdef HAVE_TEGRA_OPTIMIZATION  
    4.     if(tegra::createLaplacePyr(img, num_levels, pyr))  
    5.         return;  
    6. #endif  
    7.   
    8.     pyr.resize(num_levels + 1);  
    9.   
    10.     if(img.depth() == CV_8U)  
    11.     {  
    12.         if(num_levels == 0)  
    13.         {  
    14.             img.convertTo(pyr[0], CV_16S);  
    15.             return;  
    16.         }  
    17.   
    18.         Mat downNext;  
    19.         Mat current = img;  
    20.         pyrDown(img, downNext);  
    21.   
    22.         for(int i = 1; i < num_levels; ++i)  
    23.         {  
    24.             Mat lvl_up;  
    25.             Mat lvl_down;  
    26.   
    27.             pyrDown(downNext, lvl_down);  
    28.             pyrUp(downNext, lvl_up, current.size());  
    29.             subtract(current, lvl_up, pyr[i-1], noArray(), CV_16S);  
    30.   
    31.             current = downNext;  
    32.             downNext = lvl_down;  
    33.         }  
    34.   
    35.         {  
    36.             Mat lvl_up;  
    37.             pyrUp(downNext, lvl_up, current.size());  
    38.             subtract(current, lvl_up, pyr[num_levels-1], noArray(), CV_16S);  
    39.   
    40.             downNext.convertTo(pyr[num_levels], CV_16S);  
    41.         }  
    42.     }  
    43.     else  
    44.     {  
    45.         pyr[0] = img;  
    46.         //构建高斯金字塔  
    47.         for (int i = 0; i < num_levels; ++i)  
    48.             pyrDown(pyr[i], pyr[i + 1]);//先高斯滤波,在亚采样,得到比pyr【i】缩小一半的图像  
    49.         Mat tmp;  
    50.         for (int i = 0; i < num_levels; ++i)  
    51.         {  
    52.             pyrUp(pyr[i + 1], tmp, pyr[i].size());//插值(偶数行,偶数列赋值为0),然后高斯滤波,核是5*5。  
    53.             subtract(pyr[i], tmp, pyr[i]);//pyr[i] = pyr[i]-tmp,得到的全是高频信息  
    54.         }  
    55.     }  
    56. }  
          最终把所有层得金字塔叠加的程序:

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

    [cpp]  view plain copy
    1. void MultiBandBlender::blend(Mat &dst, Mat &dst_mask)  
    2. {  
    3.     for (int i = 0; i <= num_bands_; ++i)  
    4.         normalizeUsingWeightMap(dst_band_weights_[i], dst_pyr_laplace_[i]);  
    5.   
    6.     if (can_use_gpu_)  
    7.         restoreImageFromLaplacePyrGpu(dst_pyr_laplace_);  
    8.     else  
    9.         restoreImageFromLaplacePyr(dst_pyr_laplace_);  
    10.   
    11.     dst_ = dst_pyr_laplace_[0];  
    12.     dst_ = dst_(Range(0, dst_roi_final_.height), Range(0, dst_roi_final_.width));  
    13.     dst_mask_ = dst_band_weights_[0] > WEIGHT_EPS;  
    14.     dst_mask_ = dst_mask_(Range(0, dst_roi_final_.height), Range(0, dst_roi_final_.width));  
    15.     dst_pyr_laplace_.clear();  
    16.     dst_band_weights_.clear();  
    17.   
    18.     Blender::blend(dst, dst_mask);  
    19. }  

    展开全文
  • opencv 金子塔采样

    2012-09-20 15:53:17
    基于opencv里面函数实现的金字塔采样,可以运行一下看看采样效果。
  • opencv采样和下采样

    2019-04-21 14:54:59
    #include "stdafx.h" ...opencv2\highgui\highgui.hpp> #include<opencv\cv.h> #include<opencv2\core\core.hpp> #include<opencv2\opencv.hpp> #include "opencv2/imgproc/imgproc....


    #include "stdafx.h"
    #include<opencv2\highgui\highgui.hpp>
    #include<opencv\cv.h>
    #include<opencv2\core\core.hpp>
    #include<opencv2\opencv.hpp>
    #include "opencv2/imgproc/imgproc.hpp"
    #include<iostream>
    #include "imgproc/imgproc.hpp"
    using namespace cv;
    using namespace std;

    void on_MouseHandle(int evev, int x, int y, int flags, void* param);
    void DrawRectangle(Mat& img, Rect box);
    bool g_bDrawingBox = false;
    Rect g_rectangle;
    int main()
    {
        string name = "原图";
        Mat qq = imread("D:\\test11.jpg");
        Mat qq1,qq2;

        imshow(name, qq);
        pyrUp(qq, qq1, Size(qq.cols * 2, qq.rows * 2));//上采样
        imshow("上采样", qq1);
        pyrDown(qq, qq2, Size(qq.cols /2, qq.rows /2));//下采样
        imshow("下采样", qq2);
        waitKey(0);
        return 0;
    }


        

    展开全文
  • 《RANSAC(随机采样一致算法)原理及openCV代码实现》 原文: http://www.lai18.com/content/1046939.html  本文转自:http://blog.csdn.net/yihaizhiyan/article/details/5973729 ...

    《RANSAC(随机采样一致算法)原理及openCV代码实现》

    原文: http://www.lai18.com/content/1046939.html

    
    本文转自:http://blog.csdn.net/yihaizhiyan/article/details/5973729

    http://blog.csdn.net/Sway_2012/article/details/37765765

    http://blog.csdn.net/zouwen198317/article/details/38494149

    1.什么是RANSAC?

    RANSAC是 RAN dom  SA mple  C onsensus(随机抽样一致性)的缩写。它是从一个观察数据集合中,估计模型参数(模型拟合)的迭代方法。它是一种随机的不确定算法,每次运算求出的结果可能不相同,但总能给出一个合理的结果,为了提高概率必须提高迭代次数。

    2.算法详解

    给定两个点p1与p2的坐标,确定这两点所构成的直线,要求对于输入的任意点p3,都可以判断它是否在该直线上。初中解析几何知识告诉我们,判断一个点在直线上,只需其与直线上任意两点点斜率都相同即可。实际操作当中,往往会先根据已知的两点算出直线的表达式(点斜式、截距式等等),然后通过向量计算即可方便地判断p3是否在该直线上。 

    生产实践中的数据往往会有一定的偏差。例如我们知道两个变量X与Y之间呈线性关系,Y=aX+b,我们想确定参数a与b的具体值。通过实验,可以得到一组X与Y的测试值。虽然理论上两个未知数的方程只需要两组值即可确认,但由于系统误差的原因,任意取两点算出的a与b的值都不尽相同。我们希望的是,最后计算得出的理论模型与测试值的误差最小。大学的高等数学课程中,详细阐述了最小二乘法的思想。通过计算最小均方差关于参数a、b的偏导数为零时的值。事实上,在很多情况下,最小二乘法都是线性回归的代名词。 

    遗憾的是,最小二乘法只适合与误差较小的情况。试想一下这种情况,假使需要从一个噪音较大的数据集中提取模型(比方说只有20%的数据时符合模型的)时,最小二乘法就显得力不从心了。例如下图,肉眼可以很轻易地看出一条直线(模式),但算法却找错了。 



    RANSAC算法的输入是一组观测数据(往往含有较大的噪声或无效点),一个用于解释观测数据的参数化模型以及一些可信的参数。RANSAC通过反复选择数据中的一组随机子集来达成目标。被选取的子集被假设为局内点,并用下述方法进行验证: 

    有一个模型适应于假设的局内点,即所有的未知参数都能从假设的局内点计算得出。
    用1中得到的模型去测试所有的其它数据,如果某个点适用于估计的模型,认为它也是局内点。
    如果有足够多的点被归类为假设的局内点,那么估计的模型就足够合理。
    然后,用所有假设的局内点去重新估计模型(譬如使用最小二乘法),因为它仅仅被初始的假设局内点估计过。
    最后,通过估计局内点与模型的错误率来评估模型。
    上述过程被重复执行固定的次数,每次产生的模型要么因为局内点太少而被舍弃,要么因为比现有的模型更好而被选用。

    整个过程可参考下图: 



    3.代码实现

    随机一致性采样RANSAC是一种鲁棒的模型拟合算法,能够从有外点的数据中拟合准确的模型。

    RANSAC过程中用到的参数

    N-- 拟合模型所需要的最少的样本个数

    K--算法的迭代次数

    t--用于判断数据是否是内点

    d--判定模型是否符合使用于数据集,也就是判断是否是好的模型

    RANSAC算法过程

    1 for K 次迭代

    2 从数据中均匀随机采样N个点

    3 利用采样的N个点拟合你个模型

    4 for 对于除采样点外的每一个样本点

    5 利用t检测样本点到模型的距离,如果小于t则认为是一致,否则认为是外点

    6 end

    7 如果有d或者更多的一致点,则认为拟合的模型是好的

    8 end

    9 使用拟合误差作为标准,选择最好的拟合模型

    迭代次数的计算

    假设 r = 内点个数/所有点的个数

    则:

    p0 = pow(r, N) 表示采样的N个点全为内点,也就是是一次有效采样的概率

    p1 = 1 - pow(r, N) 表示采样的N个点中至少有一个外点,即一次无效采样的概率

    p2 = pow(p1, K) 表示K次无效采样的概率

    假设p表示K次采样中至少一次采样是有效采样,则有1-p = pow(p1, K), 两边取对数

    则有 K = log(1- p )/log(1-p1).

    附一份来自google 的RANSAC的代码框架

    [cpp]  
    view plain copy print ?

    #ifndef FVISION_RANSAC_H_ 
    #define FVISION_RANSAC_H_ 

    #include <fvision/utils/random_utils.h> 
    #include <fvision/utils/misc.h> 

    #include <vector> 
    #include <iostream> 
    #include <cassert> 

    namespace fvision { 

    class RANSAC_SamplesNumber { 
    public: 
    RANSAC_SamplesNumber(int modelSampleSize) { 
    this->s = modelSampleSize; 
    this->p = 0.99; 

    ~RANSAC_SamplesNumber(void) {} 

    public: 
    long calcN(int inliersNumber, int samplesNumber) { 
    double e = 1 - (double)inliersNumber / samplesNumber; 
    //cout<<"e: "<<e<<endl; 
    if (e > 0.9) e = 0.9; 
    //cout<<"pow: "<<pow((1 - e), s)<<endl; 
    //cout<<log(1 - pow((1 - e), s))<<endl; 
    long N = (long)(log(1 - p) / log(1 - pow((1 - e), s))); 
    if (N < 0) return (long) 1000000000
    else return N; 


    private: 
    int s; //samples size for fitting a model 
    double p; //probability that at least one of the random samples if free from outliers 
    //usually 0.99 
    }; 

    //fit a model to a set of samples 
    template <typename M, typename S> 
    class GenericModelCalculator { 
    public: 
    typedef std::vector<S> Samples; 
    virtual M compute(const Samples& samples) = 0; 

    virtual ~GenericModelCalculator<M, S>() {} 

    //the model calculator may only use a subset of the samples for computing 
    //default return empty for both 
    virtual const std::vector<int>& getInlierIndices() const { return defaultInlierIndices; }; 
    virtual const std::vector<int>& getOutlierIndices() const { return defaultOutlierIndices; }; 

    // if the subclass has a threshold parameter, it need to override the following three functions 
    // this is used for algorithms which have a normalization step on input samples 
    virtual bool hasThreshold() const { return false; } 
    virtual void setThreshold(double threshold) {} 
    virtual double getThreshold() const { return 0; } 

    protected: 
    std::vector<int> defaultInlierIndices; 
    std::vector<int> defaultOutlierIndices; 
    }; 

    //evaluate a model to samples 
    //using a threshold to distinguish inliers and outliers 
    template <typename M, typename S> 
    class GenericErrorCaclculator { 
    public: 
    virtual ~GenericErrorCaclculator<M, S>() {} 

    typedef std::vector<S> Samples; 

    virtual double compute(const M& model, const S& sample) const = 0; 

    double computeAverage(const M& model, const Samples& samples) const { 
    int n = (int)samples.size(); 
    if (n == 0) return 0; 
    double sum = 0; 
    for (int i = 0; i < n; i++) { 
    sum += compute(model, samples[i]); 

    return sum / n; 


    double computeInlierAverage(const M& model, const Samples& samples) const { 
    int n = (int)samples.size(); 
    if (n == 0) return 0; 
    double sum = 0; 
    double error = 0; 
    int inlierNum = 0; 
    for (int i = 0; i < n; i++) { 
    error = compute(model, samples[i]); 
    if (error <= threshold) { 
    sum += error; 
    inlierNum++; 


    if (inlierNum == 0) return  1000000
    return sum / inlierNum; 


    public: 

    /** set a threshold for classify inliers and outliers 
    */ 
    void setThreshold(double v) { threshold = v; } 

    double getThreshold() const { return threshold; } 

    /** classify all samples to inliers and outliers 

    */ 
    void classify(const M& model, const Samples& samples, Samples& inliers, Samples& outliers) const { 
    inliers.clear(); 
    outliers.clear(); 
    Samples::const_iterator iter = samples.begin(); 
    for (; iter != samples.end(); ++iter) { 
    if (isInlier(model, *iter)) inliers.push_back(*iter); 
    else outliers.push_back(*iter); 



    /** classify all samples to inliers and outliers, output indices 

    */ 
    void classify(const M& model, const Samples& samples, std::vector<int>& inlierIndices, std::vector<int>& outlierIndices) const { 
    inlierIndices.clear(); 
    outlierIndices.clear(); 
    Samples::const_iterator iter = samples.begin(); 
    int i = 0; 
    for (; iter != samples.end(); ++iter, ++i) { 
    if (isInlier(model, *iter)) inlierIndices.push_back(i); 
    else outlierIndices.push_back(i); 



    /** classify all samples to inliers and outliers 

    */ 
    void classify(const M& model, const Samples& samples, 
    std::vector<int>& inlierIndices, std::vector<int>& outlierIndices, 
    Samples& inliers, Samples& outliers) const { 

    inliers.clear(); 
    outliers.clear(); 
    inlierIndices.clear(); 
    outlierIndices.clear(); 
    Samples::const_iterator iter = samples.begin(); 
    int i = 0; 
    for (; iter != samples.end(); ++iter, ++i) { 
    if (isInlier(model, *iter)) { 
    inliers.push_back(*iter); 
    inlierIndices.push_back(i); 

    else { 
    outliers.push_back(*iter); 
    outlierIndices.push_back(i); 




    int calcInliersNumber(const M& model, const Samples& samples) const { 
    int n = 0; 
    for (int i = 0; i < (int)samples.size(); i++) { 
    if (isInlier(model, samples[i])) ++n; 

    return n; 


    bool isInlier(const M& model, const S& sample) const { 
    return (compute(model, sample) <= threshold); 


    private: 
    double threshold; 
    }; 

    /** generic RANSAC framework 
    * make use of a model calculator and an error calculator 
    * M is the model type, need to support copy assignment operator and default constructor 
    * S is the sample type. 

    * Interface: 
    * M compute(samples); input a set of samples, output a model. 
    * after compute, inliers and outliers can be retrieved 

    */ 
    template <typename M, typename S> 
    class Ransac : public GenericModelCalculator<M, S> { 
    public: 
    typedef std::vector<S> Samples; 

    /** Constructor 

    * @param pmc a GenericModelCalculator object 
    * @param modelSampleSize how much samples are used to fit a model 
    * @param pec a GenericErrorCaclculator object 
    */ 
    Ransac(GenericModelCalculator<M, S>* pmc, int modelSampleSize, GenericErrorCaclculator<M, S>* pec) { 
    this->pmc = pmc; 
    this->modelSampleSize = modelSampleSize; 
    this->pec = pec; 
    this->maxSampleCount = 500; 
    this->minInliersNum =  1000000

    this->verbose = false; 


    const GenericErrorCaclculator<M, S>* getErrorCalculator() const { return pec; } 

    virtual ~Ransac() { 
    delete pmc; 
    delete pec; 


    void setMaxSampleCount(int n) { 
    this->maxSampleCount = n; 


    void setMinInliersNum(int n) { 
    this->minInliersNum = n; 


    virtual bool hasThreshold() const { return true; } 

    virtual void setThreshold(double threshold) { 
    pec->setThreshold(threshold); 


    virtual double getThreshold() const { 
    return pec->getThreshold(); 


    public: 
    /** Given samples, compute a model that has most inliers. Assume the samples size is larger or equal than model sample size 
    * inliers, outliers, inlierIndices and outlierIndices are stored 

    */ 
    M compute(const Samples& samples) { 
    clear(); 

    int pointsNumber = (int)samples.size(); 

    assert(pointsNumber >= modelSampleSize); 

    long N =  100000
    int sampleCount = 0; 
    RANSAC_SamplesNumber ransac(modelSampleSize); 

    M bestModel; 
    int maxInliersNumber = 0; 

    bool stop = false; 
    while (sampleCount < N && sampleCount < maxSampleCount && !stop) { 

    Samples nsamples; 
    randomlySampleN(samples, nsamples, modelSampleSize); 

    M sampleModel = pmc->compute(nsamples); 
    if (maxInliersNumber == 0) bestModel = sampleModel; //init bestModel 

    int inliersNumber = pec->calcInliersNumber(sampleModel, samples); 
    if (verbose) std::cout<<"inliers number: "<<inliersNumber<<std::endl; 

    if (inliersNumber > maxInliersNumber) { 
    bestModel = sampleModel; 
    maxInliersNumber = inliersNumber; 
    N = ransac.calcN(inliersNumber, pointsNumber); 
    if (maxInliersNumber > minInliersNum) stop = true; 


    if (verbose) std::cout<<"N: "<<N<<std::endl; 

    sampleCount ++; 


    if (verbose) std::cout<<"sampleCount: "<<sampleCount<<std::endl; 

    finalModel = computeUntilConverge(bestModel, maxInliersNumber, samples); 

    pec->classify(finalModel, samples, inlierIndices, outlierIndices, inliers, outliers); 

    inliersRate = (double)inliers.size() / samples.size(); 

    return finalModel; 


    const Samples& getInliers() const { return inliers; } 
    const Samples& getOutliers() const { return outliers; } 

    const std::vector<int>& getInlierIndices() const { return inlierIndices; } 
    const std::vector<int>& getOutlierIndices() const { return outlierIndices; } 

    double getInliersAverageError() const { 
    return pec->computeAverage(finalModel, inliers); 


    double getInliersRate() const { 
    return inliersRate; 


    void setVerbose(bool v) { 
    verbose = v; 


    private: 
    void randomlySampleN(const Samples& samples, Samples& nsamples, int sampleSize) { 
    std::vector<int> is = ranis((int)samples.size(), sampleSize); 
    for (int i = 0; i < sampleSize; i++) { 
    nsamples.push_back(samples[is[i]]); 



    /** from initial model, iterate to find the best model. 

    */ 
    M computeUntilConverge(M initModel, int initInliersNum, const Samples& samples) { 
    if (verbose) { 
    std::cout<<"iterate until converge...."<<std::endl; 
    std::cout<<"init inliers number: "<<initInliersNum<<std::endl; 


    M bestModel = initModel; 
    M newModel = initModel; 

    int lastInliersNum = initInliersNum; 

    Samples newInliers, newOutliers; 
    pec->classify(initModel, samples, newInliers, newOutliers); 
    double lastInlierAverageError = pec->computeAverage(initModel, newInliers); 

    if (verbose) std::cout<<"init inlier average error: "<<lastInlierAverageError<<std::endl; 

    while (true && (int)newInliers.size() >= modelSampleSize) { 

    //update new model with new inliers, the new model does not necessarily have more inliers 
    newModel = pmc->compute(newInliers); 

    pec->classify(newModel, samples, newInliers, newOutliers); 

    int newInliersNum = (int)newInliers.size(); 
    double newInlierAverageError = pec->computeAverage(newModel, newInliers); 

    if (verbose) { 
    std::cout<<"new inliers number: "<<newInliersNum<<std::endl; 
    std::cout<<"new inlier average error: "<<newInlierAverageError<<std::endl; 

    if (newInliersNum < lastInliersNum) break; 
    if (newInliersNum == lastInliersNum && newInlierAverageError >= lastInlierAverageError) break; 

    //update best model with the model has more inliers 
    bestModel = newModel; 

    lastInliersNum = newInliersNum; 
    lastInlierAverageError = newInlierAverageError; 


    return bestModel; 


    void clear() { 
    inliers.clear(); 
    outliers.clear(); 
    inlierIndices.clear(); 
    outlierIndices.clear(); 


    private: 
    GenericModelCalculator<M, S>* pmc; 
    GenericErrorCaclculator<M, S>* pec; 
    int modelSampleSize; 

    int maxSampleCount; 
    int minInliersNum; 

    M finalModel; 

    Samples inliers; 
    Samples outliers; 

    std::vector<int> inlierIndices; 
    std::vector<int> outlierIndices; 

    double inliersRate; 

    private: 
    bool verbose; 

    }; 


    #endif // FVISION_RANSAC_H_ 

    实例2

    #include <math.h> 
    #include "LineParamEstimator.h" 

    LineParamEstimator::LineParamEstimator(double delta) : m_deltaSquared(delta*delta) {} 
    /*****************************************************************************/ 
    /* 
    * Compute the line parameters [n_x,n_y,a_x,a_y] 
    * 通过输入的两点来确定所在直线,采用法线向量的方式来表示,以兼容平行或垂直的情况 
    * 其中n_x,n_y为归一化后,与原点构成的法线向量,a_x,a_y为直线上任意一点 
    */ 
    void LineParamEstimator::estimate(std::vector<Point2D *> &data, 
    std::vector<double> ¶meters) 

    parameters.clear(); 
    if(data.size()<2) 
    return; 
    double nx = data[1]->y - data[0]->y; 
    double ny = data[0]->x - data[1]->x;// 原始直线的斜率为K,则法线的斜率为-1/k 
    double norm = sqrt(nx*nx + ny*ny); 

    parameters.push_back(nx/norm); 
    parameters.push_back(ny/norm); 
    parameters.push_back(data[0]->x); 
    parameters.push_back(data[0]->y); 

    /*****************************************************************************/ 
    /* 
    * Compute the line parameters [n_x,n_y,a_x,a_y] 
    * 使用最小二乘法,从输入点中拟合出确定直线模型的所需参量 
    */ 
    void LineParamEstimator::leastSquaresEstimate(std::vector<Point2D *> &data, 
    std::vector<double> ¶meters) 

    double meanX, meanY, nx, ny, norm; 
    double covMat11, covMat12, covMat21, covMat22; // The entries of the symmetric covarinace matrix 
    int i, dataSize = data.size(); 

    parameters.clear(); 
    if(data.size()<2) 
    return; 

    meanX = meanY = 0.0; 
    covMat11 = covMat12 = covMat21 = covMat22 = 0; 
    for(i=0; i<dataSize; i++) { 
    meanX +=data[i]->x; 
    meanY +=data[i]->y; 

    covMat11 +=data[i]->x * data[i]->x; 
    covMat12 +=data[i]->x * data[i]->y; 
    covMat22 +=data[i]->y * data[i]->y; 


    meanX/=dataSize; 
    meanY/=dataSize; 

    covMat11 -= dataSize*meanX*meanX; 
    covMat12 -= dataSize*meanX*meanY; 
    covMat22 -= dataSize*meanY*meanY; 
    covMat21 = covMat12; 

    if(covMat11<1e-12) { 
    nx = 1.0; 
    ny = 0.0; 

    else { //lamda1 is the largest eigen-value of the covariance matrix 
    //and is used to compute the eigne-vector corresponding to the smallest 
    //eigenvalue, which isn't computed explicitly. 
    double lamda1 = (covMat11 + covMat22 + sqrt((covMat11-covMat22)*(covMat11-covMat22) + 4*covMat12*covMat12)) / 2.0; 
    nx = -covMat12; 
    ny = lamda1 - covMat22; 
    norm = sqrt(nx*nx + ny*ny); 
    nx/=norm; 
    ny/=norm; 

    parameters.push_back(nx); 
    parameters.push_back(ny); 
    parameters.push_back(meanX); 
    parameters.push_back(meanY); 

    /*****************************************************************************/ 
    /* 
    * Given the line parameters [n_x,n_y,a_x,a_y] check if 
    * [n_x, n_y] dot [data.x-a_x, data.y-a_y] < m_delta 
    * 通过与已知法线的点乘结果,确定待测点与已知直线的匹配程度;结果越小则越符合,为 
    * 零则表明点在直线上 
    */ 
    bool LineParamEstimator::agree(std::vector<double> ¶meters, Point2D &data) 

    double signedDistance = parameters[0]*(data.x-parameters[2]) + parameters[1]*(data.y-parameters[3]); 
    return ((signedDistance*signedDistance) < m_deltaSquared); 

    #include <math.h>
    #include "LineParamEstimator.h"
    
    LineParamEstimator::LineParamEstimator(double delta) : m_deltaSquared(delta*delta) {}
    /*****************************************************************************/
    /*
     * Compute the line parameters  [n_x,n_y,a_x,a_y]
     * 通过输入的两点来确定所在直线,采用法线向量的方式来表示,以兼容平行或垂直的情况
     * 其中n_x,n_y为归一化后,与原点构成的法线向量,a_x,a_y为直线上任意一点
     */
    void LineParamEstimator::estimate(std::vector<Point2D *> &data, 
    																	std::vector<double> ¶meters)
    {
    	parameters.clear();
    	if(data.size()<2)
    		return;
    	double nx = data[1]->y - data[0]->y;
    	double ny = data[0]->x - data[1]->x;// 原始直线的斜率为K,则法线的斜率为-1/k
    	double norm = sqrt(nx*nx + ny*ny);
    	
    	parameters.push_back(nx/norm);
    	parameters.push_back(ny/norm);
    	parameters.push_back(data[0]->x);
    	parameters.push_back(data[0]->y);		
    }
    /*****************************************************************************/
    /*
     * Compute the line parameters  [n_x,n_y,a_x,a_y]
     * 使用最小二乘法,从输入点中拟合出确定直线模型的所需参量
     */
    void LineParamEstimator::leastSquaresEstimate(std::vector<Point2D *> &data, 
    																							std::vector<double> ¶meters)
    {
    	double meanX, meanY, nx, ny, norm;
    	double covMat11, covMat12, covMat21, covMat22; // The entries of the symmetric covarinace matrix
    	int i, dataSize = data.size();
    
    	parameters.clear();
    	if(data.size()<2)
    		return;
    
    	meanX = meanY = 0.0;
    	covMat11 = covMat12 = covMat21 = covMat22 = 0;
    	for(i=0; i<dataSize; i++) {
    		meanX +=data[i]->x;
    		meanY +=data[i]->y;
    
    		covMat11	+=data[i]->x * data[i]->x;
    		covMat12	+=data[i]->x * data[i]->y;
    		covMat22	+=data[i]->y * data[i]->y;
    	}
    
    	meanX/=dataSize;
    	meanY/=dataSize;
    
    	covMat11 -= dataSize*meanX*meanX;
            covMat12 -= dataSize*meanX*meanY;
    	covMat22 -= dataSize*meanY*meanY;
    	covMat21 = covMat12;
    
    	if(covMat11<1e-12) {
    		nx = 1.0;
    	        ny = 0.0;
    	}
    	else {	    //lamda1 is the largest eigen-value of the covariance matrix 
    	           //and is used to compute the eigne-vector corresponding to the smallest
    	           //eigenvalue, which isn't computed explicitly.
    		double lamda1 = (covMat11 + covMat22 + sqrt((covMat11-covMat22)*(covMat11-covMat22) + 4*covMat12*covMat12)) / 2.0;
    		nx = -covMat12;
    		ny = lamda1 - covMat22;
    		norm = sqrt(nx*nx + ny*ny);
    		nx/=norm;
    		ny/=norm;
    	}
    	parameters.push_back(nx);
    	parameters.push_back(ny);
    	parameters.push_back(meanX);
    	parameters.push_back(meanY);
    }
    /*****************************************************************************/
    /*
     * Given the line parameters  [n_x,n_y,a_x,a_y] check if
     * [n_x, n_y] dot [data.x-a_x, data.y-a_y] < m_delta
     * 通过与已知法线的点乘结果,确定待测点与已知直线的匹配程度;结果越小则越符合,为
     * 零则表明点在直线上
     */
    bool LineParamEstimator::agree(std::vector<double> ¶meters, Point2D &data)
    {
    	double signedDistance = parameters[0]*(data.x-parameters[2]) + parameters[1]*(data.y-parameters[3]); 
    	return ((signedDistance*signedDistance) < m_deltaSquared);
    }


    RANSAC寻找匹配的代码如下:

    [cpp]  
    view plain copy print ?





    /*****************************************************************************/ 
    template<class T, class S> 
    double Ransac<T,S>::compute(std::vector<S> ¶meters, 
    ParameterEsitmator<T,S> *paramEstimator , 
    std::vector<T> &data, 
    int numForEstimate) 

    std::vector<T *> leastSquaresEstimateData; 
    int numDataObjects = data.size(); 
    int numVotesForBest = -1; 
    int *arr = new int[numForEstimate];// numForEstimate表示拟合模型所需要的最少点数,对本例的直线来说,该值为2 
    short *curVotes = new short[numDataObjects]; //one if data[i] agrees with the current model, otherwise zero 
    short *bestVotes = new short[numDataObjects]; //one if data[i] agrees with the best model, otherwise zero 


    //there are less data objects than the minimum required for an exact fit 
    if(numDataObjects < numForEstimate) 
    return 0; 
    // 计算所有可能的直线,寻找其中误差最小的解。对于100点的直线拟合来说,大约需要100*99*0.5=4950次运算,复杂度无疑是庞大的。一般采用随机选取子集的方式。 
    computeAllChoices(paramEstimator,data,numForEstimate, 
    bestVotes, curVotes, numVotesForBest, 0, data.size(), numForEstimate, 0, arr); 

    //compute the least squares estimate using the largest sub set 
    for(int j=0; j<numDataObjects; j++) { 
    if(bestVotes[j]) 
    leastSquaresEstimateData.push_back(&(data[j])); 

    // 对局内点再次用最小二乘法拟合出模型 
    paramEstimator->leastSquaresEstimate(leastSquaresEstimateData,parameters); 

    delete [] arr; 
    delete [] bestVotes; 
    delete [] curVotes; 

    return (double)leastSquaresEstimateData.size()/(double)numDataObjects; 
    }

     

    相关文章推荐

    大型网站的灵魂——性能
    http://www.lai18.com/content/318116.html

    RSA,DSA等加解密算法介绍
    http://www.lai18.com/content/321599.html

    权限项目总结(一)权限设计
    http://www.lai18.com/content/409780.html

    快速排序里的学问:随机化快排
    http://www.lai18.com/content/423495.html

    快速排序里的学问:枢纽元选择与算法效率
    http://www.lai18.com/content/423496.html

    快速排序里的学问:霍尔快排的实现
    http://www.lai18.com/content/423497.html

    快速排序里的学问:霍尔与快速排序
    http://www.lai18.com/content/423498.html

    快速排序里的学问:快速排序的过程
    http://www.lai18.com/content/423499.html

    快速排序里的学问:信息熵
    http://www.lai18.com/content/423501.html

    快速排序里的学问:再看看称球问题
    http://www.lai18.com/content/423502.html

    快速排序里的学问:从猜数字开始
    http://www.lai18.com/content/423503.html

    漫谈递归:从汇编看尾递归的优化
    http://www.lai18.com/content/423517.html

    漫谈递归:PHP里的尾递归及其优化
    http://www.lai18.com/content/423518.html

    漫谈递归:补充一些Continuation的知识
    http://www.lai18.com/content/423519.html

    漫谈递归:尾递归与CPS
    http://www.lai18.com/content/423520.html

    漫谈递归:从斐波那契开始了解尾递归
    http://www.lai18.com/content/423521.html

    展开全文
  • 7-2 OpenCV算法解析

    2020-09-11 20:13:25
    7-2 OpenCV算法解析目录1 OpenCVOpenCV大坑之BGROpenCV 性能OpenCV常见算法2 最小二乘法2.1 线性回归什么是线性回归?2.2 线性回归--最小二乘法(Least Square Method)2.2.1 最小二乘法2.2.2 三种范数:2.2.3 最小...
  • OpenCV 图像采样 插值 几何变换

    千次阅读 2010-10-09 14:11:00
    OpenCV 图像采样 插值 几何变换
  • 参考: https://blog.csdn.net/jay463261929/article/details/53818611
  • opencv中插值算法详解

    千次阅读 多人点赞 2020-03-08 23:00:21
    做图像处理的同学应该经常都会用到图像的缩放,我们都知道图片存储的时候其实就是一个矩阵,所以在对图像进行缩放操作的时候,也就是在对矩阵进行操作,如果想要将图片放大,这里我们就需要用到过采样算法来扩大矩阵...
  • OpenCV之LBP算法学习

    千次阅读 2017-12-20 16:22:08
    目录 目录 前言 LBP算法概述 LBP算法原理 原始LBP特征描述及计算方法 ...LBP算法概述 ...LBP指局部二值模式,英文全称:...LBP常应用于人脸识别和目标检测中,在OpenCV中有使用LBP特征进行人脸识别的接口,也有用LBP特征
  •  在计算视觉领域,尺度空间被象征性的表述为一个图像金字塔,向下降采样一般用高斯金字塔。其中,输入图像函数反复与高斯函数的核卷积并反复对其进行二次抽样,这种方法主要用于SIFT算法的实现,但每层图像依赖于
  • Python-opencv3 SIFT算法做特征匹配

    万次阅读 热门讨论 2018-04-13 17:51:39
    最近接触一个项目:根据设计师定出的psd格式文件(photoshop),生成不同尺寸的海报。这里面牵扯到了尺度不变而对特征做变换的问题...SIFT(Scale-invariant feature transform),也叫尺度不变特征变换算法,是Davi...
  • #include #include #include using namespace std;...//基于旧版本的opencv的LBP算法opencv1.0 // 3 x 3 矩阵如下所示 // [ 1, 2, 3] // [ 8, ij,4] // [ 7, 6, 5] void LBP(IplImage *src, IplImage *dst)
  • 在网上找了好久都没找到基于opencv的金字塔模板匹配算法代码,我就自己把金字塔和模板匹配的代码结合了一下,代码基于opencv2.48.
  • vs2019使用opencv实现ViBe算法 参考代码网址:https://github.com/upcAutoLang/BackgroundSplit-OpenCV/tree/master/src/ViBe 链接 ViBe参考文章:https://www.jianshu.com/p/48baa72c6e5f 链接 (1)Vibe.h #...
  • opencv之SURF算法原理及关键点检测

    万次阅读 多人点赞 2017-04-18 10:51:54
    1.概述在基础篇里面讲模板匹配的时候已经介绍过,图像匹配主要有...SIFT是一种鲁棒性好的尺度不变特征描述方法,但SIFT算法计算数据量大、时间复杂度高、算法耗时长。针对上述缺点许多研究者对SIFT算法做了不同的改进,
  • OpenCV图像插值算法 1.1 简介   在图像处理中,平移变换、旋转变换以及放缩变换是一些基础且常用的操作。这些几何变换并不改变图象的象素值,只是在图象平面上进行象素的重新排列。在一幅输入图象[u,v][u,v][u,...
  • OpenCV图像滤波算法总结(Python)

    千次阅读 2019-08-05 23:30:36
    一.图像的滤波处理 ...过滤可以移除图像中的噪音、提取感兴趣的可视特征、允许图像重采样等等。 **频域分析 :**将图像分成从低频到高频的不同部分。低频对应图像强度变化小的区域,而高频是...
  • 基于OpenCV的SIFT算法的实现

    千次阅读 2014-03-30 15:34:01
    首先要在你的电脑上安装OpenCV的环境,因为实现SIFT的算法我们用到的是OpenCV的库文件中的有关函数来完成的,因此,首先必须安装OpenCV到你的编译环境中。住:我的电脑中安装的是VS2010,因此我会着重的讲一下在VS...
  • OpenCV 的LBPH算法原理

    千次阅读 2019-04-04 16:38:20
    1 背景及理论基础 人脸识别是指将一个需要识别的人脸和人脸库中的某个人脸对应起来(类似于指纹识别),目的是完成识别功能,该术语需要和...从OpenCV2.4开始,加入了新的类FaceRecognizer,该类用于人脸识别,...
  • 本文实例讲述了Python基于opencv的图像压缩算法。分享给大家供大家参考,具体如下: 插值方法: CV_INTER_NN – 最近邻插值, CV_INTER_LINEAR – 双线性插值 (缺省使用) CV_INTER_AREA – 使用象素关系重采样。当图像...
  • Opencv3中SURF算法学习

    千次阅读 2021-07-20 22:35:24
    Opencv学习之SURF算法; harr特征以及积分图像; Hessian矩阵; sift算法; 图像金字塔、高斯金字塔、差分金字塔(DOG金字塔)、尺度空间; 非极大值抑制;
  • Opencv学习之SURF算法

    万次阅读 多人点赞 2017-06-27 14:52:36
    Opencv学习之SURF算法 SURF(加速版的具有鲁棒性的特征,SpeededUp Robust Features),SURF是尺度不变特征变换算法(SIFT算法)的加速版。SURF最大的特征在于采用了harr特征以及积分图像的概念。 SURF原理: ...
  • AI学习笔记之图像滤波器、OpenCV算法解析、深度学习与神经网络图像滤波器图像噪声噪声的产生信噪比高斯噪声椒盐噪声其他噪声图像滤波滤波的目的滤波的要求各种滤波器均值滤波中值滤波最大最小值滤波双边滤波限幅滤波...
  • 1.图像上采样和降采样 (1)图像金字塔概念 我们在图像处理中常常会调整图像大小,最常见的就是放大(zoomin)和缩小(zoomout),尽管几何变换也可以实现图像放大和缩小,但是这里我们介绍图像金字塔。 一个图像...
  • 关于OpenCV的上采样和下采样

    千次阅读 2018-12-18 10:13:05
    采样(subsampled)(或称为缩小图像 或降采样(downsampled))的主要目的有两个:1、使得图像符合显示区域的大小;2、生成对应图像的缩略图(最直观的理解,所以深度学习领域对被卷积核之后的特征图叫下采样,从...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 8,043
精华内容 3,217
关键字:

opencv采样算法