像差与zernike多项式

------------------------像差------------------------

一、像差的概念：

像差是指光学系统中的成像缺陷。几何光学上把像差（几何像差）分为单色光像差和色光像差，前者包括球差、彗差、像散、场曲和畸变，后者包括位置色差和倍率色差；而物理光学上把像差称之为波前像差或波阵面像差，即是点光源发出的球面波经光学系统后形成的波形与理想球面波之间的距离。波前像差可以通过Zernike多项式周期表或球差、彗差等几何像差来表达。

1、球差是指轴上点光源发出的光线经屈光系统后，近轴光线与边缘光线像点之间的距离。存在球差的光学系统所形成的像是对称的弥散圆。

2、彗差是指轴外点光源发出的光线经屈光系统后，上光线和下光线的交点离开主光线的距离。存在彗差的光学系统形成的像是不对称的弥散斑。

3、像散是子午面上的像点和弧矢面上的像点的距离。子午面为通过光轴的平面，而弧矢面为垂直于子午面并通过主光线的平面。
4、场曲为平面物体通过光学系统后形成的矢状弯曲。人眼的视网膜是球形向后弯曲状，正好能补偿眼屈光系统产生的这种成像缺陷。
5、畸变为方形物体通过光学系统后周边各点产生了不同棱镜像移所致。

6、位置色差即轴位色差，白光中不同波长的光线经光学系统后形成像点的距离，短波长的交点近于长波长的交点。
7、倍率色差某一物体经光学系统成像后不同波长的光线在物像大小上的差异。
人眼是一复杂的光学系统，存在波前像差。波前像差分为低阶像差和高阶像差。按照Zernike多项式周期表，1-2阶为低阶像差，3阶以上为高阶像差。屈光不正属于低阶像差。

二、像差产生的原因：

人眼产生像差的原因包括各屈光面固有的成像缺陷、调节时的动态变化和各屈光面间的相互影响等多方面。
1、角膜角膜前表面不是理想的球面，而是非球面，然而中央4毫米区域却近似球形，所以产生球差。角膜顶点处较陡，边缘部较扁平。但顶点并不总在角膜的几何中心，往往偏下，不规则角膜的顶点偏离几何中心可达2mm以上。角膜各部分的厚度和曲率半径在各测量点上并不一致，这些就是角膜的不对称性和表面不规则性。
2、晶状体晶状体前表面较平坦，可抵消80%的角膜球差，但晶状体前表面也不平滑。随年龄增加，晶状体增厚，核发生硬化，各部位屈光指数不一致。晶状体的调节变化，除晶状体屈光力发生改变外，还可有X、Y、Z轴的变化。晶状体亦存在不对称性和表面不规则性。
3、其他玻璃体的变性、液化、混浊、后脱离等。泪膜的不均匀和不稳定，如干眼症或用药等影响。房水的改变。高度近视患者的视网膜形态变化等。
4、角膜和晶状体的光学中心不一致；与入瞳中心不一致。
5、光轴和视轴本身的偏差。
6、瞳孔除随光线的强弱发生改变，人群中存在相当大的生理差异。瞳孔增大，像差明显增加。入瞳中心并不在角膜的几何中心的对应点上。

--------------------zernike多项式-------------------

1. 矩的概念

图像识别的一个核心问题是图像的特征提取，即用一组简单的数据来描述整个图像，这组数据越简单越有代表性越好。良好的特征不受光线、噪点、几何形变的干扰，图像不变矩就是其中一个。
矩是概率与统计中的一个概念，是随机变量的一种数字特征。设x为随机变量，c为常数，k为正整数。则

为x关于c点的k阶矩。
1. c=0时表示x的k阶原点矩.
2. c=E(x)时表示x的k阶中心矩。
一阶原点矩就是期望。一阶中心矩μ1=0，二阶中心矩μ2就是X的方差Var(X)。在统计学上，高于4阶的矩极少使用。μ3可以去衡量分布是否有偏。μ4可以去衡量分布（密度）在均值附近的陡峭程度如何。
针对于一幅图像，我们把像素的坐标看成是一个二维随机变量(X,Y)，那么一幅灰度图像可以用二维灰度密度函数来表示，因此可以用矩来描述灰度图像的特征。
不变矩(Invariant Moments)是一处高度浓缩的图像特征，具有平移、灰度、尺度、旋转不变性。M.K.Hu在1961年首先提出了不变矩的概念。1979年M.R.Teague根据正交多项式理论提出了Zernike矩。

2. Zernike矩

