图像匹配

基于像素的匹配

1、归一化积相关灰度匹配：

 

      模板图像 以窗口滚动的方式 在源图像中 扫一遍。

 

     具体运算公式如下：

 

R（i,j） = dSigmaST / (dSigmaT * dSigmaS)  对应上公式； R（i,j）=[0,1]

 

M ，N 模板大小

 

对于公式的解释：

 

                          dSimgmaST --  在 在原图（i ，j） 位置 模板图像每个像素与对应原图像素的积 的和）

 

                            dSigmaT  -- 模板图像每个像素的积 的和

 

                            dSigmaS  -- 在原图（i,j）位置，模板图像对应的原图的每个像素的积的和

                           

 

R（i,j）最大的位置就是最匹配的位置。

 

 

/* 
    函数：NormalizeGrayMatch
    功能：归一化灰度值匹配
    参数：src 原图
             template 模板
             point 匹配的位置左上角
    返回值： 0 --- 未找打匹配对象
                  其他 --- 返回匹配位置相似度
*/
double NormalizeGrayMatch(IplImage  * src , IplImage  * temp , CvPoint  & point)
{
     int nSrcWidth  = src - >width ; 
     int nSrcHeight  = src - >height;
     int nTwidth  = temp - >width ; 
     int nTheight  = temp - >height ; 


     //计算模板像素灰度值 dSigmaT
     double dSigmaT  =  0 ; 
     unsigned  char piexl  =  0 ; 
     for ( int i  =  0 ; i  < nTheight ;  ++ i)
    {
         for( int j  =  0 ; j < nTwidth ;  ++j)
        {
            piexl  =  *(temp - >imageData  + i  * temp - >widthStep  + j ) ; 
            dSigmaT  += ( double) piexl  * piexl ; 
        }
    }

     double R  =  0 ; 
     double dSigmaS ;
     double dSigmaST ;
     unsigned  char piexlT,piexlS;
     double dMaxR  =  - 1 ; 
     int nMaxHeight  =  - 1 ;
     int nMaxWidht  =  - 1 ;
     for ( int i  =  0 ; i  < nSrcHeight  - nTheight  + 1 ;  ++ i )
    {
         for ( int j  =  0 ; j  < nSrcWidth  - nTwidth  +  1 ;  ++ j)
        {

            dSigmaST  =  0 ;
            dSigmaS  =  0 ;
             // 计算dSigmaST dSigmaS
             for ( int k  =  0 ; k  < nTheight ;  ++ k)
            {
                 for ( int l  =  0 ; l  < nTwidth ;  ++ l)
                {
                     //模板像素
                    piexlT  = *( temp - >imageData  + k  * temp - >widthStep  + l ); 
                     //源图像像素
                    piexlS  =  *(src - >imageData  + (k  + i )  * src - >widthStep  + (l  + j) ) ; 

                    dSigmaST  += ( double) piexlS  * piexlT ;
                    dSigmaS  += ( double) piexlS  * piexlS ; 
                }
            }

            R  = dSigmaST  / ( sqrt(dSigmaS)  * sqrt(dSigmaT) ) ;
             if (R  > dMaxR )
            {
                nMaxWidht  = j;
                nMaxHeight  = i ;
                dMaxR =R ; 
            }

        }
    }

     if (dMaxR  !=  - 1)
    {
        point .x  = nMaxWidht ; 
        point .y  = nMaxHeight;
         return dMaxR ;
    }

     return  0 ; 
}

 

 

2. 序贯相似性检测法匹配 SSDA  -- Similarity Sequential Dectection Algorithm

 

 前面的归一化积相关匹配算法计算量很大，原因在于搜索窗口在源图像上进行滑动，每滑动一次就要做一次匹配运算，除了匹配的点外在其他匹配点做了‘无用功’，导

 

致了匹配算法的计算量上升。所以，一旦发现所在的参考位置为非匹配点，就丢弃不再计算，立刻换到新的参考点计算，可以大大加速匹配过程。

 

SSDA算法过程：

 

1.定义绝对误差：

 

 

其中：

 

 模板对应原图（子图）的灰度平均值

 

      模板的灰度平均值

 

2.取一个不变阈值

 

