求Zernike多项式的偏导数

目的：根据Zernike多项式的偏导数公式撰写代码

由于前一节介绍了使用Matlab代码实现Zernike多项式，而本次的目的是用C语言实现对Zernike多项式求导，得到数值解（Matlab版暂时还没有整理好）。在使用夏克-哈特曼波前探测从波前斜率中重构出原始波前需要用到Zernike多项式的偏导数，该方法也称为Zernike模式法重建【1，2】，与之对应的是区域法重构波前【3】。

本文中的代码部分包括Zernike多项式，Zernike多项式在 x 和y 方向上的偏导数。需要修改的参数：宏定义的Zernike系数的个数(NumofZerCoe) 和 二维数组的分辨率(或称为像素数：Pixels)。

1 理论部分介绍

1.1 Zernike多项式公式【4】

其中，

1.2 Zernike多项式的偏导数【5】

有了Zernike多项式偏导数的表示式就可以用代码实现它了。

重要说明：Zernike多项式定义在单位圆（非离散下）是正交的，但是用代码实现它的时候设置的网格是离散的，严格来说Zernike多项式并不是正交的。那么，Zernike多项式的偏导数呢，相关文献也给出了结论：即Zernike多项式的偏导数之间并不是正交的，本人能力有限给不了证明哦，看过的文献也不记得了哈哈，等以后注意到了再回来修改。

2. Matlab代码实现

function [Zernike, DZernike_Dx, DZernike_Dy] = DzernikeDxDy(j, res)

    x           = linspace(-1,1,res);
    [x,y]       = meshgrid(x,x);
    [theta,rho] = cart2pol(x,y);% 由(x,y)换算(r,theta)
    [n,m]       = Noll_to_Zernike(j);
    if m == 0
       deltam = 1;
    else
       deltam = 0;
    end
    Norm        = sqrt(2*(n+1)/(1+deltam));% 归一化因子
    Rnm_rho     = zeros(res, res);         % 径向多项式初始化
    DRnm_Drho   = zeros(res, res);         % 径向多项式对极径偏导数的初始化
    for s = 0:(n-abs(m))/2
        Rnm_rho =   Rnm_rho + (-1)^s.* ...
                    prod(1:(n-s))* ...
                    rho.^(n-2*s)/ ...
                    (prod(1:s)*...
                    prod(1:((n+abs(m))/2-s))* ...
                    prod(1:((n-abs(m))/2-s))); % 径向多项式
  
 
        DRnm_Drho = DRnm_Drho + (n-2*s)*(-1)^s.*...
                    prod(1:(n-s))*...
                    rho.^(n-2*s-1)/...
                    (prod(1:s)*...
                    prod(1:((n+abs(m))/2-s))*...
                    prod(1:((n-abs(m))/2-s))); % 径向多项式 对极径的偏导数
    end
    
    DT_Dtheta     = zeros(res, res);
    if m < 0
        DT_Dtheta =  abs(m) * cos( abs(m) * theta );
    else 
        DT_Dtheta = -abs(m) * sin( abs(m) * theta ) ;
    end
            
    Zernike        = zeros(res, res); % Zernike多项式初始化
    DZernike_Dx    = zeros(res, res); % Zernike多项式在x方向上的初始化
    DZernike_Dy    = zeros(res, res); % Zernike多项式在y方向上的初始化
    if m < 0
        T           = sin(abs(m) .* theta);
        Zernike     = -Norm.*Rnm_rho .* sin(m .* theta); % m<0时候的zernike多项式
        DT_Dtheta   = abs(m) * cos( abs(m) * theta );
        
        DZernike_Dx = Norm .* ( DRnm_Drho .* T .* cos(theta) - ...
                     (Rnm_rho ./rho) .* DT_Dtheta .* sin(theta));
                 
        DZernike_Dy = Norm .* ( DRnm_Drho .* T .* sin(theta) + ...
                     (Rnm_rho ./rho) .* DT_Dtheta .* cos(theta));       
    else
        T           = cos( abs(m) .* theta);
        Zernike     = Norm.*Rnm_rho .*cos( m .* theta);  % m>0时候的zernike多项式
        DT_Dtheta   = -abs(m) * sin( abs(m) * theta ) ;
        
        DZernike_Dx = Norm .* ( DRnm_Drho .* T .* cos(theta) - ...
                     (Rnm_rho ./rho) .* DT_Dtheta .* sin(theta));
                 
        DZernike_Dy = Norm .* ( DRnm_Drho .* T .* sin(theta) + ...
                     (Rnm_rho ./rho) .* DT_Dtheta .* cos(theta));        
    end
    Zernike     = Zernike .* (rho<=1);
    DZernike_Dx = DZernike_Dx .* (rho<=1);
    DZernike_Dy = DZernike_Dy .* (rho<=1);
