matlab练习程序（泽尼克多项式拟合）

泽尼克多项式是一个正交多项式，分为奇偶两类。

奇多项式：

偶多项式：

其中：

这里fai为方位角，范围[0-2pi]；p为径向距离，范围[0,1]；n-m大于等于0；

如果n-m=0，则R=0。

根据不同的m和n值，可以得到不同的多项式，用j表示不同的多项式，通常称为Noll序列：

n,m	0,0	1,1	1,−1	2,0	2,−2
j	1	2	3	4	5
n,m	2,2	3,−1	3,1	3,−3	3,3
j	6	7	8	9	10
n,m	4,0	4,2	4,−2	4,4	4,−4
j	11	12	13	14	15
n,m	5,1	5,−1	5,3	5,−3	5,5
j	16	17	18	19	20

Noll序列前15项泽尼克多项式为：

1	$1$
2	$2\rho cos\theta $
3	$2\rho sin\theta $
4	$\sqrt{3}(2\rho ^{2}-1)$
5	$\sqrt{6}\rho ^{2}sin2\theta $
6	$\sqrt{6}\rho ^{2}cos2\theta $
7	$\sqrt{3}(3\rho ^{3}-2\rho )sin\theta $
8	$\sqrt{3}(3\rho ^{3}-2\rho )cos\theta $
9	$\sqrt{8}\rho ^{3}sin3\theta $
10	$\sqrt{8}\rho ^{3}cos3\theta $
11	$\sqrt{5}(6\rho^{4}-6\rho ^{2}+1)$
12	$\sqrt{10}(4\rho^{4}-3\rho ^{2})cos2\theta $
13	$\sqrt{10}(4\rho^{4}-3\rho ^{2})sin2\theta $
14	$\sqrt{10}\rho^{4}cos4\theta $
15	$\sqrt{10}\rho^{4}sin4\theta $

下面用15阶泽尼克多项式拟合了一个高斯函数。

matlab代码如下：

clear all;close all;clc;

[x,y]=meshgrid(-100:100);
sigma=50;
target_img = 10000*(1/(2*pi*sigma^2))*exp(-(x.^2+y.^2)/(2*sigma^2));

B = target_img(:);
srcA = zeros(length(B),15);
for i=1:15
    [x,y] = meshgrid(-1:0.01:1);
    p = sqrt(x.^2+y.^2);
    theta = atan2(y,x);
    
    img = zernike_polynomials(i,p,theta);
    img(p>1) = nan;
    subplot(3,5,i);
    imagesc(img);
    srcA(:,i) = img(:);
    ax=gca;
    ax.XAxis.Visible='off';
    ax.YAxis.Visible='off';
end
A = srcA;

B(isnan(srcA(:,1))) = [];
A(isnan(srcA(:,1)),:) = [];

x = inv(A'*A)*A'*B;

fit_img = zeros(size(target_img));
for i=1:15
    fit_img = fit_img + x(i)*reshape(srcA(:,i),size(target_img));
end
figure;
mesh(target_img);
hold on;
surf(fit_img);

function re = zernike_polynomials(n,p,fai)
switch n
    case 1
        re = ones(size(p));
    case 2
        re = 2*p.*sin(fai);
    case 3
        re = 2*p.*cos(fai);
    case 4
        re = sqrt(6)*p.^2.*sin(2*fai);
    case 5
        re = sqrt(3)*(2*p.^2-1);
    case 6
        re = sqrt(6)*p.^2.*cos(2*fai);
    case 7
        re = sqrt(8)*p.^3.*sin(3*fai);
    case 8
        re = sqrt(8)*(3*p.^3-2*p).*sin(fai);
    case 9
        re = sqrt(8)*(3*p.^3-2*p).*cos(fai);
    case 10
        re = sqrt(8)*p.^3.*cos(3*fai);
    case 11
        re = sqrt(10)*p.^4.*sin(4*fai);
    case 12
        re = sqrt(10)*(4*p.^4-3*p.^2).*sin(2*fai);
    case 13
        re = sqrt(5)*(6*p.^4-6*p.^2+1);
    case 14
        re = sqrt(10)*(4*p.^4-3*p.^2).*cos(2*fai);
    case 15
        re = sqrt(10)*p.^4.*cos(4*fai);
end
end

结果如下:

泽尼克图像：

拟合结果对比，浅色曲面为原高斯曲面，深色为拟合曲面：

参考了维基百科。