3.扫面原图每个像素点待匹配点对应的模板子图，根据1中的公式据算绝对误差，当累加值超过阈值时停止累加，停止此模板子图的扫描，记录带匹配点的位置和累加次数

 

4，循环 3 直到扫描完全图

 

5，累加次数最少的像素点就为最佳匹配点

 

对于序贯相似相算法有很多可以优化的地方， 比如 第三步 扫描子图像素的时候 可以用 隔行 列 扫描 ， 第二部 阈值可以改为自动阈值等等，不过这里的优化有能

 

怎么样呢 还是很慢呀，还是期待后面的算法吧。

 

/* 
        函数：SSDAMatch 
        参数：src --- 源图像 in
                  temp -- 匹配模板 in
                  point -- 匹配位置 out 
                  threshold -- 阈值 in
       返回值：累加次数
       限制： 8位灰度图
*/
int SSDAMatch(IplImage  * src , IplImage  * temp , CvPoint  & point ,  int threshold)
{
     int nSrcWidth  = src - >width ; 
     int nSrcHeight  = src - >height;
     int nTwidth  = temp - >width ; 
     int nTheight  = temp - >height ; 
     int ntempsize  = nTwidth  * nTheight ;  //模板大小

     //计算模板像素灰度值 dSigmaT
     double dSigmaT  =  0 ; 
     unsigned  char piexl  =  0 ; 
     for ( int i  =  0 ; i  < nTheight ;  ++ i)
    {
         for( int j  =  0 ; j < nTwidth ;  ++j)
        {
            piexl  =  *(temp - >imageData  + i  * temp - >widthStep  + j ) ; 
            dSigmaT  += ( double) piexl ; 
        }
    }
    dSigmaT  /=ntempsize ; 


     double dSigmaS  =  0 ; 
     double dbr  = 0;  //误差
     long lr  =  0;  //误差累积次数
     long maxR  = 0;  //最大累积次数
     int nMaxHeight  =  0 ;  //最大累计次数 对应匹配位置 (左上角)
     int nMaxWidht  =  0 ;
     for ( int i  =  0 ; i  < nSrcHeight  - nTheight  + 1 ;  ++i)
    {
         for ( int j  =  0 ; j  < nSrcWidth  - nTwidth  +  1 ;  ++ j)
        {
             //计算dSigmaS
            dSigmaS  =  0 ;
            dbr  =  0 ;
            lr  =  0 ; 
             for ( int k  =  0 ; k  < nTheight ;  ++k)
            {
                 for( int l  =  0 ; l  < nTwidth ;  ++l)
                {
                    dSigmaS  +=( unsigned  char )  *(src - >imageData  + (k  +i)  * src - >widthStep  + (l  + j)) ;
                }
            }
            dSigmaS  /=ntempsize ; 

             //计算误差 一旦超过阈值则抛弃不再计算
             for ( int k  =  0 ; k  < nTheight ;  ++ k)
            {
                 for ( int l  =  0 ; l  < nTwidth ;  ++l)
                {
                    dbr  += abs( ( unsigned  char ) *(src - >imageData  + (k +i)  * src - >widthStep  + (l +j) ) -dSigmaS  - 
                        ( unsigned  char ) *(temp - >imageData  + k  * temp - >widthStep  + l)  +dSigmaT );
                    lr  ++ ;
                     if (dbr  > = threshold) 
                         break;
                }
                 if (dbr > =threshold)
                     break;
            }

             //取达到threshold 累加最多的 位置
             if ( lr  >maxR )
            {
                maxR  = lr ; 
                nMaxHeight  = i;
                nMaxWidht  = j ;
            }
        }
    }
        point.x  = nMaxWidht ; 
        point.y  = nMaxHeight ; 
         return lr;
}

 

 

基于特征的匹配

 

利用灰度信息匹配的方法主要缺陷是计算量过大，对图像的灰度变换很敏感，尤其是非线性的光照变化，此外，对目标的旋转、变形以及遮挡也比较敏感，为

 

了克服这些缺点，可以利用图像的特征进行匹配，由于图像的特征点比像素点要少很多，大大减少了匹配过程的计算量，同时，特征点的匹配度量值对位置的

 

变化比较敏感，可以大大的提高匹配的精度。而且，特征点的提取过程可以减少噪声的影响，对灰度变化、图像形变以及遮挡等有较好的适应能力。

 