end

 调用方式如下：

clc; clear all; close all
num         = 54;                    % Zernike模式的个数
nPx         = 60；                   % 分辨率     
Zerpol      = zeros(nPx, nPx, num);
DZernike_Dx = zeros(nPx, nPx, num);
DZernike_Dy = zeros(nPx, nPx, num);
for i= 1:num
    [Zerpol(:,:,i), DZernike_Dx(:,:,i), DZernike_Dy(:,:,i)] = DzernikeDxDy(i, nPx);
end

function [n,m] = Noll_to_Zernike(j)
    if(j<1)
        error('Noll indices start at 1.');
    end

    n  = 0;
    m  = 0;
    j1 = j-1;
    while (j1 > n)
        n  = n + 1;
        j1 = j1 - n;
        m  = (-1)^j * (mod(n,2) + 2*floor((j1+mod(n+1,2))/2));
    end
end

参考文献：

[1] WANG J Y, MARKEY J K. Modal compensation of atmospheric turbulence phase distortion [J]. JOSA, 1978, 68(1): 78-87.

[2] CUBALCHINI R. Modal wave-front estimation from phase derivative measurements [J]. JOSA, 1979, 69(7): 972-7.

[3] SOUTHWELL W H. Wave-front estimation from wave-front slope measurements [J]. JOSA, 1980, 70(8): 998-1006.

[4] NOLL R J. Zernike polynomials and atmospheric turbulence [J]. JOsA, 1976, 66(3): 207-11.

[5] Efficient and stable numerical method for evaluation of Zernike polynomials and their Cartesian derivatives

3.C代码如下

# include <stdio.h>
# include <math.h>
# include <malloc.h>
#define	NumofZerCoe	30	// Zernike系数的个数，从piston开始
#define	Pixels		60	// 像素数
#define	NR_END		0
#define	NRANSI
#define	FREE_ARG	char*

typedef struct XYRT
{/*  
 在结构体中 定义直角坐标(X,Y)和极坐标(R,T)，然后用结构体进行传参、节省函数的形式参数的个数
 */
	double X[Pixels][Pixels];
	double Y[Pixels][Pixels];
	double R[Pixels][Pixels];
	double T[Pixels][Pixels];
}st_XT, *st_pXT;
/* * * * * * * * * * * * * * * *
 * 函数的声明
 * * * * * * * * * * * * * * * */
int		power_int(int a, int n);
double		power_double(double a, int n);
double		factoria_double(int value);
double		DZernikeDtheta(int n, int m, double Radial, double Theta);
double		DZernikeDradial(int n, int m, double Radial, double Theta);
double		NormCoeficients(int n, int m);
double		DRadialPolynomialsDr(int n, int m, double);
double		ZernikeRadialPolynomials(int n, int m, double);
double		**matrix(long nrl, long nrh, long ncl, long nch);
void		Noll_to_Zernike(int j, int *n, int *m);
void		ZernikePolynomials(int n, int m, st_pXT CartesiantoandPolar, double **Zernikepoly);
void		DZernikeDxDyFromCartesian(int n, int m, st_pXT CartesiantoandPolar, double ** DzernikeDx, double ** DzernikeDy);
void		free_matrix(double **m, long nrl, long nrh, long ncl, long nch);
void		ZernikeMask(st_pXT CartesiantoandPolar, double **mask);
void		Transpose(double **X, double **X_T);
void		meshgrid(double *x, st_pXT Coor);
void		linspace(double x1, double x2, int n, double *y);
void		InitMatrix(double **Ma, int Long, int Width);