Zernike矩能够很容易地构造图像的任意高阶矩，并能够使用较少的矩来重建图像。Zernike矩是基于Zernike多项式的正交化函数，虽然其计算比较复杂，但是Zernide矩在图像旋转和低噪声敏感度方面具有较大的优越性。由于Zernike矩具有图像旋转不变性，而且可以构造任意高阶矩，所以被广泛应用对目标进行识别中。

计算一副图像的Zernike矩时，必须将图像的中心移到坐标的原点，将图像的像素点映射到单位圆内，由于Zernike矩具有旋转不变性，可以将 作为图像的不变特征，其中图像的低频特征有p值小的提取，高频特征由p值高的 提取。Zernike矩可以构造任意高阶矩，由于Zernike矩只具有旋转不变性，不具有平移和尺度不变性，所以要提前对图像进行归一化，采用标准矩的方法来归一化一副图像，标准矩定义为：

 ，由标准矩我们可以得到图像的"重心"，

我们将图像的"重心"移动到单位圆的圆心（即坐标的原点），便解决了平移问题。 表征了图像的"面积"，归一图像的尺度无非就是把他们的大小变为一致的，（这里的大小指的是图像目标物的大小，不是整幅图像的大小，"面积"也是目标物的"面积"）。

所以，对图像进行变换 就可以达到图像尺寸一致的目的。

综合上面结果，对图像进行 变换，最终图像 的Zernike矩就是平移，尺寸和旋转不变的。

3.Zernike多项式

以下内容大部分摘自：http://blog.csdn.net/wrj19860202/article/details/6334275

Zernike在1934年引入了一组定义在单位圆 上的复值函数集{ }，{ }具有完备性和正交性，使得它可以表示定义在单位圆盘内的任何平方可积函数。其定义为：

 表示原点到点 的矢量长度； 表示矢量 与 轴逆时针方向的夹角， 是实值径向多项式：

称为Zernike多项式，Zernike多项式满足正交性：

 为克罗内克符号， 

 是 的共轭多项式，由于Zernike多项式的正交完备性，所以在单位圆内的任何图像 都可以唯一的用下面式子来展开：

式子中 就是Zernike矩，其定义为：

式子中 和 采用的是不同的坐标系， 采用直角坐标，而 采用的极坐标系，在计算的时候要进行坐标转换。对于离散的数字图像，可将积分形式改为累加形式：

， 

Noll index ({\displaystyle j})Radial degree ({\displaystyle n})Azimuthal degree ({\displaystyle m}){\displaystyle Z_{j}}Classical name100{\displaystyle 1}Piston211{\displaystyle 2\rho \cos \theta }Tip (lateral position) (X-Tilt)31−1{\displaystyle 2\rho \sin \theta }Tilt (lateral position) (Y-Tilt)420{\displaystyle {\sqrt {3}}(2\rho ^{2}-1)}Defocus (longitudinal position)52−2{\displaystyle {\sqrt {6}}\rho ^{2}\sin 2\theta }Oblique astigmatism622{\displaystyle {\sqrt {6}}\rho ^{2}\cos 2\theta }Vertical astigmatism73−1{\displaystyle {\sqrt {8}}(3\rho ^{3}-2\rho )\sin \theta }Vertical coma831{\displaystyle {\sqrt {8}}(3\rho ^{3}-2\rho )\cos \theta }Horizontal coma93−3{\displaystyle {\sqrt {8}}\rho ^{3}\sin 3\theta }Vertical trefoil1033{\displaystyle {\sqrt {8}}\rho ^{3}\cos 3\theta }Oblique trefoil1140{\displaystyle {\sqrt {5}}(6\rho ^{4}-6\rho ^{2}+1)}Primary spherical1242{\displaystyle {\sqrt {10}}(4\rho ^{4}-3\rho ^{2})\cos 2\theta }Vertical secondary astigmatism134−2{\displaystyle {\sqrt {10}}(4\rho ^{4}-3\rho ^{2})\sin 2\theta }Oblique secondary astigmatism1444{\displaystyle {\sqrt {10}}\rho ^{4}\cos 4\theta }Oblique quadrafoil154−4{\displaystyle {\sqrt {10}}\rho ^{4}\sin 4\theta }Vertical quadrafoil