3. 不变矩匹配法     TM算法 具有平移、旋转、尺寸不变性

 

 

                                                                                                                                                                                p+q>=2

 

 

 

归一化公式：

 

算法过程：计算 分别计算模板和原图的7个不变矩 ，根据归一化公式得出相似度。

 

 

double momentMatch( unsigned  char  * src ,  unsigned  char  * temp , int nwidth ,  int nheight , int nwidthstep ) 
{
     //原图和模板重心矩
     int nSBarycenterX , nSBarycenterY;
     int nTBasrycenterX,nTBarycenterY;
    CalBarycenter(src , nwidth , nheight , nwidthstep ,nSBarycenterX , nSBarycenterY ) ; 
    CalBarycenter(temp , nwidth , nheight ,nwidthstep , nTBasrycenterX , nTBarycenterY);

     //原图和模板二阶三阶规格化中心矩
     double Su00 ,Su02 ,Su20 ,Su11 ,Su30 ,Su12 ,Su21 ,Su03;
     double Tu00 ,Tu02 ,Tu20 ,Tu11 ,Tu30 ,Tu12 ,Tu21 ,Tu03;
    Su00  = CalCenterMoment(src , nwidth ,nheight ,nwidthstep ,nSBarycenterX,nSBarycenterY, 0, 0) ; 
    Su02  = CalCenterMoment(src , nwidth ,nheight , nwidthstep , nSBarycenterX,nSBarycenterY, 0, 2) /pow(Su00 ,  2);
    Su20  = CalCenterMoment(src , nwidth , nheight , nwidthstep , nSBarycenterX,nSBarycenterY, 2, 0) /pow(Su00 ,  2);
    Su11  = CalCenterMoment(src , nwidth , nheight , nwidthstep , nSBarycenterX,nSBarycenterY, 1, 1) /pow(Su00, 2);
    Su30  = CalCenterMoment(src , nwidth , nheight ,nwidthstep , nSBarycenterX,nSBarycenterY , 3 , 0 ) /pow(Su00 ,  2. 5);
    Su12  = CalCenterMoment(src , nwidth , nheight , nwidthstep , nSBarycenterX,nSBarycenterY , 1 , 2 ) /pow(Su00 ,  2. 5);
    Su21  = CalCenterMoment(src , nwidth , nheight , nwidthstep , nSBarycenterX,nSBarycenterY ,  2 , 1 ) /pow(Su00 ,  2. 5);
    Su03  = CalCenterMoment(src , nwidth , nheight , nwidthstep , nSBarycenterX,nSBarycenterY , 0, 3) /pow(Su00 ,  2. 5);

    Tu00  = CalCenterMoment(temp , nwidth ,nheight ,nwidthstep ,nTBasrycenterX , nTBarycenterY, 0, 0) ; 
    Tu02  = CalCenterMoment(temp , nwidth ,nheight , nwidthstep , nTBasrycenterX , nTBarycenterY, 0, 2) /pow(Su00 ,  2);
    Tu20  = CalCenterMoment(temp , nwidth , nheight , nwidthstep , nTBasrycenterX , nTBarycenterY, 2, 0) /pow(Su00 ,  2);
    Tu11  = CalCenterMoment(temp , nwidth , nheight , nwidthstep , nTBasrycenterX , nTBarycenterY, 1, 1) /pow(Su00, 2);
    Tu30  = CalCenterMoment(temp , nwidth , nheight ,nwidthstep , nTBasrycenterX , nTBarycenterY , 3 , 0 ) /pow(Su00 ,  2. 5);
    Tu12  = CalCenterMoment(temp , nwidth , nheight , nwidthstep , nTBasrycenterX , nTBarycenterY , 1 , 2 ) /pow(Su00 ,  2. 5);
    Tu21  = CalCenterMoment(temp , nwidth , nheight , nwidthstep , nTBasrycenterX , nTBarycenterY ,  2 , 1 ) /pow(Su00 ,  2. 5);
    Tu03  = CalCenterMoment(temp , nwidth , nheight , nwidthstep , nTBasrycenterX , nTBarycenterY, 0, 3) /pow(Su00 ,  2. 5);


     //原图和模板不变矩
     double Sa[ 7] , Ta[ 7];
    Sa[ 0]  = Su02  + Su20 ; 
    Sa[ 1]  = pow(Su20  - Su02 ,  2)  +  4  * pow(Su11 ,  2 );
    Sa[ 2]  = pow(Su30  -  3 *Su12 ,  2)  + pow( 3  * Su12  -Su03 , 2);
    Sa[ 3]  = pow(Su30 +Su12 ,  2)  + pow(Su21  + Su03 , 2);
    Sa[ 4]  = (Su30  - 3 *Su12)  * (Su30  + Su12 )  * (pow(Su30  + Su12 , 2)  -  3 *pow(Su21  + Su03 , 2))  + 
                ( 3  * Su21  -Su03)  *(Su21  + Su03)  *( 3 * pow(Su03  + Su12 ,  2)  - pow(Su21  + Su03 , 2));
    Sa[ 5]  = (Su20  - Su02) *(pow(Su30  +Su12 , 2)  - pow(Su21  + Su03 , 2))  +  4 *Su11  *(Su30  + Su12) *(Su21  +Su03);
    Sa[ 6]  = ( 3 *Su21  - Su03)  *(Su30  + Su12)  *(pow(Su30  + Su12 , 2)  -  3  *pow(Su21  + Su03 , 2))  + (Su30  - 3 *Su12)  *(Su21  +Su03)  *
                ( 3  *pow(Su30  + Su12 , 2)  - pow(Su21  + Su03 , 2));

    Ta[ 0]  = Tu02  + Tu20 ; 
    Ta[ 1]  = pow(Tu20  - Tu02 ,  2)  +  4  * pow(Tu11 ,  2 );
    Ta[ 2]  = pow(Tu30  -  3 *Tu12 ,  2)  + pow( 3  * Tu12  -Tu03 , 2);
    Ta[ 3]  = pow(Tu30 +Tu12 ,  2)  + pow(Tu21  + Tu03 , 2);
    Ta[ 4]  = (Tu30  - 3 *Tu12)  * (Tu30  + Tu12 )  * (pow(Tu30  + Su12 , 2)  -  3 *pow(Tu21  + Tu03 , 2))  + 
        ( 3  * Tu21  -Tu03)  *(Tu21  + Tu03)  *( 3 * pow(Tu03  + Tu12 ,  2)  - pow(Tu21  + Tu03 , 2));
    Ta[ 5]  = (Tu20  - Tu02) *(pow(Tu30  +Tu12 , 2)  - pow(Tu21  + Tu03 , 2))  +  4 *Tu11  *(Tu30  + Tu12) *(Tu21  +Tu03);
    Ta[ 6]  = ( 3 *Tu21  - Tu03)  *(Tu30  + Tu12)  *(pow(Tu30  + Tu12 , 2)  -  3  *pow(Tu21  + Tu03 , 2))  + (Tu30  - 3 *Tu12)  *(Tu21  +Tu03)  *
        ( 3  *pow(Tu30  + Tu12 , 2)  - pow(Tu21  + Tu03 , 2));


     double r  =  0;
     double dSigmaST  = 0;
     double dSigmaS  =  0;
     double dSigmaT  =  0 ;
     for ( int i  =  0 ; i  <  7 ;  ++i)
    {
        dSigmaST  += Ta[i]  * Sa[i] ;
        dSigmaS  +=pow(Sa[i] , 2);
        dSigmaT  +=pow(Ta[i] , 2);
    }
     return r  = dSigmaST  / sqrt( dSigmaS  * dSigmaT) ;
}

