泽尼克波前拟合 matlab-编程知识

tic;                                             %tic()函数可以用来启动一个计时器m
zernikePhase = unwrappedPhaseMask((img_rows/2-0.3*img_rows+1):(img_rows/2+0.3*img_rows),(img_cols/2-0.3*img_cols+1):(img_cols/2+0.3*img_cols));
[zernike_rows, zernike_cols] = size(zernikePhase);
newmask = zeros(zernike_rows, zernike_cols);
newmask(isnan(zernikePhase))=1;% 2、Fringe zernike序列
indices = [0 1  1 2 2  2 3  3 4 3  3 4  4 5  5 6 4  4 5  5 6  6 7  7 8 5  5 6  6 7  7 8  8 9  9 10 12;...0 1 -1 0 2 -2 1 -1 0 3 -3 2 -2 1 -1 0 4 -4 3 -3 2 -2 1 -1 0 5 -5 4 -4 3 -3 2 -2 1 -1 0 0];
indices = indices.';zernikeDPI = size(zernikePhase);%图像分辨率
zernikeRows = zernikeDPI(1);%行数
zernikeCols = zernikeDPI(2);%列数zernikeMats = zernike_mats(zernikePhase, indices, newmask);    %计算相位数据在这些 Zernike 多项式上的投影，即计算 Zernike 模式。
a_coeffs = zernike_moments(zernikePhase, indices, newmask);         % zernike36项系数 a_coeffs  函数来计算相位数据的 Zernike 系数。zernike1to3 = zeros(zernikeRows,zernikeCols);   %计算并从相位数据中去除前三项（通常是倾斜和偏焦）的影响
zernike4to36 = zeros(zernikeRows,zernikeCols);
% 1-3项
for i =1:3zernike1to3 = zernike1to3 + a_coeffs(i).*zernikeMats(:,:,i);  % zernike前三项拟合的平面
end
% 4-36项 ---这个暂时没用上
for i =4:37zernike4to36 = zernike4to36 + a_coeffs(i).*zernikeMats(:,:,i);  % zernike 4-36项拟合的平面
endzernikeResult = zernikePhase - zernike1to3;

启动计时器：使用tic函数开始记录一个操作过程的时间。
提取相位图的中心区域：从一个更大的相位图像中根据给定的行数img_rows和列数img_cols，裁剪出中心的一个区域，裁剪比例为宽和高的60%。
创建新的掩码：基于裁剪后的相位图，创建一个新的掩码newmask，标记出相位图中的NaN值，因为这些值在计算中可能会引起问题。
定义Zernike多项式的指标：这一系列的指标用于计算与相位图相关的Zernike多项式，它们对应于不同的失真类型。
初始化Zernike矩阵和系数：通过调用zernike_mats函数构建一个与Zernike多项式相对应的矩阵数组，并通过zernike_moments计算相位数据的Zernike系数，这些系数用于描述波前的失真。
分离和去除前三项影响：前三项Zernike多项式（通常对应于倾斜和偏焦）通常会从相位数据中去除，因为它们表示整个波前的全局移动和变形，而不是更细节的特定像差。这是通过累加与前三个系数相乘的Zernike矩阵实现的。
计算剩余的Zernike多项式影响：尽管在这段代码中没有被用到，计算第4到第36项多项式的平面由一个循环完成，这可以用于分析更高级的波前失真。
结果：最终结果zernikeResult是原始相位图减去前三个Zernike多项式的影响，它代表去除全局移动和变形后的波前误差。
在这段代码末尾，很可能调用toc函数来停止计时器并打印出整个处理过程所需的时间。

这个过程通常在光学测量和仪器校准中使用。注意，此代码示例依赖于外部定义的各种函数，如unwrappedPhaseMask（用于处理或裁剪相位图），zernike_mats，zernike_moments，和zernike（它们用于计算Zernike多项式和其矩阵以及系数）。所有这些函数都需要事先正确实现，才能使此代码正常工作。

zernike_mats