/* * * * * * * * * * * * * * * *
 * 阶乘的运算，result = value ! 
 * 输入：int 类型的 value
 * 输出：阶乘后的结果- double类型
 * * * * * * * * * * * * * * * */
double factoria_double(int value) {
	double result = 1.0;
	while (value != 0) {
		result = result * value;
		value--;
	}
	return result;
}

/* * * * * * * * * * * * * * * *
 * 幂指数的运算，result = a^n
 * 输入：double 类型的 a, 
 * 输出：double类型
 * * * * * * * * * * * * * * * */
double power_double(double a, int n)
{
	double result = 1.0;
	for (int i = 0; i < n; i++) {
		result = result * a;
	}
	return result;
}

/* * * * * * * * * * * * * * * *
 * 幂指数的运算，result = a^n
 * 输入：int 类型的 a,
 * 输出：int类型
 * * * * * * * * * * * * * * * * */
int power_int(int a, int n){	
	double result = 1.0;
	for (int i = 0; i < n; i++) {
		result = result * a;
	}
	return result;
}

/* * * * * * * * * * * * * * * *
 * 函数目的：根据Zernike多项式 的序号j 求 径向级数(n) 和 角向级数(m)
 * 输入参数：序号j， 指针n 和 m，
 * 无返回值；可以直接对指针操作，把求得的结果放在 n 和 m 中
 * * * * * * * * * * * * * * * */
void Noll_to_Zernike(int j, int *n, int *m){
	if (j < 1)
		printf(" Noll indices j must be more than 1 !");
	else{
		int nn	= 0;
		int mm	= 0;
		int j1	= j - 1;
		int ss	= 0;
		while (j1 > nn){
			nn		= nn + 1;
			j1		= j1 - nn;
			double s	= (j1 + ((nn + 1) % 2)) / 2.0;
			ss		= floor(s);
			mm		= power_int(-1, j) * ((nn % 2) + 2 * ss);
		}
		*m	= mm;
		*n	= nn;
	}
}

/* * * * * * * * * * * * * * * * *
 * 函数目的：和Matlba 中 linspace 用法一样
 * 输入参数：x1 :起始点； y1:终止点； n：采样点数； *y：返回的结果
 * 输出：可以直接对指针操作，就不设置了返回值，
 * * * * * * * * * * * * * * * * */
void linspace(double x1, double x2, int n, double *y){
	double d	= (x2 - x1) / (n - 1);
	for (int i	= 0; i < n; i++){
		*(y+i)	= x1 + i * d;
	}
}

/* * * * * * * * * * * * * * * * *
 * 函数目的：和Matlba 中 meshgrid 用法一样
 * 输入参数：*x
 * 输出：可以直接对指针操作，就不设置了返回值，
 * * * * * * * * * * * * * * * * */
void meshgrid(double *x, st_pXT CartesiantoandPolar){
	for (int i = 0; i < Pixels; i++) {
		for (int j = 0; j < Pixels; j++) {
			CartesiantoandPolar->X[i][j] = *(x + j);
			CartesiantoandPolar->Y[j][i] = CartesiantoandPolar->X[i][j];
		}
	}
}

/* * * * * * * * * * * * * * * * *
 * 函数目的：二维数组的转置
 * 输入参数：转置前的参数：**X，转置后的参数： **X_T 
 * 输出：可以直接对指针操作，就不设置了返回值
 * * * * * * * * * * * * * * * * */
void Transpose(double **X, double **X_T){	
	int i, j;
	double temp;
	for (i = 0; i < Pixels; i++) {
		for (j = 0; j < Pixels; j++) {
			X_T[j][i] = X[i][j];
		}
	}
}

/* * * * * * * * * * * * * * * * *
 * 函数目的：求出单位圆的 pupil mask
 * 输入参数：结构体中 直角坐标和极坐标(-1,1范围内)；用指针的指针mask接收最后的结果
 * 输出：可以直接对指针操作，就不设置了返回值
 * 直角坐标转为极坐标,，也要得到相应的mask，mask是一个圆，其内部元素只有0和1
 * * * * * * * * * * * * * * * * */