/*
    函数：CalBarycenter
    功能：计算重心矩
    参数：pdata -- 图像数据 in
              nwidth -- 宽 in
              nheight -- 高 in
              nwidthstep -- 步长 in
              nBarycenterX -- 重心坐标 out
              nBarycenterY   
*/
void CalBarycenter( unsigned  char  * pdata ,  int nwidth ,  int nheight , int nwidthstep ,  int  &nBarycenterX ,  int  &nBarycenterY)
{
         double m00 , m01 ,m10; 

        m00  =  0 ;
        m01  =  0 ;
        m10  =  0 ;
         for ( int i  =  0 ; i  < nheight ;  ++i)
        {
             for ( int j  =  0 ; j  < nwidth ;  ++ j)
            {
                m00  += *(pdata  + i  * nwidthstep  + j) ; 
                m01  += *(pdata  + i  * nwidthstep  + j)  * j ;
                m10  += *(pdata  + i  * nwidthstep  + j)  * i ; 
            }
        }

        nBarycenterX  =( int) (m10  / m00  + 0. 5);
        nBarycenterY  = ( int)(m01 / m00  +  0. 5);
}

/*
    函数：CalCenterMoment
    功能：计算中心矩
    参数：pdata --- 图像数据 in
              nwidth -- 宽 in
              nheight -- 高 in
              nwidthstep -- 步长 in
              nBarycenterX -- 重心矩 in
              nBarycenterY 
              ip -- 阶数 in
              jq
      返回值：中心距值
*/
double CalCenterMoment( unsigned  char  * pdata ,  int nwidth ,  int nheight , int nwidthstep , 
                                     double nBarycenterX ,  double nBarycenterY, int ip, int jq)
{
     double Upq  =  0 ; 

     for ( int i  =  0 ; i  < nheight ;  ++i)
    {
         for ( int j = 0; j  <nwidth ;  ++ j)
        {
            Upq  += *(pdata  + i  * nwidthstep  + j)  + pow(j  -nBarycenterX , ip)  * pow(i  - nBarycenterY , jq) ;
        }
    }

     return Upq ; 
}

 

 