function zernikeMatrices = zernike_mats(im,zernike_indices,mask)% Returns an array of matrices corresponding to the Zernike polynomials% of degrees n and m as defined by the two columns of zernike_indices,% where the first column contains values of n and the second column% contains values of m. The dimensions of these matrices will be the% same as that of im.% This code is intended to be used as a helper function for% zernike_moments and zernike_recreation.%% Example:%     %     % Create some function z(x,y) over a 125x75 grid.%%     x=linspace(-1,1,75);%     y=linspace(1,-1,125);%     [x,y] = meshgrid(x,y);%     z = x.^2+y;%     figure,imagesc(elliptical_crop(z,1));%     colormap jet;%     title('z(x,y)');%%     % Specify the indices of Zernike polynomials of interest. In this%     % case, the indices correspond to the "piston", "x-tilt", "y-tilt",%     % and "power" Zernike polynomials.%%     indices = [0 0; 1 -1; 1 1; 2 0];%     zern_mats = zernike_mats(z,indices);%     for i = 1:size(zern_mats,3)%         figure, imagesc(elliptical_crop(zern_mats(:,:,i),1));%         colormap jet;%         title(indices(i,:));%     end%       %     % The images displayed show these four Zernike polynomials over%     % 125x75 grids.%% Functions required for use: zernike, zernike_radial, elliptical_crop % See also: zernike_moments, zernike_recreation% % Evan Czako, 8.14.2019% -------------------------------------------dim=size(im);x=linspace(-1,1,dim(2));y=linspace(1,-1,dim(1));[x,y] = meshgrid(x,y);[t,r] = cart2pol(x,y);zernikeMatrices = [];for i = 1:size(zernike_indices,1)zernikeMatrices(:,:,i) = zernike(r,t,zernike_indices(i,1),zernike_indices(i,2));zernikeMatrices(:,:,i) = elliptical_crop(zernikeMatrices(:,:,i),1,mask);endend

这段MATLAB函数zernike_mats的代码目的是创建一个矩阵阵列，这些矩阵对应于Zernike多项式，这些多项式通过zernike_indices索引数组中的n和m的值来定义。输入图像im的尺寸将决定这些Zernike矩阵的尺寸。代码首先创建一个基于输入图像尺寸的笛卡尔坐标网格，然后将这个网格转换为极坐标系，之后，根据zernike_indices给出的每对n和m值，调用zernike函数计算相应的Zernike多项式，并将结果存储到zernikeMatrices数组中。最后，利用elliptical_crop函数将每个Zernike多项式矩阵裁剪到椭圆形的遮罩。

当然，让我来逐行解释这个函数的每个步骤：

1. `dim=size(im);`：这行代码计算输入图像`im`的尺寸。`dim(1)`表示图像的行数，`dim(2)`表示图像的列数。

2. `x=linspace(-1,1,dim(2));`：创建一个在-1到1之间的等距向量`x`，元素数量与图像的列数相同。这用于创建后续的网格坐标。

3. `y=linspace(1,-1,dim(1));`：创建一个在1到-1之间的等距向量`y`，元素数量与图像的行数相同。注意`y`是从1开始减到-1，与`x`相反。

创建范围在-1到1之间的原因是Zernike多项式通常在单位圆上定义，也就是说，它们是在圆内部进行计算的，圆的半径为1。单位圆内部的点有一个很重要的特性：它们的半径（距离原点的距离）不会超过1。因此，通过使用-1到1的范围来创建x和y坐标，我们可以确保所有生成的点都位于或接近单位圆的内部。这样，通过笛卡尔到极坐标的转换，我们得到的点都是有效的Zernike多项式的输入参数。使用这个范围，我们可以覆盖单位圆内的每一个点，包括圆心（0,0），圆的边界（半径为1的点），以及圆内的所有其他点。Zernike多项式是通过半径（r）和角度（θ）来定义的，所以通过适当地选择坐标范围，我们可以确保Zernike多项式能够在单位圆上正确计算，以此来表示波前的形状或光学系统的像差。在实际应用中，Zernike多项式被广泛用于分析和校正光学系统中的像差，如望远镜、显微镜或眼科检查设备。单位圆的范围选择就成为了一个自然的约定，因为它提供了一个标准化的方法来表示和计算这些多项式。

4. `[x,y] = meshgrid(x,y);`：使用`x`向量和`y`向量生成一个二维的网格。这个网格中的每一个点(x,y)代表输入图像上对应位置的坐标。

meshgrid这个函数在MATLAB中是用来从两个一维数组生成一个二维坐标网格的。当你将两个向量x_1和y_1传递给meshgrid时，它会返回两个矩阵：X和Y。

以下是[X, Y] = meshgrid(x_1, y_1);如何操作的逐步解释：

输入变量x_1和y_1：

x_1：一个向量，包含网格上所有点的x坐标。
y_1：一个向量，包含网格上所有点的y坐标。
通常，x_1和y_1都是从一个最小值开始，以一定的间隔延伸到一个最大值。

输出[X, Y]：

X：一个矩阵，每一行都是x_1的复制，并且行数和y_1中元素的数量一样多。
Y：一个矩阵，每一列都是y_1的复制，并且列数和x_1中元素的数量一样多。

操作过程：