void ZernikeMask(  st_pXT CartesiantoandPolar, double **mask ){
	for (int i = 0; i < Pixels; i++){
		for (int j = 0; j < Pixels; j++){
			CartesiantoandPolar->R[i][j] = sqrt((CartesiantoandPolar->X[i][j]) * (CartesiantoandPolar->X[i][j]) + (CartesiantoandPolar->Y[i][j]) * (CartesiantoandPolar->Y[i][j])  );// 极径
			CartesiantoandPolar->T[i][j] = atan2( CartesiantoandPolar->X[i][j] ,  CartesiantoandPolar->Y[i][j]  );		// 极角
			if (CartesiantoandPolar->R[i][j] <= 1.0){
				mask[i][j] = 1.0;	
			}else
				mask[i][j] = 0;
		}
	}
}

/* * * * * * * * * * * * * * * * *
 * 函数目的：申请动态内存后 用于 初始化 二级指针函数（指针的指针），这一步必不可少
 * 输入参数：待初始化的 变量Ma(指针的指针)，二维数组的长度：Long，二维数组的长度：Width
 * 输出：可以直接对指针操作，就不设置了返回值
 * * * * * * * * * * * * * * * * */
void InitMatrix(double **Ma, int Long, int Width){
	for (int i = 0; i < Long; i++) {
		for (int j = 0; j < Width; j++) {
			Ma[i][j] = 0.0;	
		}
	}
}

/* * * * * * * * * * * * * * * * *
 * 函数目的：用指针的指针后 需要申请动态内存，不然，用二维数组作为函数的形式参数时必须要设置数组的某个维度
 * 输入参数：对于 m*n的二维数组的													行方向索引的起点 nrl(尽量用0)，行方向索引的终点  nrh( m - 1)，													列方向索引的起点 ncl(尽量用0)，列方向索引的终点 nch( n - 1)，
 * 输出：申请后的结果。
 * 用法： 为二维数组A[nrl, nrh][ncl, nch]  申请内存，
 * 一般来说，索引都是从 0 开始，此时的 nrl = ncl = 0;
 * 数组的长 和 宽 分别为 m 和 n，则		:
	nrh  = nrl + n -1 = n-1
	nch  = ncl + m - 1 = m - 1
  * * * * * * * * * * * * * * * * */
double **matrix(long nrl, long nrh, long ncl, long nch){
	long i, nrow = nrh - nrl + 1, ncol = nch - ncl + 1;
	double **m;
	m = (double **)malloc((size_t)((nrow + NR_END) * sizeof(double*)));
	if (!m){
		printf("allocation failure 1 in matrix()");
	}
	m += NR_END;
	m -= nrl;
	m[nrl] = (double *)malloc((size_t)((nrow*ncol + NR_END) * sizeof(double)));
	if (!m[nrl]){
		printf("allocation failure 2 in matrix()");
	}
	m[nrl] += NR_END;
	m[nrl] -= ncl;
	for (i = nrl + 1; i <= nrh; i++)
		m[i] = m[i - 1] + ncol;
	return m;
}

/* * * * * * * * * * * * * * * * *
 * 函数目的：申请动态内存后必须释放内存。
 * 输入参数：待释放内存的变量m(指针的指针)，对于 m*n的二维数组的																											行方向索引的起点 nrl(尽量用0)，行方向索引的终点  nrh( m - 1)																												列方向索引的起点 ncl(尽量用0)，列方向索引的终点 nch( n - 1)，
 * 输出：可以直接对指针操作，就不设置了返回值
  * * * * * * * * * * * * * * * * */
void free_matrix(double **m, long nrl, long nrh, long ncl, long nch){
/* free a double matrix allocated by matrix() , 用法同 **matrix()   */
	free((FREE_ARG)(m[nrl] + ncl - NR_END));
	free((FREE_ARG)(m + nrl - NR_END));
}


