模板匹配之zernike矩

zernike矩具有旋转不变性，广泛用于目标识别模板匹配，抗噪能力强。

zernike矩的原理就不过多赘述，上代码。

   
   //求n的阶乘
   long TemplateMatch::factorial(int n)
   {
   	int result=1;
   	while(n>1)
   	{
   		result *= n;
   		n-=1;
   	}
   	return result;
   }
   
   
   //zernike矩径向多项式
   double TemplateMatch::RadialPolynomial(double p, int n, int m)
   {
   	double radial=0;
   	int s;
   	double c;
   	
   	for (s=0; s<=(n-m)/2; s++)
   	{
   		c=pow((double)(-1),(double)(s))*factorial(n-s)/((double)(factorial(s)*factorial((n+m)/2-s)*factorial((n-m)/2-s)));
   		radial=radial+c*pow(p,n-2*s);
   	}
   	return radial;
   }
   
   
   //n阶m次zernike矩
   double TemplateMatch::ZernikeMoments(double *im,int imW,int imH, int n, int m)
   {
   	#define im(ROW,COL) im[imW*(ROW)+(COL)]
   	double zr=0;
   	double zi=0;
   	int cnt=0;
   	int x,y;
   	double p;
   	double radial=0;
   	double theta=0;
   
   	for (y=0; y<imH; y++)
   	{
   		for(x=0; x<imW; x++)
   		{
   			p=sqrt((double)((2*x-imW+1)*(2*x-imW+1)+(imH-1-2*y)*(imH-1-2*y)))/imW;
   			if (p<=1)
   			{
   				radial=RadialPolynomial(p,n,m);
   				theta=atan2((double)(imH-1-2*y),(double)(2*x-imW+1));
   				zr+=im(y,x)*radial*cos(m*theta);
   				zi+=im(y,x)*radial*sin(m*theta);
   				cnt++;
   			}
   		}
   	}
   	double A=(n+1)*sqrt((double)(zr*zr+zi*zi))/cnt;
   	return A;
   }