meshgrid函数会取x_1向量，垂直复制，创建多行，从而生成矩阵X。
然后它会取y_1向量，水平复制，创建多列，生成矩阵Y。

结果网格：

结果矩阵X和Y创建了一个坐标系统，X包含了网格上每个点的水平坐标（x坐标），Y则包含每个点的垂直坐标（y坐标）。这样你就可以在由x_1和y_1定义的范围内对函数进行网格评估。
% 定义向量 x_1 = [1 2 3];
% x坐标 y_1 = [4 5];
% y坐标
% 创建网格
[X, Y] = meshgrid(x_1, y_1);
% 得到的矩阵：
%X为：
1 2 3
1 2 3
Y为：
4 4 4
5 5 5

x ：

y ：

5. `[t,r] = cart2pol(x,y);`：将笛卡尔坐标(x,y)转换为极坐标(t,r)，其中`t`是角度(theta)，`r`是半径(radius)。

`cart2pol`是MATLAB中的一个函数，用于将笛卡尔坐标（x,y）转换为极坐标（r,t）。在这个转换过程中，

`x` 和 `y` 分别是笛卡尔坐标系中的横纵坐标，

而 `r` 和 `t` 是极坐标系中的半径（或长度）和角度。

转换过程如下：

1. **输入变量`x`和`y`：** - `x`：在笛卡尔坐标系中的横坐标。 - `y`：在笛卡尔坐标系中的纵坐标。

2. **输出`r`和`t`：** - `r`：从原点（0,0）到点（x,y）的距离，可以通过公式\[r = \sqrt{x^2 + y^2}\]计算得出。 - `t`：从正x轴到点（x,y）的线段相对于原点（0,0）的角度（以弧度为单位），可以用\[t = \text{atan2}(y,x)\]来计算。`atan2`函数可以处理x和y的符号以准确确定象限。

3. **操作过程：** - 对于每对输入的`(x,y)`坐标，`cart2pol`函数首先计算出距离`r`，这是通过将`x`的平方与`y`的平方相加，然后取平方根得到的。 - 接着，计算出角度`t`，这是通过用`(y,x)`调用`atan2`函数来完成的，这样可以根据`y`和`x`的符号确定正确的角度。

4. **结果：** - 结果是将每个`(x,y)`点转换为其对应的极坐标`(r,t)`，其中`r`表示点到原点的距离，而`t`表示从正x轴到该点的线段的角度。角度`t`是以弧度表示的，范围从`-π`到`π`。这个转换允许你将在x-y平面上定义的点转换为以原点为中心，以角度和距离表示位置的极坐标系中的点。这对于某些类型的数学和工程计算特别有用，比如那些涉及周期性或旋转对称性的计算。

6. `zernikeMatrices = [];`：初始化一个空数组，用来存储计算出的Zernike多项式矩阵。

7. 接下来的`for`循环遍历`zernike_indices`矩阵的每一行： - `for i = 1:size(zernike_indices,1)`：`for`循环开始，`size(zernike_indices,1)`获取`zernike_indices`行数，循环会迭代这么多次。

8. `zernikeMatrices(:,:,i) = zernike(r,t,zernike_indices(i,1),zernike_indices(i,2));`：在循环的每次迭代中，函数调用`zernike`，它通过输入的极坐标和Zernike多项式的n和m的值（来自`zernike_indices`矩阵的第i行）来计算对应的Zernike多项式矩阵，然后将结果矩阵赋值给`zernikeMatrices`的第i层。

9. `zernikeMatrices(:,:,i) = elliptical_crop(zernikeMatrices(:,:,i),1,mask);`：使用`elliptical_crop`函数对计算出来的Zernike多项式矩阵进行椭圆形的裁剪，其中使用`mask`来指定要裁剪的区域，`1`通常是这个函数的默认参数，表示需要裁剪的程度。以上这些步骤就创建了一系列与图像尺寸相同的Zernike多项式矩阵，可以进一步用于图像重构或波前分析。这些矩阵能帮助我们理解和补偿光学系统中的像差。

function cropped_im = elliptical_crop(im,crop_frac,mask)

% Performs elliptical cropping of greyscale image im. Extent of crop is
% defined by crop_frac. crop_frac of 1 corresponds to an ellipse whose
% major and minor axes are the dimensions of im. crop_frac of 0
% corresponds to an ellipse of 0 area.
%
%
% Example:
%
% % Display image 'trees.tif' with various croppings.
% trees = imread('trees.tif');
% figure, imagesc(trees)
% title('Or iginal image');
% figure, imagesc(elliptical_crop(trees,1));
% title('crop frac = 1');
% figure, imagesc(elliptical_crop(trees,0.7));
% title('crop frac = 0.7');
% figure, imagesc(elliptical_crop(trees,0.4));
% title('crop frac = 0.4');
% figure, imagesc(elliptical_crop(trees,0.1));
% title('crop frac = 0.1');
%
%
% Evan Czako, 8.14.2019.
% -------------------------------------------

