匹配 算法 总结归纳
基于像素的匹配
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 ;
}
- 序贯相似性检测法匹配 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;
}
基于特征的匹配
利用灰度信息匹配的方法主要缺陷是计算量过大，对图像的灰度变换很敏感，尤其是非线性的光照变化，此外，对目标的旋转、变形以及遮挡也比较敏感，为
了克服这些缺点，可以利用图像的特征进行匹配，由于图像的特征点比像素点要少很多，大大减少了匹配过程的计算量，同时，特征点的匹配度量值对位置的
变化比较敏感，可以大大的提高匹配的精度。而且，特征点的提取过程可以减少噪声的影响，对灰度变化、图像形变以及遮挡等有较好的适应能力。
-
不变矩匹配法 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
最小均方误差的原理是求解
其中
