Разложение волнового фронта по полиномам Цернике

[eddy_em]

Итак, нагулявшись в отпуске я таки решил вернуться к своим баранам.
Прежде, чем начать писать что-то для разложения нормалей к волновому фронту по упомянутым раньше функциям Жао, я решил попробовать обычное разложение. Для этого с сайта mathworks был сворован исходник функции zernfun2 (которая вычисляет полиномы), а также deco (которая выполняет декомпозицию).
Функцию deco пришлось маленько подправить (она вычисляла коэффициенты на основе двойного численного интегрирования, что крайне медленно). Но получилось довольно-таки сносно.

Итак, вот что у меня получилось с декомпозицией:

% deco1.m - decomposition
N=101; % dimension of coordinates grid !!!MUST be odd!!!
p = 0:100; % indexes of Zernike functions to be used in the decomposition (later should determine automatically)
x = linspace(-1,1,N); y = linspace(-1,1,N);
[x,y] = meshgrid(x,y); % make simple equidistant coordinate grid
[t,r] = cart2pol(x,y); % polar coordinates corresponding to cartesian grid
idx = find(r>1); % indexes of external elements
x(idx) = NaN; y(idx) = NaN; t(idx) = NaN; r(idx) = NaN; % to be on a safe side, clear bad data
F = sin(x*4).*cos(y*10)/4-tan(x.^2+y.^2);
%F = x.^3-y.^2; % any surface to decompose
surf(x,y,F); shading interp 
title('Function to be decomposed'); axis square; print -dpng composed.png; close
n = length(p); 
z = zernfun2(p,r(:),t(:),'norm');
zz = reshape(z,N,N,length(p)); 
for k = 1:n
    z = zz(:,:,k);
    subplot(ceil(n^.5),ceil(n^.5),k)
    surf(x,y,z), shading interp
    set(gca,'XTick',[],'YTick',[])
    axis square
    title(['Z_' num2str(p(k))]);
end
print -dpng zernike.png; close
fprintf("Make decomposition\n"); fflush(1);
Q = [];
idx = find(~isnan(F)); % indexes of inner elements
for i = 1:n
    fprintf(1, "K = %d\n", i); fflush(1);
    aSum = (zz(:,:,i).*F)(idx);
    aNorm = (zz(:,:,i).^2)(idx);
    Q(i) = sum(sum(aSum)) / sum(sum(aNorm));
end
fprintf("Q = "), disp(Q); fflush(1);
Z = zeros(N);
for k = 1:n
    Z=Z+Q(k)*zz(:,:,k);
    fprintf(1, " K = %d, std(Z-F) = %f\n", k, std((Z-F)(idx))); fflush(1);
    subplot(ceil(n^.5),ceil(n^.5),k)
    surf(x,y,Z), shading interp
    set(gca,'XTick',[],'YTick',[])
    axis square
    if k==1
        title(['Z_' num2str(p(k))]);
    else
        title(['Z_0 to Z_' num2str(p(k))]);
    end
end
print -dpng decompose.png; close
surf(x,y,Z-F); shading interp;
title('Residual error'); axis square
print -dpng error.png; close

N и p - изменяющиеся параметры. Так же как и функция F.
Для начала я взял функцию F=x^3-y^2 (она на рисунке выше, функция эта и была в примере). Так как функция простая, здесь вполне хватит 16 первых полиномов. Вот они:

После разложения получилось вот что:

Каждый N-й график - это результат суммирования полиномов с 0 по N-1, помноженных на соответствующие коэффициенты.
Получилось вполне неплохо, несмотря на то, что интегрирование по сути используется методом прямоугольников. Вот - остаточная ошибка:

Теперь попробуем более сложную функцию F = sin(x*4).*cos(y*10)/4-tan(x.^2+y.^2):

Разложение по 16 полиномам получилось плохое. Коэффициенты: Q = -1.0893e+00 -1.1108e-16 1.1683e-02 -1.9055e-18 -7.3561e-01 1.5183e-15 -2.8639e-17 -1.4932e-16 6.6768e-03 -1.3585e-02 8.0573e-18 1.3411e-17 -1.1139e-01 1.6325e-16 1.6105e-02 -1.9671e-18 -4.7822e-17, среднеквадратичное отклонение:

 K = 1, std(Z-F) = 0.438467
 K = 2, std(Z-F) = 0.438467
 K = 3, std(Z-F) = 0.438417
 K = 4, std(Z-F) = 0.438417
 K = 5, std(Z-F) = 0.139009
 K = 6, std(Z-F) = 0.139009
 K = 7, std(Z-F) = 0.139009
 K = 8, std(Z-F) = 0.139009
 K = 9, std(Z-F) = 0.138957
 K = 10, std(Z-F) = 0.138750
 K = 11, std(Z-F) = 0.138750
 K = 12, std(Z-F) = 0.138750
 K = 13, std(Z-F) = 0.122085
 K = 14, std(Z-F) = 0.122085
 K = 15, std(Z-F) = 0.122385
 K = 16, std(Z-F) = 0.122385
 K = 17, std(Z-F) = 0.122385

Вот такая вот вышла остаточная ошибка:

Оно и понятно: первые 16 полиномов имеют довольно низкую частоту (это как и с вейвлетами - аналогия схожая). Поэтому попробуем увеличить количество полиномов.
Увеличение до 36 ничего путного не дало: среднеквадратичное отклонение упало лишь до 0.09.
Увеличение до сотни дало лучшие результаты:

ошибка уменьшилась до 0.0288 на N=77 и возрасла до 0.0324 на N=100.

Итак, несмотря на такое огрубление вычисления коэффициентов декомпозиции, результат получается вполне приличным. Остается лишь переписать zernfun2 для вычисления "полиномов Жао".