通常使用幂级数展开式的形式来描述光学系统的像差。由于泽尼克多项式和光学检测中观测到的像差多项式的形式是一致的，因而常用zernike描述波前特性。但并不意味着泽尼克多项式就是用来拟合检测数据的最佳多项式形式。在某些情况下，用泽尼克多项式来描述波前数据具有很大的局限性。例如，当需要考虑空气扰动的时候，泽尼克多项式几乎没有什么价值。同样地，我们也无法找到一组合适的泽尼克多项式来描述单点金刚石车削加工（single point diamond turning process）中的制造误差。为了准确地描述圆锥面光学元件（conical optical elements）的对准误差，必须对泽尼克多项式进行修正。盲目地使用泽尼克多项式来表达检测数据只会导致糟糕的结果。

泽尼克多项式是由无穷数量的多项式完全集组成的，它有两个变量，ρ和θ，它在单位圆内部是连续正交的。需要注意的是，泽尼克多项式仅在单位圆的内部连续区域是正交的，通常在单位圆内部的离散的坐标上是不具备正交性质的。泽尼克多项式具有三个和其他正交多项式集不一样的性质。

1.泽尼克多项式Z(ρ, θ)可以被化解为径向坐标ρ和角度坐标θ的函数，其形式如下： Z (ρ, θ) = R ( ρ ) G ( θ )， 这里，关于角度的函数G(θ)是一个以2π弧度为周期的连续函数，并且满足当坐标系旋转α角度之后，其形式不发生改变，也就是旋转不变性。

2.泽尼克多项式的第二个性质是径向函数R ( ρ ) （Radial Function）必须是ρ的n次多项式，并且不包含幂次低于m次的ρ方项。 

3.第三个性质是当m为偶数时R(ρ)也为偶函数，m为奇数时，R(ρ)也为奇函数。 径向多项式R ( ρ )可以看作是雅可比多项式（Jacobi polynomials）的特例。

zernike矩C++程序：http://blog.csdn.net/wrj19860202/article/details/6334275