if crop_frac < 0 || crop_frac > 1
error('crop_frac must have value between 0 and 1')
end

cropped_im = im;
center_x = (size(im,2)+1)/2;
center_y = (size(im,1)+1)/2;
radius_x = (size(im,2)-center_x)*crop_frac;
radius_y = (size(im,1)-center_y)*crop_frac;

for row = 1:size(im,1)
for col = 1:size(im,2)
if ((sqrt((row-center_y)^2/radius_y^2 + (col-center_x)^2/radius_x^2) > 1)||(mask(row,col)==1))
cropped_im(row,col) = nan;
end
end
end

% Necessary because of potential DIV/0 behavior
if radius_x == 0 || radius_y == 0
cropped_im = nan(size(cropped_im));
end

end

这段代码是一个用于对图像进行椭圆形剪裁的函数，具体步骤如下：

zernikeMatrices(:,:,i) = elliptical_crop(zernikeMatrices(:,:,i),1,mask);

这行代码是函数调用. 它将椭圆形剪裁应用于3D矩阵zernikeMatrices中的一个切片，并保存回原来的位置.

函数定义部分: function cropped_im = elliptical_crop(im,crop_frac,mask)

定义一个名为 elliptical_crop 的函数，该函数接受三个参数：待剪裁的图像 im，剪裁比例 crop_frac，和 mask. 函数返回被剪裁后的图像 cropped_im.

函数开始按照crop_frac(裁剪比例)进行裁剪. 如果crop_frac不在0到1之间，代码则抛出错误信息.

center_x = (size(im,2)+1)/2;
center_y = (size(im,1)+1)/2;

这两行计算图像中心的坐标点.

radius_x = (size(im,2)-center_x)*crop_frac;
radius_y = (size(im,1)-center_y)*crop_frac;

这两行计算根据剪裁比例定义的椭圆的x轴和y轴半径.

接下来的两个循环遍历图像的所有像素. 如果一个像素位于椭圆之外或者mask定义的区域内，那么将此像素设置为NaN.

如果任一半径为零，因为会导致0除以0的情况，那么这个函数会将整个图像全部设置为NaN.

这个函数最后将处理过後的图像返回.

这个函数允许你根据crop_frac的值来定义原始图像的一个椭圄区域进行剪裁，数据的裁剪区域之外的部分的值将会被设为NaN，而在mask区域内的值同样会被设为NaN. 换句话说，你可以选择保留图像的中心区域并去掉边缘和不需要的区域，具体的比例由你设定的crop_frac决定.

这样，这个函数最终输出与输入图像具有相同尺寸的一系列Zernike多项式矩阵，每个矩阵都与给定的n和m值对应。这个过程对于波前分析和图像重建等应用是非常有用的，因为Zernike多项式通常被用于描述光学系统中的波前失真。

其中上面那个代码用到了下面这个函数 zernike(r,t,n,m)

function zern = zernike(r,t,n,m)% Returns a matrix corresponding to the Zernike polynomial of radial% degree n and azimuthal degree m.% Matrices r and t correspond to the radii and angles of points on an % x-y grid satisfying 1>x>-1 and 1>y>-1; they are obtained via % cart2pol. See example for clarification.% Must satisfy n - m >= 0, n - m even, and n >= 0. n and m are % both integers.% Function elliptical_crop should be used in this context to ignore all % data outside of the unit circle.%% Example 1:%%     % Display the 'power' Zernike polynomial (n = 2, m = 0) over a %     % 100x100 grid.%%     x=linspace(-1,1,100);%     y=linspace(1,-1,100);%     [x,y] = meshgrid(x,y);%     [t,r] = cart2pol(x,y);%     z_power = zernike(r,t,2,0);%     figure, imagesc(elliptical_crop(z_power,1));%     colormap jet;%% Example 2:%%     % Display all Zernike polynomials up to radial degree 4 over a%     % 100x100 grid.% %     for n = 0:4%         for m = -n:2:n%             x=linspace(-1,1,100);%             y=linspace(1,-1,100);%             [x,y] = meshgrid(x,y);%             [t,r] = cart2pol(x,y);%             zern_mat = zernike(r,t,n,m);%             figure, imagesc(elliptical_crop(zern_mat,1));%             title(strcat('n = ',{' '},string(n),', m = ',{' '},string(m)));%             colormap jet;%         end%     end%%% Functions required for use: elliptical_crop, zernike_radial% See also: zernike_moments, zernike_recreation% % Evan Czako, 8.14.2019% -------------------------------------------if mod(n-m,2) == 1error('n-m must be even');endif n < 0error('n must both be positive')endif floor(n) ~= n || floor(m) ~= merror('n and m must both be integers')endif m < 0zern = -zernike_radial(r,n,-m).*sin(m*t);elsezern = zernike_radial(r,n,m).*cos(m*t);end
end