4.距离变换匹配法

 

距离变换是一种常见的二值图像处理算法，用来计算图像中任意位置到最近边缘点的距离

 

 

欧几米德空间距离：

 

其中a=（x1 , y1） b =(x2 ,y2)

 

设R 为二维图像空间集合 S 为R中的边缘点集合 R中任意一点 r，的距离变换为

 

 

如果该点是边缘点则 距离就为0 ，如果不是边缘点则找与之最近的边缘点距离，就需要做最近邻域搜索这样计算量很大。

 

所以有一种近似方法作为替代 即查看点8邻域内的情况如下表所示：

 

 

8邻域内 P 到各个点的距离 ， 如果8邻域内找不到 边缘点 则距离为1；

 

g(x) = 0     if=0

          0.3   x=1

          0.7   x=sqrt(2)

          1      x>sqrt(2)

 

 

对于两幅二值图像A（模板图），B（模板子图） 的匹配误差度量为：

 

 

A，B 为图像中的边缘点集合。   a,b分别为A，B中的任意一点   Na Nb为A,B的点个数。

 

Pmatch = [0,1] ；这个公式表示：在模板对应图中如果遍历所有边缘点，累加 边缘点位置对应的模板图的该点的距离变换。

 

可知如果两个图像完全一样 则 Pmatch =0 ；

 

 

 

算法过程： 求模板图的距离变换 ， 模板图边缘点个数

 

                    在待匹配图中 滑动窗口（大小为模板图大小） ，在模板对应图中如果遍历所有边缘点，累加 边缘点位置对应的模板图的该点的距离变换

 

                    dSigmaST  ， Pmatch = （dSigmaST + | Na - Nb | ） /(Na + Nb) ;

 

                    Pmatch最小时为最佳匹配点.

               

            

/* 
    函数：DisMatch
    功能：距离变换匹配 
    参数：src -- 原图 in
              temp -- 模板图 in
              point -- 最佳匹配点 out
    返回值：匹配误差
*/