/* * * * * * * * * * * * * * * * *
 * 函数目的：Zernike多项式的径向表达式。
 * 输入参数：径向级数(n) ，角向级数(m)，极坐标 极径(Radial)
 * 输出：    径向多项式：RadialPolynomials
 * * * * * * * * * * * * * * * * */
double ZernikeRadialPolynomials(int n, int m, double Radial) {
	double RadialPolynomials = 0.0;
	for (int s = 0; s <= (n - abs(m)) / 2; s++) {
		RadialPolynomials =	RadialPolynomials +
					power_double(-1, s) * 
					factoria_double(n - s) /
					factoria_double(s) /
					factoria_double((n +  abs(m)) / 2 - s) /
					factoria_double((n  -  abs(m)) / 2 - s) *
					power_double(Radial, n - 2 * s);
	}
	return RadialPolynomials;
}

/* * * * * * * * * * * * * * * * *
 * 函数目的：Zernike多项式的径向表达式 对 极径 的偏导数
 * 输入参数：径向级数(n) ，角向级数(m)，极坐标(Radial)，
 * 输出：    径向多项式对r 的偏导数 DrzDr
 * * * * * * * * * * * * * * * * */
double DRadialPolynomialsDr(int n, int m, double Radial){
	double DrzDr = 0.0;
	for (int s = 0; s <= (n - abs(m)) / 2; s++) {
		DrzDr =	DrzDr + 
				power_double(-1, s) * 
				factoria_double(n - s) * (n - 2 * s) /
				factoria_double(s) /
				factoria_double((n + abs(m)) / 2 - s) /
				factoria_double((n -  abs(m)) / 2 - s) *
				power_double(Radial, n - 2 * s - 1);
	}
	return DrzDr;
}

/* * * * * * * * * * * * * * * * *
 * 函数目的：Zernike多项式中的常数项
 * 输入参数：径向级数(n) ，角向级数(m)
 * 输出：    coefficients. 
 * 这个函数没有使用。
 * * * * * * * * * * * * * * * * */
double NormCoeficients(int n, int m) {
	// 求的是 归一化后的系数
	double coefficients;
	if (m != 0)
		coefficients = sqrt(2 * (n + 1));
	else
		coefficients = sqrt(2 * (n + 1) / 2);
	return coefficients;
}

/* * * * * * * * * * * * * * * * * *
 * 函数目的：Zernike多项式 对 极径r 的偏导数
 * 输入参数：径向级数(n) ，角向级数(m)，直角坐标和极角Theta
 * 输出：    Zernike多项式 对极径 的偏导数
 * * * * * * * * * * * * * * * * */
double DZernikeDradial(int n, int m, double Radial, double Theta) {
	double DrzDr	= DRadialPolynomialsDr(n, m, Radial);
	double DzDr	= 0;
	if (m > 0) {
		DzDr = DrzDr * cos( (double)m * Theta );
	}else if (m < 0) {
		DzDr = DrzDr * sin( (double)(-m)  *Theta );
	}else {
		DzDr = DrzDr;
	}
	return DzDr;
}

/* * * * * * * * * * * * * * * * * 
 * 函数目的：Zernike多项式 对 极角θ 的偏导数
 * 输入参数：径向级数(n) ，角向级数(m)，极坐标-- 极径Radial 和 极角Theta
 * 输出：    zernike多项式对极角的偏导数DzernikeDt
 * * * * * * * * * * * * * * * * */
double DZernikeDtheta(int n, int m, double Radial, double Theta){	
	double RadialPolynomials = ZernikeRadialPolynomials( n, m, Radial);
	double DzernikeDt = 0.0;
	if (m > 0) {
		DzernikeDt = -(double)m * RadialPolynomials * sin((double)m * Theta);
	}else if (m < 0) {
		DzernikeDt = -(double)m * RadialPolynomials * cos((double)(-m) * Theta);
	}else {
		DzernikeDt = 0;
	}
	return DzernikeDt;
}