这个函数zernike的作用是生成对应于特定Zernike多项式的矩阵。Zernike多项式是一组正交和完备的多项式集，用于在单位圆内描述波前或光学系统中的像差。这里给出的函数能够根据给定的径向度（n）、方位角度（m）、半径（r）和角度（t）的值，计算对应的Zernike多项式值。每个Zernike多项式都由一对(n, m)值唯一确定，其中n代表径向阶数，m代表方位（角度）阶数。该函数详细说明如下：

函数参数：
- r: 单位圆内点的半径值。这是通过笛卡尔坐标到极坐标的转换得到的。
- t: 单位圆内点的角度值，即极角，也是通过坐标转换得到的。
- n: Zernike多项式的径向阶数。
- m: Zernike多项式的方位（角度）阶数。
参数条件：
- n - m 必须是偶数。这是因为Zernike多项式的定义要求阶数之差为偶数。
- n 必须是非负整数。因为径向阶数不能是负的。
- n 和 m 都必须是整数，这保证了多项式具有正确的形式。
函数功能：
- 函数首先检查给定的n和m是否符合Zernike多项式的基本条件（n - m是偶数，n是非负的，n和m都是整数）。如果不符合，会抛出错误。
- 当m为负数时，根据定义调用zernike_radial函数计算径向部分，并乘以sin(m*t)来得到最终的Zernike多项式。这对应于Zernike多项式的奇数模式。
- 当m为正数或零时，同样首先调用zernike_radial函数计算径向部分，再乘以cos(m*t)，得到对应的Zernike多项式。这对应于Zernike多项式的偶数模式。
结果：
- 返回值zern是计算出的Zernike多项式的矩阵，它基于给出的(r, t, n, m)参数。这个矩阵可在之后用于图像处理、波前分析等应用。
附加说明：
- 函数注释中提到的elliptical_crop函数，用于确保只处理位于单位圆内的数据，因为Zernike多项式只在单位圆内定义。
- zernike_radial函数负责计算Zernike多项式的径向部分，该部分只依赖于半径r和阶数n、m。

总的来说，这个zernike函数提供了一种计算特定Zernike多项式值的通用方法，这些值可以进一步用于描述和分析光学系统中的像差。

上面的函数又调用了zernike_radial(r,n,m)

具体函数实现如下：

function radial = zernike_radial(r,n,m)% Returns an array corresponding to the radial component of the Zernike% polynomial of radial degree n and azimuthal degree m evaluated at all% points specified by matrix r.% 'r' is a matrix containing distances from the center of itself,% corresponding to an x-y grid satisfying 1>x>-1 and 1>y>-1. See% example for clarification.% Must satisfy n - m >= 0, n - m even, n >= 0, and m >= 0. n and m are% both integers.% Function elliptical_crop should be used in this context to ignore all % data outside of the unit circle.% This function is utilized in function zernike to generate full % Zernike polynomials.% Utilizes recursive algorithm described by Chong et al.% Pattern Recognition, 36;731-742 (2003).%% % Example:%%     % Display the radial component of 'power' Zernike (n = 2, m = 0)  %     % over a 100x100 grid.%     x=linspace(-1,1,100);%     y=linspace(1,-1,100);%     [x,y] = meshgrid(x,y);%     [t,r] = cart2pol(x,y);%     zRad = zernike_radial(r,2,0);%     figure, imagesc(elliptical_crop(zRad,1));%     colormap jet;%%% Functions required for use: elliptical_crop% See also: zernike, zernike_moments, zernike_recreation% % Evan Czako, 8.14.2019% -------------------------------------------if mod(n-m,2) == 1error('n-m must be even');endif n < 0 || m < 0error('n and m must both be positive in radial function')endif floor(n) ~= n || floor(m) ~= merror('n and m must both be integers')endif n == mradial = r.^n;elseif n - m == 2radial = n*zernike_radial(r,n,n)-(n-1)*zernike_radial(r,n-2,n-2);elseH3 = (-4*((m+4)-2)*((m+4)-3)) / ((n+(m+4)-2)*(n-(m+4)+4));H2 = (H3*(n+(m+4))*(n-(m+4)+2)) / (4*((m+4)-1))  +  ((m+4)-2);H1 = ((m+4)*((m+4)-1) / 2)  -  (m+4)*H2  +  (H3*(n+(m+4)+2)*(n-(m+4))) / (8);radial = H1*zernike_radial(r,n,m+4) + (H2+H3 ./ r.^2).*zernike_radial(r,n,m+2);% Fill in NaN values that may have resulted from DIV/0 in prior% line. Evaluate these points directly (non-recursively) as they% are scarce if present.if sum(sum(isnan(radial))) > 0[row, col] = find(isnan(radial));c=1;while c<=length(row)x = 0;for k = 0:(n-m)/2%((-1)^k*factorial(n-k))/(factorial(k)*factorial((n+m)/2-k)*factorial((n-m)/2-k))*0^(n-2*k);x = x + ((-1)^k*factorial(n-k))/(factorial(k)*factorial((n+m)/2-k)*factorial((n-m)/2-k))*0^(n-2*k);endradial(row(c),col(c)) = x;c=c+1;endendendend