double DisMatch(IplImage  * src,IplImage  *temp,CvPoint  &point )
{
     int nSwidth  = src - >width ; 
     int nSheight  =src - >height;
     int nSwidthstep  = src - >widthStep;

     int nTwidth  = temp - >width;
     int nTheight  = temp - >height;
     int nTwidthstep  = temp - >widthStep;

     unsigned  char piexl;

     //将模板图像长宽加2 方便8邻域计算
     unsigned  char  * pTdata  = ( unsigned  char  *)malloc( (nTwidthstep + 2 )  *(nTheight  + 2) );
    memset(pTdata ,  0, (nTwidthstep + 2 )  *(nTheight  + 2)  *  sizeof( unsigned  char));
     for ( int n  =  1 ; n  < nTheight  + 1;  ++n)
        memcpy( pTdata  +n  *(nTwidthstep  + 2) + 1 , temp - >imageData  + n  * nTwidthstep , nTwidthstep  * sizeof( unsigned  char) );

     //计算模板距离变换
     double  * pTDist  = ( double  *)malloc(nTwidth  * nTheight  *  sizeof( double));
     //模板T边缘点个数
     int Nb  = 0;
     //8邻域
     unsigned  char u11,u12,u13,u21,u23,u31,u32,u33;
     for( int i  =  1 ; i  < nTheight  + 1 ;  ++i)
    {
         for ( int j  =  1 ; j  <nTwidth  + 1 ; ++j)
        {
            piexl  =  *( pTdata  + i  * (nTwidthstep + 2)  +j);
             //如果该点是边缘点 dist=0；
             if (piexl  == 255)
            {
                pTDist[(i - 1)  *nTwidth  + (j - 1)]  = 0;
                Nb ++;
            }
             else
            {
                 //否则看8邻域 膨化加权
                u11  =  *(pTdata  + (i - 1) *(nTwidthstep + 2) +j - 1);
                u12  =  *(pTdata  +(i - 1) *(nTwidthstep + 2) +j);
                u13  = *(pTdata  +(i - 1) *(nTwidthstep + 2) +j + 1);
                u21  = *(pTdata  +i *(nTwidthstep + 2) +j - 1);
                u23  = *(pTdata  +i *(nTwidthstep + 2) +j + 1);
                u31  = *(pTdata  + (i + 1) *(nTwidthstep + 2) +j - 1);
                u32  = *(pTdata  +(i + 1) *(nTwidthstep + 2) +j);
                u33  = *(pTdata  +(i + 1) *(nTwidthstep + 2) +j + 1);
                 if (u12 == 255  ||u21 == 255  ||u23 == 255  ||u32 == 255)
                    pTDist[(i - 1)  *nTwidth  + (j - 1)]  = 0. 3;
                 else  if(u11 == 255  || u13 == 255  ||u31 == 255  ||u33 == 255)
                    pTDist[(i - 1)  *nTwidth  + (j - 1)]  = 0. 7;
                 else
                    pTDist[(i - 1)  *nTwidth  + (j - 1)]  = 1;
            }
        } //for j
    }  //for i

     //匹配
     double dbMatch  =  0 ;  //匹配误差
     //最小匹配误差
     double dbMinMatch  =  1;
     double dSigmaST  =  0 ;
     //边缘点个数
     int Na  =  0;
     //最佳匹配点
     int nMatchWidth , nMatchHeight;

     for( int i  =  0 ; i  <nSheight  -nTheight + 1;  ++i)
    {
         for ( int j = 0; j <nSwidth  -nTwidth + 1; ++j)
        {
            dSigmaST  =  0;
            dbMatch  =  0;
            Na  = 0;
             //模板滑动到（i，j）点对应的原图
             for ( int k = 0;k <nTheight ;  ++k)
            {
                 for ( int l = 0 ; l <nTwidth ; ++l)
                {
                    piexl  =  *(src - >imageData  + (i  +k) *nSwidthstep  +(j +l));
                     if (piexl == 255)
                    {
                        dSigmaST  +=pTDist[k  *nTwidth  +l];
                        Na ++;
                    }
                } //for l
            } //for l

             //计算匹配误差
            dbMatch  = ( dSigmaST  + abs(Na -Nb))  /(Na  +Nb);
             if (dbMatch  < dbMinMatch)
            {
                dbMinMatch  = dbMatch;
                nMatchHeight  = i ;
                nMatchWidth  = j;
            }

        }  // for j
    }  //for i

    point.x  = nMatchWidth;
    point.y  = nMatchHeight;
    free(pTDist);
    free(pTdata);

     return dbMinMatch;
}
            
            

 

5.最小均方误差匹配法

   

 最小均方误差匹配方法是利用图像中的对应 特征点，通过解特征点的变换防长来计算图像间的变换参数。

 

 基本原理： 最小均方误差匹配方法是以模板中的特征点构造矩阵X ，图像子图中的特征点构造矩阵Y ，求解矩阵X 到矩阵Y 的变换矩阵， 其中误差

 

 最小的位置为最佳匹配位置

 

 图像间的仿变换方程：

 

 

                       原点为中心旋转               平移

 

  仿变换参数为

 

 

构建图像矩阵 X

 

Y

 

最小均方误差的原理是求解 

 

其中