Вычисление нормалей к волновому фронту (для модели гартманограмм)
Продолжу цикл своих публикаций относительно мытарств с обработкой гартманнограмм.
Ранее я уже описывал, как получить модели пред- и зафокальных гартманнограмм. Теперь пора перейти к следующему этапу - определению координат центров пятен и получению значений нормалей к проекции волнового фронта, например, на предфокальное изображение (оно не зеркалировано, так что его использовать удобнее для расчета отклонений формы зеркала от идеальной).
Пока что особые точности нам не нужны (мы просто хотим проверить несколько моделей восстановления формы поверхности волнового фронта по значениям нормалей и выбрать наиболее точный и простой - а потом уже будем повышать точность), будем вычислять центры пятен как обычные центры тяжести. Для этого построим бинарную маску, чтобы пронумеровать все пятна, а затем уже вычислим центр тяжести каждого пятна. Чтобы маска была чуть больше наших пятен, отфильтруем изображение гауссовым фильтром 9х9 со среднеквадратичным отклонением, равным трем. Бинаризацию будем проводить по уровню интенсивности больше 0.5.
Вот такая простая функция получается для вычисления простых центров тяжести наших пятен:
function [ xc, yc ] = find_centers(Img)
% [ xc, yc ] = find_centers(Img)
% Нахождение центров пятен бинарного изображения Img
% xc, yc - массивы координат
%
f = fspecial('gaussian', 9, 3);
tmp = imfilter(Img, f);
[ Labels, Nlabels ] = bwlabel(tmp > 0.5);
if(Nlabels < 1)
printf(stderr, "\nError! There's no spots!!!\n\n");
return
endif
[yy, xx] = ndgrid(1:size(Img,1),1:size(Img,2));
xc = []; yc = [];
for i = [1 : Nlabels]
idxs = find(Labels == i);
Is = sum(Img(idxs));
xc = [ xc sum(Img(idxs) .* xx(idxs)) / Is ];
yc = [ yc sum(Img(idxs) .* yy(idxs)) / Is ];
end
endfunction
Полученные координаты - просто два массива, в которых точки идут по порядку сканирования функцией bwlabel - т.е. начиная с начала координат и перемещаясь слева направо и снизу вверх. Нам же необходимо сопоставить точки зафокального изображения соответствующим точкам предфокального. В своей fitsview для решения этой задачи я вынужден был не только пронумеровать точки, но и привести их в единую систему координат (т.е. компенсировать смещение изображения, возникающее из-за неточности наведения на объект, а также вращение поля, возникающее из-за неточности установки P2). Здесь же (в модели) изображения отцентрированы, вращение отсутствует - т.е. процесс нумерации сильно упрощается.
Для того, чтобы получить значение нормалей и их координаты (в метрах), нам необходимо знать координаты (в пикселях) пятен на пред- и зафокальной гартманнограммах, расстояние между гартманнограммами и размер пикселя на изображении. На выходе будем иметь проекции x и y нормалей (проекция z нам не нужна, т.к. x^2+y^2+z^2 = 1 - это же нормаль), соответствующие координаты и (на всякий случай - для тестов) соответствие номеров точек из входных массивов данных.
Основная функция поиска нормалей проста: мы получаем индексы одних и тех же пятен на пред- и зафокальном изображении, вычисляем на основе этих индексов смещения лучей и нормализуем эти смещения, получая в результате искомые нормали.
function [Nx Ny X Y Ipre Ipost] = compute_normals(Xpre, Ypre, Xpost, Ypost, dZ, CCDpixsize)
%
% [Nx Ny X Y Idx] = compute_normals(Xpre, Ypre, Xpost, Ypost, dZ, CCDpixsize)
% вычисляет нормали к волновому фронту в точках
% с координатами Xpre(Idx), Ypre(Idx)
% выход:
% Nx, Ny - проекции нормалей (Nx^2+Ny^2+Nz^2 = 1)
% X,Y - координаты точек (в метрах) на изображении pre
% Ipre, Ipost - индексы (по порядковому номеру) точек в pre и post
% вход:
% Xpre, Ypre - координаты центров пятен на предфокальном снимке (для него
% и вычисляются нормали)
% Xpost, Ypost - координаты на зафокальном снимке
% dZ - расстояние между снимками
% CCDpixsize - размер пикселя в метрах
%
% Точки pre и post могут быть в любом порядке (функция находит соответствие
% между точками)
%
if(nargin < 6) CCDpixsize = 12.5e-6; endif
Npts = max(size(Xpre)); % количество точек
if(Npts != max(size(Xpost)))
fprintf(1, "Error: different number of points in pre- and postfocal images!\n");
return;
endif
Ipre = findIdx(Xpre, Ypre, 1);
Ipost = findIdx(Xpost, Ypost, 0);
Nx = CCDpixsize*(Xpost(Ipost) - Xpre(Ipre)); % смещение луча по X
Ny = CCDpixsize*(Ypost(Ipost) - Ypre(Ipre)); % смещение луча по Y
Nl = sqrt(Nx.^2 + Ny.^2 + dZ^2);% длина луча
Nx ./= Nl; % нормализуем лучи
Ny ./= Nl;
X = Xpre(Ipre) * CCDpixsize;
Y = Ypre(Ipre) * CCDpixsize;
endfunction
Для вычисления индексов, сначала найдем центр тяжести изображения и переведем координаты точек в радиальные относительно этого центра. Затем отсортируем их по угловой координате и, проходя точку за точкой, будем искать перепад угловой координаты больше некоторого порогового значения. За это пороговое значение примем 0.05 радиан (это примерно четверть от нормального расстояния между лучами гартманнограммы, составляющего 11.25°), с выбором порога следует быть осторожным, т.к. слишком малый порог вызовет ложные срабатывания (особенно если качество поверхности зеркала плохое), а слишком большой может вообще не дать срабатываний (по той же причине).
Всем точкам каждого луча присваивается номер этого луча (для этого мы вводим дополнительный вектор данных, в котором каждой точке из входного массива данных присваивается ее номер). Нулевым лучем считаем первый против часовой стрелки луч, следующий после ближнего к центру маркера.
Если у нас происходит два срабатывания подряд, это означает, что мы нашли маркер (кстати, поэтому порог должен быть меньше половины расстояния между лучами - иначе можно "прозевать" маркер). Маркер нумеруется по номеру предыдущего луча.
Далее ближайший к центру маркер мы нумеруем цифрой 800, а дальний - 900 (чтобы они не мешались при общей нумерации остальных точек).
Затем, в соответствии с найденным нулевым маркером, мы перенумеровываем все наши лучи так, чтобы нумерация следовала по часовой стрелке и нулевым лучом был луч, располагающийся против часовой стрелки относительно нулевого маркера.
После такой перенумерации нам остается лишь отсортировать точки по радиальной координате и аналогичным образом, проверяя R на скачки, большие некоего порога (будем считать его равным 1/7 среднего расстояния между кольцами), нумеруем кольца. В итоге мы формируем для каждой точки (кроме уже помеченных маркеров) уникальный идентификатор - номер, состоящий из номера радиуса (десятки и единицы) и номера кольца (сотни): т.е. номер кольца (начиная от нулевого - внутреннего) мы умножаем на сто и прибавляем к нему номер радиуса. Если теперь мы отсортируем эти идентификаторы, индексы идентификаторов из неотсортированного массива и дадут порядок следования пятен, который можно использовать в дальнейшем для сопоставления одинаковых пятен на за- и предфокальном снимках друг другу.
function I = findIdx(X, Y, preFlag)
% вычисляет индексы для пятен с координатами X,Y
% preFlag == 1, если изображение предфокальное, иначе == 0
%
PHI_STEP = 0.05; % порог для выделения радиусов
%aStep = pi / 16; % 11.25градусов - угол между радиусами
N_ANGLES = 16; % количество диаметров на гартманнограмме
N_RAY_MAX = 2*N_ANGLES-1;
N_RAYS = 2*N_ANGLES;
N_RINGS = 8; % количество колец на гартманнограмме
Npts = max(size(X));
I = zeros(1, Npts); % номера точек
[Rp, Tp] = conv_to_radial(X, Y); % преобразуем в полярную СК
% СОРТИРУЕМ ПО УГЛУ
[T idx] = sort(Tp);
T0 = T(1);
step = 0; % флаг свежевыполненного перехода на угол > PHI_STEP
Num = 0;
MarkerIdx = [];
for i = 1:Npts;
dT = T(i) - T0; T0 = T(i);
if(dT > PHI_STEP)
Num++;
if(!step)
step = 1;
else
Num--;
printf("Found marker %d\n", Num - 1);
MarkerIdx = [MarkerIdx idx(i-1)];
endif
else
step = 0;
endif
I(idx(i)) = Num;
endfor
% индекс ближнего к центру маркера
m0 = MarkerIdx(1); m1 = MarkerIdx(2);
if (Rp(m0) < Rp(m1)) zMarker = m0; else zMarker = m1; endif;
mZero = I(zMarker);
if(Rp(m0) < Rp(m1))
I(m0) = 800; I(m1) = 900;
else
I(m1) = 900; I(m0) = 800;
endif
if(mZero > N_RAY_MAX) mZero -= N_RAYS
endif;
% перенумеровываем в соответствии с номером mZero, пропуская маркеры
idx = find(I <= N_RAY_MAX);
sid = mZero - I;
sid(find(sid < 0)) += N_RAYS;
I(idx) = sid(idx);
% так как работаем с моделью, угол наклона гартманограммы считаем равным нулю
% и шаг определения наклона пропускаем
% СОРТИРУЕМ ПО РАДИУСУ
[R idx] = sort(Rp);
R0 = R(1);
% погрешность для выявления колец - 1/7 среднего растояния между кольцами
Rtres = (R(end) - R0) / (N_RINGS - 1) / 7.;
Num = 0;
for i = 1:Npts;
if(I(idx(i)) > N_RAY_MAX) continue; endif % маркеры мы уже пронумеровали
dR = R(i) - R0; R0 = R(i);
if(abs(dR) > Rtres) Num += 100; endif; % следующее кольцо - увеличиваем idx
I(idx(i)) += Num;
endfor
[tmp I] = sort(I); % сортируем, чтобы I(N) указывал на N-ю точку
endfunction
Кроме указанных функций, используются еще две вспомогательных.
Получение центра тяжести набора точек:
function [Xc Yc] = get_CG(X, Y)
% вычисление центра тяжести
N = max(size(X));
Xc = sum(X)/N; Yc = sum(Y)/N;
endfunction
и преобразование прямоугольных координат в радиальные:
function [R theta] = conv_to_radial(X,Y)
% преобразует координаты точек X,Y в радиальные относительно центра тяжести
[Xc Yc] = get_CG(X,Y);
yy = Y - Yc; xx = X - Xc;
theta = atan2(yy, xx);
R = sqrt(xx.^2+yy.^2);
endfunction
Теперь проверим, правильно ли мы нашли соответствие между точками:
[Nx Ny X Y Ipre Ipost] = compute_normals(Xpre, Ypre, Xpost, Ypost, 24.148-23.9);
hold on; for i = [1:16:255]; plot3([Xpre(Ipre(i)) Xpost(Ipost(i))], [Ypre(Ipre(i)) Ypost(Ipost(i))], [0 100]); endfor; hold off
Да, действительно, изображения выстроены "в ряд":
Смотрим векторное поле (командой quiver(X, Y, Nx, Ny,3)), вектора должны быть более-менее одинаково ориентированы к центру. Все в порядке:
Продолжение следует (теперь необходимо взять хотя бы 3-4 метода восстановления волнового фронта и посмотреть, что получится).
Ранее я уже описывал, как получить модели пред- и зафокальных гартманнограмм. Теперь пора перейти к следующему этапу - определению координат центров пятен и получению значений нормалей к проекции волнового фронта, например, на предфокальное изображение (оно не зеркалировано, так что его использовать удобнее для расчета отклонений формы зеркала от идеальной).
Пока что особые точности нам не нужны (мы просто хотим проверить несколько моделей восстановления формы поверхности волнового фронта по значениям нормалей и выбрать наиболее точный и простой - а потом уже будем повышать точность), будем вычислять центры пятен как обычные центры тяжести. Для этого построим бинарную маску, чтобы пронумеровать все пятна, а затем уже вычислим центр тяжести каждого пятна. Чтобы маска была чуть больше наших пятен, отфильтруем изображение гауссовым фильтром 9х9 со среднеквадратичным отклонением, равным трем. Бинаризацию будем проводить по уровню интенсивности больше 0.5.
Вот такая простая функция получается для вычисления простых центров тяжести наших пятен:
function [ xc, yc ] = find_centers(Img)
% [ xc, yc ] = find_centers(Img)
% Нахождение центров пятен бинарного изображения Img
% xc, yc - массивы координат
%
f = fspecial('gaussian', 9, 3);
tmp = imfilter(Img, f);
[ Labels, Nlabels ] = bwlabel(tmp > 0.5);
if(Nlabels < 1)
printf(stderr, "\nError! There's no spots!!!\n\n");
return
endif
[yy, xx] = ndgrid(1:size(Img,1),1:size(Img,2));
xc = []; yc = [];
for i = [1 : Nlabels]
idxs = find(Labels == i);
Is = sum(Img(idxs));
xc = [ xc sum(Img(idxs) .* xx(idxs)) / Is ];
yc = [ yc sum(Img(idxs) .* yy(idxs)) / Is ];
end
endfunction
Полученные координаты - просто два массива, в которых точки идут по порядку сканирования функцией bwlabel - т.е. начиная с начала координат и перемещаясь слева направо и снизу вверх. Нам же необходимо сопоставить точки зафокального изображения соответствующим точкам предфокального. В своей fitsview для решения этой задачи я вынужден был не только пронумеровать точки, но и привести их в единую систему координат (т.е. компенсировать смещение изображения, возникающее из-за неточности наведения на объект, а также вращение поля, возникающее из-за неточности установки P2). Здесь же (в модели) изображения отцентрированы, вращение отсутствует - т.е. процесс нумерации сильно упрощается.
Для того, чтобы получить значение нормалей и их координаты (в метрах), нам необходимо знать координаты (в пикселях) пятен на пред- и зафокальной гартманнограммах, расстояние между гартманнограммами и размер пикселя на изображении. На выходе будем иметь проекции x и y нормалей (проекция z нам не нужна, т.к. x^2+y^2+z^2 = 1 - это же нормаль), соответствующие координаты и (на всякий случай - для тестов) соответствие номеров точек из входных массивов данных.
Основная функция поиска нормалей проста: мы получаем индексы одних и тех же пятен на пред- и зафокальном изображении, вычисляем на основе этих индексов смещения лучей и нормализуем эти смещения, получая в результате искомые нормали.
function [Nx Ny X Y Ipre Ipost] = compute_normals(Xpre, Ypre, Xpost, Ypost, dZ, CCDpixsize)
%
% [Nx Ny X Y Idx] = compute_normals(Xpre, Ypre, Xpost, Ypost, dZ, CCDpixsize)
% вычисляет нормали к волновому фронту в точках
% с координатами Xpre(Idx), Ypre(Idx)
% выход:
% Nx, Ny - проекции нормалей (Nx^2+Ny^2+Nz^2 = 1)
% X,Y - координаты точек (в метрах) на изображении pre
% Ipre, Ipost - индексы (по порядковому номеру) точек в pre и post
% вход:
% Xpre, Ypre - координаты центров пятен на предфокальном снимке (для него
% и вычисляются нормали)
% Xpost, Ypost - координаты на зафокальном снимке
% dZ - расстояние между снимками
% CCDpixsize - размер пикселя в метрах
%
% Точки pre и post могут быть в любом порядке (функция находит соответствие
% между точками)
%
if(nargin < 6) CCDpixsize = 12.5e-6; endif
Npts = max(size(Xpre)); % количество точек
if(Npts != max(size(Xpost)))
fprintf(1, "Error: different number of points in pre- and postfocal images!\n");
return;
endif
Ipre = findIdx(Xpre, Ypre, 1);
Ipost = findIdx(Xpost, Ypost, 0);
Nx = CCDpixsize*(Xpost(Ipost) - Xpre(Ipre)); % смещение луча по X
Ny = CCDpixsize*(Ypost(Ipost) - Ypre(Ipre)); % смещение луча по Y
Nl = sqrt(Nx.^2 + Ny.^2 + dZ^2);% длина луча
Nx ./= Nl; % нормализуем лучи
Ny ./= Nl;
X = Xpre(Ipre) * CCDpixsize;
Y = Ypre(Ipre) * CCDpixsize;
endfunction
Для вычисления индексов, сначала найдем центр тяжести изображения и переведем координаты точек в радиальные относительно этого центра. Затем отсортируем их по угловой координате и, проходя точку за точкой, будем искать перепад угловой координаты больше некоторого порогового значения. За это пороговое значение примем 0.05 радиан (это примерно четверть от нормального расстояния между лучами гартманнограммы, составляющего 11.25°), с выбором порога следует быть осторожным, т.к. слишком малый порог вызовет ложные срабатывания (особенно если качество поверхности зеркала плохое), а слишком большой может вообще не дать срабатываний (по той же причине).
Всем точкам каждого луча присваивается номер этого луча (для этого мы вводим дополнительный вектор данных, в котором каждой точке из входного массива данных присваивается ее номер). Нулевым лучем считаем первый против часовой стрелки луч, следующий после ближнего к центру маркера.
Если у нас происходит два срабатывания подряд, это означает, что мы нашли маркер (кстати, поэтому порог должен быть меньше половины расстояния между лучами - иначе можно "прозевать" маркер). Маркер нумеруется по номеру предыдущего луча.
Далее ближайший к центру маркер мы нумеруем цифрой 800, а дальний - 900 (чтобы они не мешались при общей нумерации остальных точек).
Затем, в соответствии с найденным нулевым маркером, мы перенумеровываем все наши лучи так, чтобы нумерация следовала по часовой стрелке и нулевым лучом был луч, располагающийся против часовой стрелки относительно нулевого маркера.
После такой перенумерации нам остается лишь отсортировать точки по радиальной координате и аналогичным образом, проверяя R на скачки, большие некоего порога (будем считать его равным 1/7 среднего расстояния между кольцами), нумеруем кольца. В итоге мы формируем для каждой точки (кроме уже помеченных маркеров) уникальный идентификатор - номер, состоящий из номера радиуса (десятки и единицы) и номера кольца (сотни): т.е. номер кольца (начиная от нулевого - внутреннего) мы умножаем на сто и прибавляем к нему номер радиуса. Если теперь мы отсортируем эти идентификаторы, индексы идентификаторов из неотсортированного массива и дадут порядок следования пятен, который можно использовать в дальнейшем для сопоставления одинаковых пятен на за- и предфокальном снимках друг другу.
function I = findIdx(X, Y, preFlag)
% вычисляет индексы для пятен с координатами X,Y
% preFlag == 1, если изображение предфокальное, иначе == 0
%
PHI_STEP = 0.05; % порог для выделения радиусов
%aStep = pi / 16; % 11.25градусов - угол между радиусами
N_ANGLES = 16; % количество диаметров на гартманнограмме
N_RAY_MAX = 2*N_ANGLES-1;
N_RAYS = 2*N_ANGLES;
N_RINGS = 8; % количество колец на гартманнограмме
Npts = max(size(X));
I = zeros(1, Npts); % номера точек
[Rp, Tp] = conv_to_radial(X, Y); % преобразуем в полярную СК
% СОРТИРУЕМ ПО УГЛУ
[T idx] = sort(Tp);
T0 = T(1);
step = 0; % флаг свежевыполненного перехода на угол > PHI_STEP
Num = 0;
MarkerIdx = [];
for i = 1:Npts;
dT = T(i) - T0; T0 = T(i);
if(dT > PHI_STEP)
Num++;
if(!step)
step = 1;
else
Num--;
printf("Found marker %d\n", Num - 1);
MarkerIdx = [MarkerIdx idx(i-1)];
endif
else
step = 0;
endif
I(idx(i)) = Num;
endfor
% индекс ближнего к центру маркера
m0 = MarkerIdx(1); m1 = MarkerIdx(2);
if (Rp(m0) < Rp(m1)) zMarker = m0; else zMarker = m1; endif;
mZero = I(zMarker);
if(Rp(m0) < Rp(m1))
I(m0) = 800; I(m1) = 900;
else
I(m1) = 900; I(m0) = 800;
endif
if(mZero > N_RAY_MAX) mZero -= N_RAYS
endif;
% перенумеровываем в соответствии с номером mZero, пропуская маркеры
idx = find(I <= N_RAY_MAX);
sid = mZero - I;
sid(find(sid < 0)) += N_RAYS;
I(idx) = sid(idx);
% так как работаем с моделью, угол наклона гартманограммы считаем равным нулю
% и шаг определения наклона пропускаем
% СОРТИРУЕМ ПО РАДИУСУ
[R idx] = sort(Rp);
R0 = R(1);
% погрешность для выявления колец - 1/7 среднего растояния между кольцами
Rtres = (R(end) - R0) / (N_RINGS - 1) / 7.;
Num = 0;
for i = 1:Npts;
if(I(idx(i)) > N_RAY_MAX) continue; endif % маркеры мы уже пронумеровали
dR = R(i) - R0; R0 = R(i);
if(abs(dR) > Rtres) Num += 100; endif; % следующее кольцо - увеличиваем idx
I(idx(i)) += Num;
endfor
[tmp I] = sort(I); % сортируем, чтобы I(N) указывал на N-ю точку
endfunction
Кроме указанных функций, используются еще две вспомогательных.
Получение центра тяжести набора точек:
function [Xc Yc] = get_CG(X, Y)
% вычисление центра тяжести
N = max(size(X));
Xc = sum(X)/N; Yc = sum(Y)/N;
endfunction
и преобразование прямоугольных координат в радиальные:
function [R theta] = conv_to_radial(X,Y)
% преобразует координаты точек X,Y в радиальные относительно центра тяжести
[Xc Yc] = get_CG(X,Y);
yy = Y - Yc; xx = X - Xc;
theta = atan2(yy, xx);
R = sqrt(xx.^2+yy.^2);
endfunction
Теперь проверим, правильно ли мы нашли соответствие между точками:
[Nx Ny X Y Ipre Ipost] = compute_normals(Xpre, Ypre, Xpost, Ypost, 24.148-23.9);
hold on; for i = [1:16:255]; plot3([Xpre(Ipre(i)) Xpost(Ipost(i))], [Ypre(Ipre(i)) Ypost(Ipost(i))], [0 100]); endfor; hold off
Да, действительно, изображения выстроены "в ряд":
Смотрим векторное поле (командой quiver(X, Y, Nx, Ny,3)), вектора должны быть более-менее одинаково ориентированы к центру. Все в порядке:
Продолжение следует (теперь необходимо взять хотя бы 3-4 метода восстановления волнового фронта и посмотреть, что получится).