Восстановление волнового фронта по гартманнограмме методом Саусвелла

[eddy_em]

Наконец-то нашел я время вернуться к обработке гартманнограмм.
Гартманнограммы я сгенерировал при помощи уже описанной модели, из нее получил значения нормалей к волновому фронту, преобразовал их в полярные координаты и далее стал пытаться восстановить волновой фронт так, чтобы получить нечто похожее на оригинальный.
Оригинальный волновой фронт имеет следующий вид:

Сначала — вкратце расскажу о простейших (и → грубейших) методах восстановления волнового фронта.

Линейное интегрирование

Это — простейший (и наименее точный) метод восстановления волнового фронта. В некоторой точке, являющейся начальной (например, на краю волнового фронта), величина волнового фронта считается равной нулю. Далее осуществляется переход вдоль некоей линии, соединяющей изображения пятен на гартманнограмме, к следующей точке. В ней волновой фронт вычисляется как линейная экстраполяция наклона фронта относительно предыдущей точки по значению наклона в этой точке:

где s - проекция расстояния между точками Ф_i и Ф_i-1 на ось ζ (которая может быть и нелинейной — например, профиль вдоль окружности). Компоненты Ф^ζ вычисляются для двух ортогональных осей α и β, а затем суммированием получается линейное приближение к волновому фронту:

Ф_i = Ф_i^α + Ф_i^β.

Восстановление Саусвелла

Саусвелл (W. Southwell) предложил простой метод зонального восстановления волнового фронта. Метод заключается в том, что волновой фронт в каждой точке вычисляется на основе данных о наклоне фронта во всех соседних точках:

где Ф_n,m — величина волнового фронта в пятне (n,m);
Δ_i,j — расстояние между пятнами (n,m) и (n+i,m+j),
I_{n+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)

С каждой итерацией наблюдается все большее и большее расхождение.
Похоже, этот способ уж совсем никуда не годится.

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

Comments

[vlkamov.livejournal.com]

А тупо решить систему линейных уравнений ?

Правда их там триллион, но сейчас даже в телефонах гигафлопы.