加载中…
个人资料
  • 博客等级:
  • 博客积分:
  • 博客访问:
  • 关注人气:
  • 获赠金笔:0支
  • 赠出金笔:0支
  • 荣誉徽章:
正文 字体大小:

zernike的多项式计算方法

(2009-03-17 21:48:16)
标签:

杂谈

分类: VC++图像处理
zernike的多项式计算方法

直接计算法:

就是把Rnmr通过定义计算

递归计算法:

Rnmr转换成计算Bnm,再对Bnm累积求和得到Rnmr

递归:

Rnmr = B(n,m,m)*r^m   + B(n,m,m+2) *r^(m+2) +... B(n,m,n)*r^n

n>m,and n-m is even

B(n,n,n) = 1

so,   B(n,n-2,n)= B(n,n,n) * (n+m)/(n-m+2)

B(n,m,n) = B(n,n,n) * ((n+m)/(n-m+2)) ^ (n-m)/2   , eg, (n-(n-2) )/2 =1

 

B(n,m,n-2) = (-1)*B(n,m,n)* (n+m)*(n-m) / (n+n)*(n-n+2)

B(n,m,n-4) = (-1)*B(n,m,n-2)* (n-2+m)(n-2-m) / (n-2+n)(n-(n-2)+2)

 

//递归计算法

double getRnmr(int n,int m,double r)
{
   
    double Rnmr = 0;
double Bnnn=1;
int times   = (int)(n-m)/2;
    int numbers = times+1;

    double mother = (double)(n+m)/(n-m+2);
     double Bnmn =   pow((double)mother,(double)times);

   double* Bnmk = new double[numbers];
    
     Bnmk[0] = Bnmn;

      int s=(int)n;

//Bnmk[0] is Bnmn , Bnmk[1] is Bnm(n-2) ... Bnmk[numbers-1] is Bnmm respectively

//Rnmr = Bnmk[0]*r^n + Bnmk[1]*r^(n-2)+...+Bnmk[numbers-1]*r^(m)

   for(int i=0;i<(numbers-1);i++)
   {
    double coeff = ((s+m)*(s-m)) / ((s+n)*(n-s+2));
    Bnmk[i+1]   = (-1)*(Bnmk[i])*coeff;
    s=s-2;
   }
   int temp = numbers-1;
     
   //Compute Rnmr
   for(int k=(int)m;k<=n;k=k+2)
   {
     Rnmr = Rnmr + (Bnmk[temp])*(pow((double)r,(double)k));

          temp--;
   }
    

   delete[] Bnmk;
         

   return Rnmr;

}

 

//直接计算法

double factorial(long n)
{
if(n < 0)
   return(0.0) ;
if(n == 0)
   return(1.0) ;
else
   return(n * factorial(n-1)) ;
}

double zernikeR(int n, int l, double r)
{
int m ;
double sum = 0.0 ;

if( ((n-l) % 2) != 0 )
    {
      cout <<"zernikeR(): improper values of n,l\n" ;
      return(0) ;
    }

for(m = 0; m <= (n-l)/2; m++)
    {
      sum += (pow((double)-1.0,(double)m)) * ( factorial(n-m) ) /
( factorial(m) * (factorial((n - 2*m + l) / 2)) *
   (factorial((n - 2*m - l) / 2)) ) *
( pow((double)r, (double)(n - 2*m)) );
    }

return(sum) ;

0

阅读 收藏 喜欢 打印举报/Report
  

新浪BLOG意见反馈留言板 欢迎批评指正

新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 产品答疑

新浪公司 版权所有