这是MATLAB中定义的一个函数，名为`zernike_radial`。它用于计算并返回Zernike多项式的径向分量的值。Zernike多项式在波前分析中常用于表征光学系统的像差。

输入参数介绍：

- `r`：该矩阵包含其自身中心的距离，对应于满足1>x>-1和1>y>-1的x-y网格上的所有点。

- `n`：径向次数，非负整数，定义了多项式在径向上的复杂度。

- `m`：方位次数，非负整数，定义了多项式的角向变化，必须小于等于`n`，并且`n-m`是偶数。

函数开始时，会检查输入的`n`和`m`是否满足特定条件，比如`n - m` 是否为偶数，`n`和`m`都是非负整数以及都是整数。

当`n==m`时，径向分量是简单的`r^n`。

当`n - m == 2`时，会使用递归关系`n*zernike_radial(r,n,n)-(n-1)*zernike_radial(r,n-2,n-2)`来计算径向分量。

当`n - m > 2`时，函数使用Chong等人在2003年提出的一个递归算法来计算径向分量。在这个递归关系中，使用了`H1`，`H2`和`H3`这三个辅助变量来简化计算。如果在计算过程中遇到0作为除数的情况产生`NaN`值，函数将直接计算这些特定点的径向分量，避免递归。

这个函数还提供了一种直接计算`NaN`值的方式。当调用这个函数来生成整个Zernike多项式时，通常需要搭配使用`elliptical_crop`函数来忽略单位圆外的数据点。

示例部分展示了如何显示一个100x100网格上Zernike多项式径向分量`n=2`，`m=0`（`Z^2_0`）的情形。使用`meshgrid`创建坐标网格，`cart2pol`转换为极坐标，然后调用`zernike_radial`计算径向分量，并使用`imagesc`和颜色图`jet`来显示结果。

注意：这个注释是为了给使用这个函数的人提供帮助，并且假定用户具备MATLAB编程背景和光学或者波前分析领域的知识。

这个函数中递归算法的具体实现涉及到Zernike多项式的径向分量的递归计算，尤其是当无法直接使用简单公式时。在该递归计算过程中，为了避免除以零而产生的NaN值，采取了特殊措施。下面详细解释这一过程：

首先，让我们更深入地看看H1，H2，和H3这三个辅助变量的计算方法，以及它们如何用于计算径向分量radial的递归过程：

辅助变量的计算：
- H1、H2和H3这三个变量通过一系列的计算得到，其计算依赖于输入的n和m值。这三个变量的计算涉及多个阶段，包括利用n和m的变化来调整系数，这有助于导出递归计算的公式。
递归计算径向分量：
- 径向分量radial的递归计算依赖于这些辅助变量。具体地说，递归公式会结合使用H1，H2，和除以r^2的H3，以此来计算给定n，m值下的径向分量。这个过程可能涉及到对更高n，m值的zernike_radial函数的调用。
避免除以零的策略：
- 在计算过程中，当存在除以零的风险时（即r为零时，r^2也为零，导致公式中的除法操作不被定义），可能会产生NaN值。为了处理这个问题，函数内部首先检查是否有任何计算结果是NaN。如果存在这样的结果，那么函数会直接计算这些特定点的径向分量，而不使用递归方法。
- 直接计算这些特定点的径向分量通常涉及到避免任何可能导致除以零的操作。在这个函数中，通过遍历所有含有NaN值的元素，并使用一个针对值为零时的特殊情况下的公式来直接计算它们的径向分量。
特殊情况下的直接计算：
- 直接计算的方法是用一个for循环，计算一个特定n，m下径向分量的和。其中，利用了一个关于Zernike多项式径向分量的基本公式，该公式不涉及除法，因此避免了除以零的问题。这种方法通常在遇到r为零的情况下使用，这时递归方法可能失效。