/*计算一行的像素个数 imwidth:图像宽度 ；deep:图像深度（8位灰度图为1，24位彩色图为3） */  #define  bpl(imwidth, deep) ((imwidth*deep*8+31)/32*4)  /*获取像素值；psrcBmp:图像数据指针；nsrcBmpWidth:图像宽度，以像素为单位x,y:像素点 deep:图像的位数深度，（1表示8位的灰度图，3表示24位的RGB位图） */  COLORREF J_getpixel( const BYTE *psrcBmp, const int nsrcBmpWidth, const int x, const int y, int deep = 3)  {      if (deep == 3)      {          return RGB(*(psrcBmp + x*3 + y*bpl(nsrcBmpWidth, deep) + 2 ) ,               *(psrcBmp + x*3 + y*bpl(nsrcBmpWidth, deep) + 1 ) ,               *(psrcBmp + x*3 + y*bpl(nsrcBmpWidth, deep) +0 ));      }      else if (deep == 1)      {          return *(psrcBmp + x + y*bpl(nsrcBmpWidth, deep));      }  }    //获取标准矩(只支持8位灰度图)  void GetStdMoment(BYTE *psrcBmp ,                    int nsrcBmpWidth,                   int nsrcBmpHeight,                   double *m)  {      for ( int p = 0 ; p < 2 ; p++ )          for ( int q = 0 ; q < 2 ; q++ )          {              if( p == 1 && q == 1)                  break;              for ( int y = 0 ; y < nsrcBmpHeight ; y++ )                  for ( int x = 0 ; x < nsrcBmpWidth ; x++ )                      m[p*2+q] += (pow( (double)x , p ) * pow( (double)y , q ) * J_getpixel(psrcBmp , nsrcBmpWidth , x ,y, 1));          }  }    //阶乘  double Factorial( int n )  {      if( n < 0 )          return -1;        double m = 1;      for(int i = 2 ; i <= n ; i++)      {          m *= i;      }      return m;  }    //阶乘数，计算好方便用，提高速度  double factorials[11] = {1 , 1 , 2 , 6 , 24 , 120 , 720 , 5040 , 40320 , 362880 , 39916800};    //把图像映射到单位圆，获取像素极坐标半径  double GetRadii(int nsrcBmpWidth,         int nsrcBmpHeight,         int x0,         int y0,         int x,         int y)  {      double lefttop = sqrt(((double)0 - x0)*(0 - x0) + (0 - y0)*(0 - y0));      double righttop = sqrt(((double)nsrcBmpWidth - 1 - x0)*(nsrcBmpWidth - 1 - x0) + (0 - y0)*(0 - y0));      double leftbottom = sqrt(((double)0 - x0)*(0 - x0) + (nsrcBmpHeight - 1 - y0)*(nsrcBmpHeight - 1 - y0));      double rightbottom = sqrt(((double)nsrcBmpWidth - 1 - x0)*(nsrcBmpWidth - 1 - x0) + (nsrcBmpHeight - 1 - y0)*(nsrcBmpHeight - 1 - y0));        double maxRadii = lefttop;      maxRadii < righttop ? righttop : maxRadii;      maxRadii < leftbottom ? leftbottom : maxRadii;      maxRadii < rightbottom ? rightbottom : maxRadii;        double Radii = sqrt(((double)x - x0)*(x - x0) + (y - y0)*(y - y0))/maxRadii;      if(Radii > 1)      {          Radii = 1;      }      return Radii;  }    //把图像映射到单位圆，获取像素极坐标角度  double GetAngle(int nsrcBmpWidth,                  int nsrcBmpHeight,                  int x,                  int y)  {      double o;        double dia = sqrt((double)nsrcBmpWidth*nsrcBmpWidth + nsrcBmpHeight*nsrcBmpHeight);      int x0 = nsrcBmpWidth / 2;      int y0 = nsrcBmpHeight / 2;      double x_unity = (x - x0)/(dia/2);       double y_unity = (y - y0)/(dia/2);        if( x_unity == 0 && y_unity >= 0 )          o=pi/2;      else if( x_unity ==0 && y_unity <0)          o=1.5*pi;      else          o=atan( y_unity / x_unity );      if(o*y<0)    //第三象限          o=o+pi;        return o;  }    //Zernike不变矩  J_GetZernikeMoment(BYTE *psrcBmp ,                   int nsrcBmpWidth,                  int nsrcBmpHeight,                  double *Ze )  {      double R[count][count] = {0.0};      double V[count][count] = {0.0};        double M[4] = {0.0};      GetStdMoment(psrcBmp , nsrcBmpWidth , nsrcBmpHeight , M);      int x0 = (int)(M[2]/M[0]+0.5);      int y0 = (int)(M[1]/M[0]+0.5);        for(int n = 0 ; n < count ; n++)      {          for (int m = 0 ; m < count ; m++)          {              //优化算法，只计算以下介数                if( (n == 1 && m == 0) ||                  (n == 1 && m == 1) ||                  (n == 2 && m == 0) ||                  (n == 2 && m == 1) ||                  (n == 2 && m == 2) ||                  (n == 3 && m == 0) ||                  (n == 3 && m == 1) ||                  (n == 3 && m == 2) ||                  (n == 3 && m == 3) ||                  (n == 4 && m == 0) ||                  (n == 4 && m == 1) ||                  (n == 4 && m == 2) ||                  (n == 4 && m == 3) ||                  (n == 4 && m == 4))                {                  for(int y = 0 ; y < nsrcBmpHeight ; y++)                  {                      for (int x = 0 ; x < nsrcBmpWidth ; x++)                      {                          for(int s = 0 ; (s <= (n - m)/2 ) && n >= m ; s++)                          {                              R[n][m] += pow( -1.0, s )                                      * ( n - s > 10 ? Factorial( n - s ) : factorials[ n - s ] )                                      * pow( GetRadii( nsrcBmpWidth, nsrcBmpHeight, x0, y0, x, y ), n - 2 * s )                                      / ( ( s > 10 ? Factorial( s ) : factorials[ s ] )                                      * ( ( n + m ) / 2 - s > 10 ? Factorial( ( n + m ) / 2 - s ) : factorials[ ( n + m ) / 2 - s ] )                                      * ( ( n - m ) / 2 - s > 10 ? Factorial( ( n - m ) / 2 - s ) : factorials[ ( n - m ) / 2 - s ] ) );                          }                          Ze[ n * count + m ] += R[ n ][ m ]                                              * J_getpixel( psrcBmp, nsrcBmpWidth, x ,y, 1)                                              * cos( m * GetAngle( nsrcBmpWidth, nsrcBmpHeight, x, y) );//实部                            V[n][m] += R[ n ][ m ]                                   * J_getpixel( psrcBmp, nsrcBmpWidth, x, y, 1)                                  * sin( m * GetAngle( nsrcBmpWidth, nsrcBmpHeight, x, y ) );//虚部                            R[n][m] = 0.0;                      }                  }                  *(Ze+n*count + m) = sqrt( (*(Ze+n*count + m))*(*(Ze+n*count + m)) + V[n][m]*V[n][m] )*(n+1)/pi/M[0];              }          }      }  }  

参考：

http://changzhengdr.blog.sohu.com/108180288.html

http://download.csdn.net/detail/piaoxuezhong/9791297

https://wenku.baidu.com/view/4576361614791711cc791761.html

http://blog.csdn.net/wrj19860202/article/details/6334275