![[personal profile]](https://www.dreamwidth.org/img/silk/identity/user.png){kind=small}
eddy_emНаконец-то нашел я время вернуться к обработке гартманнограмм.
Гартманнограммы я сгенерировал при помощи уже описанной модели, из нее получил значения нормалей к волновому фронту, преобразовал их в полярные координаты и далее стал пытаться восстановить волновой фронт так, чтобы получить нечто похожее на оригинальный.
Оригинальный волновой фронт имеет следующий вид:
Сначала — вкратце расскажу о простейших (и → грубейших) методах восстановления волнового фронта.
Это — простейший (и наименее точный) метод восстановления волнового фронта. В некоторой точке, являющейся начальной (например, на краю волнового фронта), величина волнового фронта считается равной нулю. Далее осуществляется переход вдоль некоей линии, соединяющей изображения пятен на гартманнограмме, к следующей точке. В ней волновой фронт вычисляется как линейная экстраполяция наклона фронта относительно предыдущей точки по значению наклона в этой точке:
где s - проекция расстояния между точками Фi и Фi-1 на ось ζ (которая может быть и нелинейной — например, профиль вдоль окружности). Компоненты Фζ вычисляются для двух ортогональных осей α и β, а затем суммированием получается линейное приближение к волновому фронту:
Саусвелл (W. Southwell) предложил простой метод зонального восстановления волнового фронта. Метод заключается в том, что волновой фронт в каждой точке вычисляется на основе данных о наклоне фронта во всех соседних точках:
где Фn,m — величина волнового фронта в пятне (n,m);
Δi,j — расстояние между пятнами (n,m) и (n+i,m+j),
In+i, m+j — интенсивность излучения в пятне (n+i,m+j).
В качестве весового коэффициента здесь используется интенсивность излучения в пятне, однако, в случае необходимости весами можно пренебречь.
Этот алгоритм — итерационный, итерации прекращаются при достижении желаемой точности (но не более количества пятен на гартманнограмме).
Итак, для начала я составил такой алгоритм:
В результате на десяти итерациях получился вот такой вот ужас:
и действительно: ведь в итоге получается, что все, вычисленное на итерации с увеличением переменной, нейтрализуется на итерации с ее уменьшением.
Как вариант, можно опустить вычисления в "негативную" сторону. Т.е. закомментировать строчки с Phi(nRplus) и Phi(nTminus). На десяти итерациях получим вот что:
После обработки при помощи функции
получилась вот такая карта поверхности зеркала (не имеющая ничего общего с оригиналом):
Может быть, сделать 100 итераций? Делаем:
Получается вот такая картинка:
Все еще непонятно что, не имеющее отношения к реальности. Однако, можно будет попытаться поковыряться в этом алгоритме: вдруг из него выйдет что-то путное. Но это оставлю на потом, если восстановление полиномами Цернике и преобразованиями Фурье тоже не сработают (в принципе, я мог допустить ошибку на этапе преобразования нормалей к волновому фронту).
Сделаем небольшие изменения в функции iterate_Phi, чтобы она считала компоненты независимо, а потом уже складывала их, получая "восстановленный" волновой фронт. Если количество итераций будет равно единице, этот метод как раз и будет чем-то вроде линейного интегрирования.
Результаты получаются уж совсем страшными:
"Восстановленный" волновой фронт
"Восстановленный" (синим) и идеальный (красным) волновой фронт
Попробуем сравнить "фронты" на 10 и на 1:
С каждой итерацией наблюдается все большее и большее расхождение.
Похоже, этот способ уж совсем никуда не годится.
Из-за того, что наша гартманнограмма имеет не прямоугольную, а концентрическую сетку, возникает уйма проблем, связанных с ее обработкой. Простейшие классические методы дают сильные погрешности из-за краевых эффектов.
Возможно, дифференциальный метод (при помощи полиномов Цернике) будет более надежным.
Гартманнограммы я сгенерировал при помощи уже описанной модели, из нее получил значения нормалей к волновому фронту, преобразовал их в полярные координаты и далее стал пытаться восстановить волновой фронт так, чтобы получить нечто похожее на оригинальный.
Оригинальный волновой фронт имеет следующий вид:
Сначала — вкратце расскажу о простейших (и → грубейших) методах восстановления волнового фронта.
Линейное интегрирование
Это — простейший (и наименее точный) метод восстановления волнового фронта. В некоторой точке, являющейся начальной (например, на краю волнового фронта), величина волнового фронта считается равной нулю. Далее осуществляется переход вдоль некоей линии, соединяющей изображения пятен на гартманнограмме, к следующей точке. В ней волновой фронт вычисляется как линейная экстраполяция наклона фронта относительно предыдущей точки по значению наклона в этой точке:
где s - проекция расстояния между точками Фi и Фi-1 на ось ζ (которая может быть и нелинейной — например, профиль вдоль окружности). Компоненты Фζ вычисляются для двух ортогональных осей α и β, а затем суммированием получается линейное приближение к волновому фронту:
Фi = Фiα + Фiβ.
Восстановление Саусвелла
Саусвелл (W. Southwell) предложил простой метод зонального восстановления волнового фронта. Метод заключается в том, что волновой фронт в каждой точке вычисляется на основе данных о наклоне фронта во всех соседних точках:
где Фn,m — величина волнового фронта в пятне (n,m);
Δi,j — расстояние между пятнами (n,m) и (n+i,m+j),
In+i, m+j — интенсивность излучения в пятне (n+i,m+j).
В качестве весового коэффициента здесь используется интенсивность излучения в пятне, однако, в случае необходимости весами можно пренебречь.
Этот алгоритм — итерационный, итерации прекращаются при достижении желаемой точности (но не более количества пятен на гартманнограмме).
Пробую метод Саусвелла
Итак, для начала я составил такой алгоритм:
function Phi = iterate_Phi(N, R, theta, dR, dT)
%
% Phi = iterate_phi(N, R, theta, dR, dT)
% восстановление фазы итерационным методом Саусвелла
% Вход:
% N - количество итераций
% R, theta - координаты пятен в полярной нормированной СК
% dR, dT - производные фазы по радиальным координатам
% (R, theta, dR, dT - результат команды conv2polar)
% Выход:
% Phi - найденная фаза
%
% для начала получим индексы точек, расположенных по соседству с данной
% методика сильно привязана к маске; если какие-то точки будут пропущены,
% будем иметь сложности с расчетами
idxRplus = 33:256; % индексы точек с бОльшим номером по R
nRplus = 1:224; % номера точек для данных индексов
idxRminus = 1:224; % индексы точек с меньшим номером по R
nRminus = 33:256;
nnRminus = 225:256; % номера точек, которые есть лишь для Rminus
idxTplus = []; % индексы точек с бОльшим номером по theta
idxTminus = []; % индексы точек с меньшим номером по theta
for i = 0:7
idxTminus = [idxTminus ([32 1:31]+i*32)];
idxTplus = [idxTplus ([2:32 1]+i*32)];
endfor
nTplus = nTminus = 1:256;
% теперь вычисляем расстояния между соседними точками
dRplus = R(nRplus) - R(idxRplus);
dRminus = R(nRminus) - R(idxRminus);
dTplus = theta(nTplus) - theta(idxTplus) ;
do
tmp = find(dTplus < 0);
dTplus(tmp) += 2*pi;
until(isempty(tmp))
dTplus .*= (R(nTplus)+R(idxTplus))/2;
dTminus = theta(nTminus) - theta(idxTminus);
do
tmp = find(dTminus > 0);
dTminus(tmp) -= 2*pi;
until(isempty(tmp))
dTminus .*= (R(nTminus)+R(idxTminus))/2;
Phi = zeros(1, 256);
% прибавляем к фазе средние отклонения
while(N-- > 0)
% проходим изнутри наружу, вычисляя точки с бОльшим R
Phi(nRminus) = Phi(idxRminus) + dRminus.*(dR(idxRminus)+dR(nRminus))/2;
% проходим снаружи вовнутрь, вычисляя точки с меньшим R
Phi(nRplus) = Phi(idxRplus) + dRplus.*(dR(idxRplus)+dR(nRplus))/2;
% проходимся по угловым проекциям
Phi(nTplus) = Phi(idxTplus) + dTplus.*(dT(idxTplus)+dT(nTplus))/2;
Phi(nTminus) = Phi(idxTminus) + dTminus.*(dT(idxTminus)+dT(nTminus))/2;
endwhile
endfunction
В результате на десяти итерациях получился вот такой вот ужас:
и действительно: ведь в итоге получается, что все, вычисленное на итерации с увеличением переменной, нейтрализуется на итерации с ее уменьшением.
Как вариант, можно опустить вычисления в "негативную" сторону. Т.е. закомментировать строчки с Phi(nRplus) и Phi(nTminus). На десяти итерациях получим вот что:
После обработки при помощи функции
function [delta_Phi img] = make_mirror_map(Phi, dR, R, phi, imsize)
%
% img = make_mirror_map(Phi, size)
% построение "карты" поверхности параболического зеркала
% Вход:
% Phi - восстановленная тем или иным способом фаза
% dR, X,Y - радиальные производные волнового фронта в точках (R,phi)
% size - размер карты (необязательный параметр, по умолчанию - 256)
% Выход:
% img - сама карта (квадратное изображение, отклонения - в единицах Phi,
% апертура расчитывается по макс. радиусу пятен)
%
if(nargin == 4) imsize = 256; endif
sz = max(size(Phi));
% на случай, если dR и (X,Y) еще содержат сведения о маркерах:
R = R(1:sz); phi = phi(1:sz);
dR = dR(1:sz);
% вычисляем нормированные координаты:
X = R .* cos(phi); Y = R .* sin(phi);
% так как уравнение параболического зеркала - z=Kr^2, dz/dr = 2Kr,
% то K = (dz/dr)/(2r) ==> z=[(dz/dr)/r]*r^2/2
Phi_ideal = median(dR./R).*R.^2/2;
% отклонение восстановленной поверхности от идеальной
% (методы, все-таки, точны до произвольной аддитивной константы)
dZ = median(Phi-Phi_ideal(1:sz))
delta_Phi = Phi - dZ - Phi_ideal; % отклонения фазы
%plot([Phi', Phi_ideal'])
printf("statistics: std = %g, dZ_max = %g, dZ(ideal)_max = %g\n", ...
std(delta_Phi), max(max(Phi)), max(max(Phi_ideal)));
[x1,y1] = meshgrid(([1:imsize] - imsize/2)* max(R)/(imsize/2)); % сетка для интерполяции
img = griddata(X(1:256), Y(1:256), delta_Phi, x1, y1);
endfunction
получилась вот такая карта поверхности зеркала (не имеющая ничего общего с оригиналом):
Может быть, сделать 100 итераций? Делаем:
Phi100 = iterate_Phi(100, R, phi, dR, dphi);
[dPhi100 img] = make_mirror_map(Phi100, dR, R, phi);
imagesc(img)
axis equal
print -color -dpng fig.png;
Получается вот такая картинка:
Все еще непонятно что, не имеющее отношения к реальности. Однако, можно будет попытаться поковыряться в этом алгоритме: вдруг из него выйдет что-то путное. Но это оставлю на потом, если восстановление полиномами Цернике и преобразованиями Фурье тоже не сработают (в принципе, я мог допустить ошибку на этапе преобразования нормалей к волновому фронту).
Пробуем линейное интегрирование
Сделаем небольшие изменения в функции iterate_Phi, чтобы она считала компоненты независимо, а потом уже складывала их, получая "восстановленный" волновой фронт. Если количество итераций будет равно единице, этот метод как раз и будет чем-то вроде линейного интегрирования.
function Phi = iterate_Phi_Linear(N, R, theta, dR, dT)
%
% Phi = iterate_Phi_Linear(N, R, theta, dR, dT)
% восстановление фазы линейной интерполяцией итерациями
% Вход:
% N - количество итераций
% R, theta - координаты пятен в полярной нормированной СК
% dR, dT - производные фазы по радиальным координатам
% (R, theta, dR, dT - результат команды conv2polar)
% Выход:
% Phi - найденная фаза
%
% для начала получим индексы точек, расположенных по соседству с данной
% методика сильно привязана к маске; если какие-то точки будут пропущены,
% будем иметь сложности с расчетами
idxRplus = 33:256; % индексы точек с бОльшим номером по R
nRplus = 1:224; % номера точек для данных индексов
idxRminus = 1:224; % индексы точек с меньшим номером по R
nRminus = 33:256;
nnRminus = 225:256; % номера точек, которые есть лишь для Rminus
idxTplus = []; % индексы точек с бОльшим номером по theta
idxTminus = []; % индексы точек с меньшим номером по theta
for i = 0:7
idxTminus = [idxTminus ([32 1:31]+i*32)];
idxTplus = [idxTplus ([2:32 1]+i*32)];
endfor
nTplus = nTminus = 1:256;
% теперь вычисляем расстояния между соседними точками
dRplus = R(nRplus) - R(idxRplus);
dRminus = R(nRminus) - R(idxRminus);
dTplus = theta(nTplus) - theta(idxTplus) ;
do
tmp = find(dTplus < 0);
dTplus(tmp) += 2*pi;
until(isempty(tmp))
dTplus .*= (R(nTplus)+R(idxTplus))/2;
dTminus = theta(nTminus) - theta(idxTminus);
do
tmp = find(dTminus > 0);
dTminus(tmp) -= 2*pi;
until(isempty(tmp))
dTminus .*= (R(nTminus)+R(idxTminus))/2;
Phi = zeros(1, 256);
% прибавляем к фазе средние отклонения
K = N;
while(K-- > 0)
% проходим изнутри наружу, вычисляя точки с бОльшим R
Phi(nRminus) = Phi(idxRminus) + dRminus.*(dR(idxRminus)+dR(nRminus))/2;
endwhile
DPRM = Phi; Phi = zeros(1, 256);
K = N;
while(K-- > 0)
% проходим снаружи вовнутрь, вычисляя точки с меньшим R
Phi(nRplus) = Phi(idxRplus) + dRplus.*(dR(idxRplus)+dR(nRplus))/2;
endwhile
DPRP = Phi; Phi = zeros(1, 256);
K = N;
while(K-- > 0)
% проходимся по угловым проекциям
Phi(nTplus) = Phi(idxTplus) + dTplus.*(dT(idxTplus)+dT(nTplus))/2;
endwhile
DPTP = Phi; Phi = zeros(1, 256);
K = N;
while(K-- > 0)
Phi(nTminus) = Phi(idxTminus) + dTminus.*(dT(idxTminus)+dT(nTminus))/2;
endwhile
Phi += DPRM + DPRP + DPTP;
endfunction
Результаты получаются уж совсем страшными:
"Восстановленный" волновой фронт
"Восстановленный" (синим) и идеальный (красным) волновой фронт
Попробуем сравнить "фронты" на 10 и на 1:
Phi10 = iterate_Phi_Linear(10, R, phi, dR, dphi);
[dPhi10 img] = make_mirror_map(Phi10, dR, R, phi);
imagesc(img)
Phi1 = iterate_Phi_Linear(1, R, phi, dR, dphi);
[dPhi1 img] = make_mirror_map(Phi1, dR, R, phi);
imagesc(img)
С каждой итерацией наблюдается все большее и большее расхождение.
Похоже, этот способ уж совсем никуда не годится.
Из-за того, что наша гартманнограмма имеет не прямоугольную, а концентрическую сетку, возникает уйма проблем, связанных с ее обработкой. Простейшие классические методы дают сильные погрешности из-за краевых эффектов.
Возможно, дифференциальный метод (при помощи полиномов Цернике) будет более надежным.