/* * * * * * * * * * * * * * * * * 
 * 函数目的：Zernike多项式 对 直角坐标x 和 y 的偏导数
 * 输入参数：径向级数(n) ，角向级数(m)，直角坐标和极坐标(CartesiantoandPola--结构体指针)，对x偏导数(DzernikeDx)，对y偏导数(DzernikeDx)
 * 输出：    可以直接对指针的指针 指向的地址 操作，就不设置了返回值
 * * * * * * * * * * * * * * * * * */
void DZernikeDxDyFromCartesian(int n, int m, st_pXT CartesiantoandPolar, double ** DzernikeDx ,  double ** DzernikeDy)
{
	/*必须说明
	我在 [-1, 1] 区间上 使用 linspace() 和 meshgrid() 函数 时，把 得到的 二维网格上的点X 和Y 弄反了，
	所以下面的函数 DzernikeDx 和 DzernikeDy 需要把 公式(理论) 中的 CartesiantoandPolar->X 和 CartesiantoandPolar->Y 相互替换。
	理论中的公式 在 那篇 英文 文献中。
	否则得到的错误结果，我也是调试了一下午才发现这个错误。我也懒得改了，还请读者注意一下即可。
	*/
	double **mask = matrix(0, Pixels - 1, 0, Pixels - 1);	// mask函数
	InitMatrix(mask, Pixels, Pixels);
	ZernikeMask(CartesiantoandPolar, mask);
	for (int i = 0; i < Pixels; i++) {
		for (int j = 0; j < Pixels; j++) {
			double DzDtheta	 = DZernikeDtheta(n, m, CartesiantoandPolar->R[i][j], CartesiantoandPolar->T[i][j]);
			double DZDradius = DZernikeDradial(n, m, CartesiantoandPolar->R[i][j], CartesiantoandPolar->T[i][j]);
			DzernikeDx[i][j] = mask[i][j] * ( CartesiantoandPolar->Y[i][j] * DZDradius  / CartesiantoandPolar ->R[i][j] )   -
					   ( CartesiantoandPolar->X[i][j]  / CartesiantoandPolar->R[i][j] / CartesiantoandPolar->R[i][j] * DzDtheta);
			DzernikeDy[i][j] = mask[i][j] * ( CartesiantoandPolar->X[i][j] * DZDradius / CartesiantoandPolar->R[i][j] ) +
					  ( CartesiantoandPolar->Y[i][j]  / CartesiantoandPolar->R[i][j] / CartesiantoandPolar->R[i][j] * DzDtheta);
			if (m != 0) {
				double Coefficients	= sqrt(2 * (n + 1));
				DzernikeDx[i][j]	= DzernikeDx[i][j] * Coefficients;
				DzernikeDy[i][j]	= DzernikeDy[i][j] * Coefficients;
			}else {
				double Coefficients	= sqrt(n + 1);
				DzernikeDx[i][j]	= DzernikeDx[i][j] * Coefficients;
				DzernikeDy[i][j]	= DzernikeDy[i][j] * Coefficients;
			}
		}
	}
	free_matrix(mask, 0, Pixels - 1, 0, Pixels - 1);
}


/* * * * * * * * * * * * * * * * * * 
 * 函数目的：求Zernike多项式 
 * 输入参数：径向级数(n) ，角向级数(m)，直角坐标和极坐标(CartesiantoandPola--结构体指针)，多项式(Zernikepoly)
 * 输出：    可以直接对指针的指针 指向的地址 操作，就不设置了返回值
 * * * * * * * * * * * * * * * * * */
