Применение Octave для вычисления центра вращения поля

[eddy_em]

Итак, применим Octave для обработки набора снимков звездного поля, сделанного для определения центра вращения поля (при остановленном компенсаторе вращения). Нужно это, чтобы определить, насколько расходится центр вращения поля с центром вращения компенсатора.

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

Чтобы отождествить схожие объекты, выберем на каждом изображении поочередно каждый объект, как опорный, и построим матрицу векторов (в полярной системе координат) от этого объекта до всех остальных. Если (xc, yc) - массив координат объектов на изображении, N - номер объекта в массиве, то построить массив полярных координат векторов, соединающих точку N со всеми остальными, можно так:

function coords = find_tree(xc, yc, N)
%
% coords = find_tree(xc, yc, N)
% для заданных массивов координат пятен, xc и yc, вычисляет вектора от точки
% с номером N до остальных точек в полярной системе координат.
% возвращает матрицу, в которой первая колонка - полярная координата
% вектора, вторая - угловая координата.
%
	if(size(xc, 1)) == 1 % проверяем, являются ли массивы координат столбцами
		xc = xc';
	endif
	if(size(yc, 1)) == 1
		yc = yc';
	endif
	if(size(xc) != size(yc) || N > size(xc, 1)) % проверяем размеры
		coords = [];
	endif
	xn = xc(N); yn = yc(N); % координаты N-й точки
	Dx = xc - xn; Dy = yc - yn; % смещения остальных точек относительно N-й
	R = sqrt(Dx.^2+Dy.^2); % полярная координата
	Phi = atan2(Dy,Dx)*180/pi; % угловая координата
	coords = [R Phi];
endfunction

Чтобы упростить построение массив опорных векторов для всех точек, напишем следующую функцию:

function [Tree xc yc] = get_tree(treshold, i)
%
% 		[Tree xc yc] = get_tree(treshold, i)
% Поиск дерева векторов расстояний между точками и координат центров этих точек
% Вход:
%	treshold - порог интенсивности для бинаризации изображения
%	i - номер кадра в серии
% Выход:
%   Tree - трехмерный массив с расстояниями между точками, у которого
%        * первая координата - номер второй точки в паре
%        * вторая координата - полярные координаты вектора расстояния между точками
%        * третья координата - номер первой точки в паре
%	xc, yc - координаты центров
%
	Tree = [];
	name = sprintf("object_%04d.fit", i);
	II = fits_read(name);
	if(isempty(II)) return; endif
	printf("Find centers of %s ... ", name); fflush(1);
	[xc, yc] = find_centers(II, treshold);
	for j = 1:size(xc,2)
		Tree(:,:,j) = find_tree(xc, yc, j);
	endfor
	printf("done (%d vectors)!\n", size(Tree, 1));
endfunction

Здесь treshold - нижний порог для выявления областей (звезд) на изображении; i - номер изображения.

Теперь если мы возьмем два последовательно снятых изображения, мы можем получить координаты звезд и векторы, соединяющие каждую звезду на изображениях с остальными звездами. Количество обнаруженных звезд может меняться как в силу изменения условий наблюдения (а также из-за стохастических процессов, если интенсивность сигнала от звезды близка к определяемому нами порогу), так и в силу вращения поля. Таким образом, нам необходимо теперь исследовать два набора векторов, чтобы выявить две группы векторов, наиболее полно соответствующих друг другу.

Пробежимся по группам первого набора, перебирая для каждой группы все группы второго набора в поиске наиболее похожей. Для оценки схожести переберем все вектора группы первого набора и найдем наиболее близкие к ним по длине в группе второго набора (с точностью до некоторого dR). Для всех найденных точек вычислим разность угловых координат. Медиана этой величины даст нам значение dPhi - примерного угла поворота второго набора относительно первого. Отсеяв теперь с все разности угловых координат точек, выходящих за пределы dPhi±dphi (где dphi - требуемая точность по угловой координате), посчитаем количество оставшихся точек. Они и будут характеризовать данную группу.

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