总的来说，这个函数的设计考虑了Zernike多项式径向分量计算中可能出现的特殊情形和计算挑战。通过使用递归算法的同时也实施了适当的错误处理和边界情况处理，确保了函数的鲁棒性和准确性。

if sum(sum(isnan(radial))) > 0 [row, col] = find(isnan(radial)); c=1; while c<=length(row) x = 0; for k = 0:(n-m)/2 %((-1)^k*factorial(n-k))/(factorial(k)*factorial((n+m)/2-k)*factorial((n-m)/2-k))*0^(n-2*k); x = x + ((-1)^k*factorial(n-k))/(factorial(k)*factorial((n+m)/2-k)*factorial((n-m)/2-k))*0^(n-2*k); end radial(row(c),col(c)) = x; c=c+1; end end

这段代码是在处理zernike_radial函数中出现的特殊情况，即计算中出现除以0的情况，从而导致某些径向分量值为NaN。其主要目的是针对这些出现NaN值的特定元素，直接计算其径向分量值，而不是通过递归方法来计算。这种直接计算方法适用于那些因为r等于0导致除以r的运算无效的情况。

原理解释：

检测NaN: 首先，通过sum(sum(isnan(radial)))计算来检查在radial矩阵中是否有NaN值。如果有，则表明存在一些由于除以0造成的无效计算。

确定位置： 使用[row, col] = find(isnan(radial))来找出这些NaN值的具体位置。这里，row和col分别存放了这些NaN值在radial矩阵中的行索引和列索引。

直接计算： 对于找到的每一个NaN位置，使用一个for循环直接计算对应位置的径向分量值。计算的公式利用到了Zernike多项式的径向分量的一个基本表达式。这个表达式是一个关于n和m的求和表达式，对所有k从0到(n-m)/2进行迭代计算，其中每一项都是利用了阶乘来计算。重要的是在这个表达式中，由于乘以了0^(n-2*k)，当n-2*k不为0时，这一项结果为0；当n-2*k等于0时，这个表达式则简化为一个与r无关的常数值计算，因而避免了除以0的问题。

更新径向分量值： 计算得到的径向分量值被赋值给radial矩阵中原本是NaN的相应位置，从而修正了之前的计算问题。

总之，这段代码通过直接计算的方法处理了计算中的特殊情况，确保所有的径向分量都能获得有效的数值，从而使zernike_radial函数可以正确完成对Zernike多项式径向分量的计算。

分析公式的各个部分：

x = x + ((-1)^k*factorial(n-k))/(factorial(k)*factorial((n+m)/2-k)*factorial((n-m)/2-k))*0^(n-2*k);

(-1)^k：这表示迭代项的符号将根据k的值交替变化。
factorial(n-k)、factorial(k)、factorial((n+m)/2-k)、factorial((n-m)/2-k)：这些是阶乘项，用于计算组合数，是求和公式的一部分。
0^(n-2*k)：这是公式中的关键，用于处理r=0时的情形。当n-2*k不为0时，这一项为0，因此整个迭代项贡献为0；当n-2*k为0时，这一项为1，意味着该迭代项将按照其它因子计算而不被忽略。这确保了在r=0的情况下，公式只计算有意义的项（即当n是2*k的倍数的情形）。

$x=x+\frac{(-1)^k\cdot (n-k))!}{k!\cdot\left(\frac{n+m}2-k\right)!\cdot\left(\frac{n-m}2-k\right)!}\cdot0^{(n-2k)}$

因此，这行代码的目的是在r=0时，根据Zernike多项式的定义直接计算中心点（或其他相应点）的径向分量值。这种计算方法特别适用于处理递归计算中出现的除0情况，如在进行图像分析或光学波前分析时计算中心点的情况。

上面这部分算是结束了。上面主要征对 zernike_mats 进行展开的讨论，后续进行zernike_moments 的讨论

zernikeMats = zernike_mats(zernikePhase, indices, newmask); %计算相位数据在这些 Zernike 多项式上的投影，即计算 Zernike 模式。
a_coeffs = zernike_moments(zernikePhase, indices, newmask); % zernike36项系数 a_coeffs 函数来计算相位数据的 Zernike 系数。