void ZernikePolynomials(int n , int m, st_pXT CartesiantoandPolar, double **Zernikepoly){	
	double **mask = matrix(0, Pixels - 1, 0, Pixels - 1);	// mask函数
	InitMatrix(mask, Pixels, Pixels);
	ZernikeMask(CartesiantoandPolar, mask);
	if (m < 0){
		double Coefficients = sqrt(2 * ( n + 1) );
		for (int i = 0; i < Pixels; i++) {
			for (int j = 0; j < Pixels; j++) {
				double RadialPolynomials = ZernikeRadialPolynomials(n, m, CartesiantoandPolar->R[i][j]);
				Zernikepoly[i][j] = mask[i][j] * Coefficients * RadialPolynomials * sin((double)abs(m) * CartesiantoandPolar->T[i][j]);
			}
		}
	}else if(m>0){
		double Coefficients = sqrt(2 * (n + 1));
		for (int i = 0; i < Pixels; i++){
			for (int j = 0; j < Pixels; j++){
				double RadialPolynomials = ZernikeRadialPolynomials(n, m, CartesiantoandPolar->R[i][j]);
				Zernikepoly[i][j] = mask[i][j] * Coefficients * RadialPolynomials* cos((double)m * CartesiantoandPolar->T[i][j] );
			}
		}
	}else{
		double Coefficients = sqrt( n + 1);
		for (int i = 0; i < Pixels; i++){
			for (int j = 0; j < Pixels; j++){
				double RadialPolynomials = ZernikeRadialPolynomials(n, m, CartesiantoandPolar->R[i][j]);
				Zernikepoly[i][j] = mask[i][j] * Coefficients * RadialPolynomials;
			}
		}
	}
	free_matrix(mask, 0, Pixels - 1, 0, Pixels - 1);
}

// 终于等来了主函数
int main(void){
	int n[NumofZerCoe] = { 0 }; // n初始化
	int m[NumofZerCoe] = { 0 }; // m初始化
	double x1 = -1.0, x2 = 1.0;
	double y[Pixels] = { 0 };   
	linspace(x1, x2, Pixels, y);											// ==> y = linspace(x1,x2,Pixels)
	st_XT CartesiantoandPolar = { 0 };									// 定义一个 含有 直角和极坐标的结构体 指针
	meshgrid(y, &CartesiantoandPolar);								// 实现：[XX,YY ] = meshgrid(y);
	double **mask = matrix(0, Pixels - 1, 0, Pixels - 1);// 单位圆的mask函数
	InitMatrix(mask, Pixels, Pixels);
	ZernikeMask(&CartesiantoandPolar, mask);	 // 这一步必须加上，才会得到极坐标 R 和 T， 不要单独分出来实现R 和T 

	/*FILE *fpWrite4 = fopen("CartesiantoandPolar.txt", "w");
	if (fpWrite4 == NULL) {
		return 0;
	}
	for (int i = 0; i < Pixels; i++) {
		for (int j = 0; j < Pixels; j++) {
			fprintf(fpWrite4, "%3.3f\t", mask[i][j]);
		}
		fprintf(fpWrite4, "\n");
	}
	fclose(fpWrite4);*/


	// 前 NumofZerCoe 阶模式 保存到  ZERPOL[][][]中，这是三维数组。
	double ZERPOLs[Pixels][Pixels][NumofZerCoe] = { 0 };	// 保存的是 前NumofZerCoe个 Zernike多项式
	double DZDXs[Pixels][Pixels][NumofZerCoe] = { 0 };	// 保存的是 前NumofZerCoe个 Zernike多项式对x的偏导数
	double DZDYs[Pixels][Pixels][NumofZerCoe] = { 0 };	// 保存的是 前NumofZerCoe个 Zernike多项式对y的偏导数

	for (int j = 0; j < NumofZerCoe; j++){
		double **Zernikepoly = matrix(0, Pixels - 1, 0, Pixels - 1);	// Zernike多项
		InitMatrix(Zernikepoly, Pixels, Pixels);

		double **DzernikeDx = matrix(0, Pixels - 1, 0, Pixels - 1);	// Zernike多项式对x的偏导数
		InitMatrix(DzernikeDx, Pixels, Pixels);

		double **DzernikeDy = matrix(0, Pixels - 1, 0, Pixels - 1);	// Zernike多项式对y的偏导数
		InitMatrix(DzernikeDy, Pixels, Pixels);

		Noll_to_Zernike(j + 1, &n[j], &m[j]);
		//printf("j = %d\t, n = %d\t, m = %d \n", 1 + j, n[j], m[j]);	// 显示 j  n  m  三个值

		ZernikePolynomials(n[j], m[j], &CartesiantoandPolar, Zernikepoly);

		// 调用函数  返回 多项式的 偏导数
		DZernikeDxDyFromCartesian(n[j], m[j], &CartesiantoandPolar, DzernikeDx, DzernikeDy);

		for (int ii = 0; ii < Pixels; ii++) {
			for (int jj = 0; jj < Pixels; jj++) {
				ZERPOLs[ii][jj][j] = Zernikepoly[ii][jj] * mask[ii][jj];
				DZDXs[ii][jj][j]		= DzernikeDx[ii][jj] * mask[ii][jj];
				DZDYs[ii][jj][j]		= DzernikeDy[ii][jj] * mask[ii][jj];
			}
		}
		free_matrix(Zernikepoly, 0, Pixels - 1, 0, Pixels - 1);
		free_matrix(DzernikeDx, 0, Pixels - 1, 0, Pixels - 1);
		free_matrix(DzernikeDy, 0, Pixels - 1, 0, Pixels - 1);
	}

	free_matrix(mask, 0, Pixels - 1, 0, Pixels - 1);

	FILE *fpWrite1 = fopen("Zernikes.txt", "w");
	if (fpWrite1 == NULL) {
		return 0;
	}
	for (int ii = 0; ii < NumofZerCoe; ii++) {
		for (int i = 0; i < Pixels; i++) {
			for (int j = 0; j < Pixels; j++) {
				fprintf(fpWrite1, "%3.3f\t", ZERPOLs[i][j][ii]);
			}
			fprintf(fpWrite1, "\n");
		}
	}
	fclose(fpWrite1);


	FILE *fpWrite2 = fopen("DzDxs.txt", "w");
	if (fpWrite2 == NULL) {
		return 0;
	}
	for (int ii = 0; ii < NumofZerCoe; ii++) {
		for (int i = 0; i < Pixels; i++) {
			for (int j = 0; j < Pixels; j++) {
				fprintf(fpWrite2, "%3.3f\t", DZDXs[i][j][ii]);
			}
			fprintf(fpWrite2, "\n");
		}
	}
	fclose(fpWrite2);

	FILE *fpWrite3 = fopen("DzDys.txt", "w");
	if (fpWrite3 == NULL) {
		return 0;
	}
	for (int ii = 0; ii < NumofZerCoe; ii++) {
		for (int i = 0; i < Pixels; i++) {
			for (int j = 0; j < Pixels; j++) {
				fprintf(fpWrite3, "%3.3f\t", DZDYs[i][j][ii]);
			}
			fprintf(fpWrite3, "\n");
		}
	}
	fclose(fpWrite3);
/*
matlab 中 加载 data.txt 后 ，用下面语句 画图
data = reshape(data', [60,60,20]);
figure(10);imagesc(data(:,:,5));colormap(jet);colorbar
*/
	return 0;
}