function coord_indexes = find_cluster_c(CC, CC1, dR, dPhi)
%
% 		coord_indexes = find_cluster_c(CC, CC1, dR, dPhi)
% Поиск соответствия индексов объектов из CC объектам из CC1
% Вход:
%   CC, CC1 - массивы векторов для двух соседних изображений
%   dR - погрешность по R (если |r1-r0| < dR, считается, что r1 == r0)
%   dPhi - аналогичная погрешность по угловой координате
% Выход:
%   coord_indexes - массив соответствия индексов точек изображений [индекс1 индекс2; ...]
%
	Ncc = size(CC, 3); Ncc1 = size(CC1, 3); % кол-во групп точек в каждом кластере
	Cluster = []; % массив по строкам: кол-во совпадений, индекс1, индекс2, угол поворота
	for jj = 1 : Ncc % цикл по группам первого набора
		Nnear = []; % массив, в который будет складываться кол-во близких точек для каждого кластера
		PhiMed = []; % медианы углов поворота 
		r0 = CC(:,1,jj); phi0 = CC(:,2,jj);
		for ii = 1 : Ncc1 % цикл по группам второго набора
			r1 = CC1(:,1,ii); phi1 = CC1(:,2,ii);
			points = {}; % объединение для поиска соответствия точек (номер точки 1; номера близких точек в 2)
			dphi_s = []; % массив, в который будут заноситься для анализа найденные углы
			for i = 1:size(r0); % цикл по векторам первого набора
				d = abs(r1-r0(i)); % разности между векторами
				idx = find(d < dR); % ищем наиболее близкие по r
				if(isempty(idx)) continue; endif; % близких нет - идем дальше
				points = [ points; {i, idx} ]; % добавляем номера похожих точек
			endfor
			for i = 1:size(points, 1)
				dphi = abs(phi1(points{i,2}) - phi0(points{i,1}));
				dphi_s = [dphi_s; dphi];
			endfor
			if(isempty(dphi_s)) continue; endif;
			dphi = median(dphi_s); % примерный угол поворота набора 2 относительно 1
			n = size(find(abs(dphi_s-dphi) < dPhi), 1);
			Nnear = [ Nnear, n ]; % формируем массив
			PhiMed = [ PhiMed, dphi ];
		endfor
		idx = find(Nnear == max(Nnear))(1); % нашли номер набора, наиболее близкого к jj-му
		Cluster = [ Cluster; Nnear(idx) jj idx PhiMed(idx) ];
	endfor
	% теперь ищем соответствие точек в наборе
	n = max(Cluster(:,1)); % максимальное соответствие
	idx = find(Cluster(:,1) == n)(1); % на всякий случай берем самое первое
	first = Cluster(idx, 2);  % номер набора на изображении 1
	second = Cluster(idx, 3); % -//- 2
	Phi = Cluster(idx, 4);    % угол поворота 2 относительно 1
	r0 = CC(:,1,first); phi0 = CC(:,2,first);
	r1 = CC1(:,1,second); phi1 = CC1(:,2,second);
	coord_indexes = []; 
	for i = 1:size(r0);     % цикл по векторам первого набора
		d = abs(r1-r0(i));  % разности между векторами
		idx = find(d < dR); % ищем наиболее близкие по r
		if(isempty(idx)) continue; endif; % близких нет - идем дальше
		dphi = abs(abs(phi1(idx) - phi0(i)) - Phi);
		if(size(idx, 1) != 1) % несколько точек
			n = find(dphi == min(dphi))(1); % номер наиболее близкой точки
			idx = idx(n);
			dphi = dphi(n);
		endif
		if(dphi > dPhi) continue; endif; % близких нет - идем дальше
		printf("p1: %d <-> p2: %d\n", i, idx);
		coord_indexes = [ coord_indexes; i idx ];
	endfor	
endfunction

(довольно длинная и "велосипедная").

Теперь по имеющимся наборам координат и индексов мы можем вычислить параметры прямых, примерно проходящих через центр вращения поля (точнее - прямых, проходящих через середину отрезка, соединяющего два изображения одного и того же объекта, перпендикулярно нему). Для этого служит функция

function [Y k] = find_line_params(x0,y0, x1,y1)
%
% 		[Y k] = find_line_params(x0,y0, x1,y1)
% Поиск параметров линии y = Y + kx, проходящей через середину отрезка
% (x0,y0) - (x1,y1), перпендикулярно этому отрезку
%
	if(y1 == y0)
		k = inf;
		Y = y0;
	else
		k = -(x1-x0)/(y1-y0);
		if(k == 0)
			Y = y0;
		else
			xx = 0.5*(x0+x1); % середина отрезка
			yy = 0.5*(y0+y1);
			Y = yy - xx * k;
		endif
	endif
endfunction

Ну и наконец, вычислив попарно центры пересечения для всех найденных прямых, получим геометрическое место положения центра вращения (по большому количеству изображений его можно будет уточнить).

Для автоматизации всего процесса служит функция

function [Xc Yc] = cluster__(i0, i1, treshold, dR, dPhi)
%
%		[Xc Yc] = cluster__(i0, i1, treshold, dR, dPhi)
% 
% Вход:
%   i1, in - начальный и конечный номера изображений
%   treshold - порог для поиска пятен
%	dR, dPhi - допуск по полярным координатам (чтобы считать точки близкими)
% Выход:
%	Xc, Yc - центр кадра
%
	n = 0; i_start = i0;
	Xc = []; Yc = [];
	do
		[CC xc yc] = get_tree(treshold, i_start);
		i_start++;
	until(!isempty(CC) || i_start > i1)
	for i = i_start:i1
		[CC1 xc1 yc1] = get_tree(treshold, i);
		if(isempty(CC1)) continue; endif
		printf("pair %d\n", ++n);
		indexes = find_cluster_c(CC, CC1, dR, dPhi);
		sz =  size(indexes,1);
		YY = []; kk = []; % опустошаем массивы с параметрами прямых
		for i = 1:sz
			i1 = indexes(i,1);
			i2 = indexes(i,2);
			[Y k] = find_line_params(xc(i1),yc(i1), xc1(i2),yc1(i2));
			YY = [YY Y]; kk = [kk k]; % добавляем очередные
		endfor
		for j = 1:sz-1
			for i = j+1:sz
				XC = -(YY(i)-YY(j))/(kk(i)-kk(j));
				YC = kk(j)*XC + YY(j);
				Xc = [Xc XC]; Yc = [Yc YC]; % накапливаем центры
			endfor
		endfor
		CC = CC1; xc = xc1; yc = yc1; % передаем параметры для следующей пары
	endfor
endfunction

которая получает аргументами номера первого и последнего анализируемого изображения (i0, i1), порог, допуск на координаты вектора при поиске соответствия объектов второго изображения объектам первого. А возвращает эта функция набор координат центров.

Накопив побольше таких координат, можно тем или иным способом выяснить положение центра вращения поля достаточно точно.

P.S. Я не учитывал здесь то, что в течение наблюдений изображение также может испытывать линейные смещения. Об этом - в следующий раз.