Разложение волнового фронта по полиномам Цернике
Итак, нагулявшись в отпуске я таки решил вернуться к своим баранам.
Прежде, чем начать писать что-то для разложения нормалей к волновому фронту по упомянутым раньше функциям Жао, я решил попробовать обычное разложение. Для этого с сайта 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 для вычисления "полиномов Жао".
Прежде, чем начать писать что-то для разложения нормалей к волновому фронту по упомянутым раньше функциям Жао, я решил попробовать обычное разложение. Для этого с сайта 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 для вычисления "полиномов Жао".