需要修改的参数：宏定义的Zernike系数的个数(NumofZerCoe) 和 二维数组的分辨率(或称为像素数：Pixels）

写入到txt中的数据 在Matlab 中加载的方式 ：

data1=importdata('data.txt');

如果写入到txt中的是三维数组的话，则

data = reshape(data', [Pixels,Pixels,NumofZerCoe]); % 转置 后 在 使用 reshape，即可读取三维数组中的数组

figure(10);imagesc(data(:,:,5));colormap(jet);colorbar

如果写入到txt中的是二维数组的话，则

data = reshape(data', [Pixels,Pixels]);

另外，常见的报错：

1）error C4996: 'fopen': This function or variable may be unsafe. Consider using fopen_s instead. To disable deprecation,

use _CRT_SECURE_NO_WARNINGS. See online help for details

解决方案，项目 =》属性 =》c/c++ =》预处理器=》点击预处理器定义，编辑，加入_CRT_SECURE_NO_WARNINGS

2）编译时没有错误，但是调试或者运行时程序就报错了！

debug调试提示：xxx.exe 中的 0x00e731d7 处未处理的异常: 0xC00000FD: Stack overflow

原来是栈溢出了！

【解决方案】扩大栈空间的大小，在vs工程里面：项目->属性->链接器->系统->堆栈保留大小（注：这里填的是字节数）如果你想把他扩大为2M的话，1024*1024*2 = 2097152(这个数字填进去)4M = 4194304