function [zernMoments] = zernike_moments(im,indices,mask)% Performs a least squares fit/decomposition of matrix im using Zernike % polynomials of radial and azimuthal degrees described by indices, % where the first column of indices contains values of n and the second % column contains values of m.% Return a column vector containing the Zernike coefficients (or % moments) that best describe matrix im.%% Example 1:% %     % Create some function z(x,y) over a 100x100 grid.%     x=linspace(-1,1,100);%     y=linspace(1,-1,100);%     [x,y] = meshgrid(x,y);%     z = x.^2+y.^3;%     figure,imagesc(elliptical_crop(z,1));%     colormap jet;%     title('z(x,y)');% %     % Specify desired indices of Zernike polynomials to be used in the%     % decomposition. In this case, include all polynomials up to and including %     % radial degree 3.% %     indices = [];%     for n = 0:3%         for m = -n:2:n%             indices = [indices; n m];%         end%     end% %     % Display the coefficients associated with each Zernike polynomial used in%     % the fit and their corresponding (n,m) indices.%     disp(['    Coeff     ', 'n        ', 'm']);%     moments = zernike_moments(z,indices);%     disp([moments indices]);%%% Example 2:%%     % Some evidence that this function works.%     % Create list of all (n,m) indices up to n = 3:%     indices = [];%     for n = 0:3%         for m = -n:2:n%             indices = [indices; n m];%         end%     end% %     % Create a 100x100 empty matrix for demonstration.%     zernike_sum_1 = zeros(100,100);% %     % 3-dimensional matrix zern_mats contains the first 10 Zernike polynomials%     % as 100x100 grids.%     zern_mats = zernike_mats(zernike_sum_1,indices);% %     % Add an equal amount of each of the first 10 Zernike polynomials to%     % zernike_sum_1. This way, the least-squares fit should return the same%     % coefficient for each of the polynomials used in the fit.%     for i = 1:size(zern_mats,3)%         zernike_sum_1 = zernike_sum_1 + zern_mats(:,:,i);%     end% %     % The coefficients corresponding to each (n,m) pair used in the fit are the%     % same. In this case, they each have a value of 1.%     disp(['    Coeff     ', 'n        ', 'm']);%     disp([zernike_moments(zernike_sum_1,indices) indices]);%%     % In case you were wondering what an equal amount of Z1-10 looks like:%     figure,imagesc(zernike_sum_1);%     colormap jet;%%% Example 3:%%     % Some more evidence that this function works.%     % Create list of all (n,m) indices up to n = 3:%     indices = [];%     for n = 0:3%         for m = -n:2:n%             indices = [indices; n m];%         end%     end% %     % Create a 100x100 empty matrix for demonstration.%     zernike_sum_2 = zeros(100,100);% %     % 3-dimensional matrix zern_mats contains the first 10 Zernike polynomials%     % as 100x100 grids.%     zern_mats = zernike_mats(zernike_sum_2,indices);% %     % Add an increasing amount of each of the first 10 Zernike polynomials to%     % zernike_sum_2. This way, the least-squares fit should return increasing %     % values for the coefficients of each polynomial used in the fit.%     for i = 1:size(zern_mats,3)%         zernike_sum_2 = zernike_sum_2 + i*zern_mats(:,:,i);%     end% % %     % The coefficients corresponding to each (n,m) pair used in the fit increase %     % incrementally from 1 to 10.%     disp(['    Coeff     ', 'n        ', 'm']);%     disp([zernike_moments(zernike_sum_2,indices) indices]);% %     % In case you were wondering what this thing looks like:%     figure,imagesc(zernike_sum_2);%     colormap jet;%%% Functions required for use: zernike_mats, zernike, zernike_radial,% elliptical_crop% See also: zernike_recreation%% Evan Czako, 8.14.2019% -------------------------------------------dim=size(im);NumCols = dim(2);NumRows = dim(1);zernikeMats = zernike_mats(im,indices,mask);
%     for i=1:36
%         zernikeMats(:,:,i)=mask.*zernikeMats(:,:,i);
%     endz_mats_reshaped = [];for i = 1:size(indices,1)z_mats_reshaped(:,i) = reshape(zernikeMats(:,:,i),NumCols*NumRows,1);endim(isnan(im))=0;   z_mats_reshaped(isnan(z_mats_reshaped))=0;image_reshaped = reshape(im,NumCols*NumRows,1);a = (z_mats_reshaped.'*z_mats_reshaped)^-1*(z_mats_reshaped.'*image_reshaped);zernMoments = a;end% function dst = Reshape(src,row,col)
% num_nan = sum(isnan(src));
% 
% 
% end

这段Matlab代码定义了一个函数 zernike_moments, 它使用Zernike多项式对一个二维矩阵（通常是一幅图像）进行最小二乘拟合，提取Zernike矩（也称为Zernike系数）。这个函数非常有用于图像分析，例如在图案识别或光波前分析中。

参数解释：