Окончание моих мытарств с центром вращения

[eddy_em]

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

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

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

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

Итак, FITS-файлы считываются уже приводимой мною функцией fits_read:

function [image header] = fits_read(filename)
%
% считываем FITS-файл и переворачиваем для правильного отображения
% Координата y отсчитывается в результате сверху!
% (если не делать flipdim, отображаться будет зеркально)
%
% ЗАВИСИМОСТИ:
%	read_fits_image (из пакета fits)
%
	image = [];
	printf("Read file %s ... ", filename); fflush(1);
	% проверяем, есть ли такой файл и является ли он файлом:
	[ s e m ] = stat(filename);
	if(e != 0 || !S_ISREG(s.mode))
		printf("No such file!\n"); fflush(1);
		return;
	endif
	[ image header ] = read_fits_image(filename);
	if(!isempty(image))
		image = flipdim(image',1); % переворачиваем, чтобы отображалось верно
		printf("OK!\n");
	else
		printf("Bad image!\n");
	endif
	fflush(1);
endfunction

Центры пятен определяются все той же элементарной пороговой бинаризацией:

function [ xc, yc ] = find_centers(Img, treshold)
% [ xc, yc ] = find_centers(Img, treshold)
% Нахождение центров пятен изображения Img с интенсивностью выше treshold
% xc, yc - массивы координат
%
	tmp = zeros(size(Img));
	mm = mean2(Img) + std2(Img)/4; % считаем простенький порог
	tres1 = max(mm, treshold);
	printf("\nTreshold: %g\n", tres1); fflush(stdout);
	tmp(find(Img > tres1)) = 1;
	tmp = medfilt2(logical(tmp));
	[ Labels, Nlabels ] = bwlabel(tmp);
	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

Не изменилась и функция вычисления относительных координат обнаруженных пятен на изображении в полярной системе координат:

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 - координаты центров
%
% ЗАВИСИМОСТИ:
%	fits_read
%	find_tree
%	find_centers
%
	Tree = []; xc = []; yc = [];
	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

Алгоритм обнаружения одних и тех же пятен на двух изображениях тоже не изменился:

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; % близких нет - идем дальше
		coord_indexes = [ coord_indexes; i idx ];
	endfor	
endfunction

кстати, если объектов на изображении очень много, то эта функция выполняется очень долго (это и понятно: таки двойной цикл), поэтому имеет смысл изменить и ее: либо воспользовавшись методом научного тыка Монте-Карло (проверить лишь несколько случайных пар в надежде, что подвыборка полностью репрезентативна), либо вообще сравнив первое и последнее изображение (в надежде, что случайные колебания системы не "выгнали" часть объектов за поле зрения в середине последовательности кадров). При переносе кода на С (с использованием openmp/pthreads, а то и CUDA) нужно будет обязательно оптимизировать и это.

Итак, вот единственные две функции, которые я изменил. Первая - получение массива (строка - условный номер объекта, столбец - номер изображения) с координатами одних и тех же объектов на всех изображениях:

function [X Y] = get_centers(i0, i1, treshold, dR, dPhi)
%
%		[X Y] = get_centers(i0, i1, treshold, dR, dPhi)
% Помещает в массивы X,Y координаты центров пятен для каждой пары
% изображений
% Вход:
%   i1, in - начальный и конечный номера изображений
%   treshold - порог для поиска пятен
%	dR, dPhi - допуск по полярным координатам (чтобы считать точки близкими)
% Выход:
%	X, Y - матрицы набора центров для пятен
%		номер строки соответствует номеру объекта для первого изображения
%		номер столбца - номер изображения
%		координаты в каждой строке - для одного и того же объета
%
% ЗАВИСИМОСТИ:
%	get_tree -> fits_read, find_tree, find_centers
%	find_cluster_c
%
	n = 1; % счетчик изображений + координата текущего столбца в X,Y
	i_start = i0; % номер первого изображения для сравнения
	CS = [];
	printf("\nimage 1\n", i);
	do
		[CC xc yc] = get_tree(treshold, i_start);
		im1 = i_start;
		i_start++;
	until(!isempty(CC) || i_start > i1)
	X = xc'; Y = yc'; % первый столбец - первое изображение
	Counters = zeros(size(X)); % и сразу заполняем нулями массив счетчиков
	for i = i_start:i1
		[CC1 xc1 yc1] = get_tree(treshold, i);
		if(isempty(CC1)) continue; endif
		printf("\nimage %d, ", ++n);
		% все сравнения идут с первым изображением
		indexes = find_cluster_c(CC, CC1, dR, dPhi);
		printf("%d points\n\n", size(indexes, 1));
		% Заполняем координаты:
		X(indexes(:,1),n) = xc1(indexes(:,2));
		Y(indexes(:,1),n) = yc1(indexes(:,2));
		% и увеличиваем счетчик для найденных пятен
		Counters(indexes(:,1))++;
	endfor
	% стираем объекты, не присутствующие на всех изображениях
	idx = find(Counters != n-1);
	X(idx,:) = []; Y(idx,:) = [];
endfunction

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

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

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

function [Xc Yc] = matr_center(X, Y)
%
%		[Xc Yc] = matr_center(X, Y)
% Вычисляет координаты центра кадра, исходя из матричных операций
% Вход:
%   X, Y - массивы с координатами пятен (из функции get_centers)
% Выход:
%	Xc, Yc - массив центров кадров
%
	Xc = []; Yc = [];
	imax = size(X,2);
	for i = 1:imax-1
		for j = i+1:imax
			% выделяем координаты объектов
			X0 = X(:,i)';
			X1 = X(:,j)';
			Y0 = Y(:,i)';
			Y1 = Y(:,j)';
			% получаем примерный угол поворота изображения (вектора - от центра тяжести)
			vec0X = X0 - mean(X0); vec0Y = Y0 - mean(Y0);
			vec1X = X1 - mean(X1); vec1Y = Y1 - mean(Y1);
			% Следующий код можно раскомментировать, если используется выборка из малого количества изображений
			%phi0 = atan2(vec0Y, vec0X); phi1 = atan2(vec1Y, vec1X);
			%% прибавляем к отрицательным 2пи, чтобы избежать переходов -pi+a -> pi-b
			%phi0(find(phi0<0)) += 2*pi; phi1(find(phi1<0)) += 2*pi;
			%dphi = median(phi1 - phi0); % медианный угол поворота изображения
			%R = [cos(dphi) -sin(dphi); sin(dphi) cos(dphi)]; % матрица поворота
			%RR1 =  R * [vec0X;vec0Y]; % вектора, которые получились бы при повороте системы 0
			%dR = sqrt((vec1X-RR1(1,:)).^2+(vec1Y-RR1(2,:)).^2); % расхождения между реальными и теоретическими
			%idx = find(dR < median(dR)+std(dR)); % отсекаем явные промахи
			%if(size(idx,2) < 3) % слишком много промахов
				%printf("Image %d: too much bad data\n", i);
				%continue;
			%endif
			%X0 = X0(idx); X1 = X1(idx);Y0 = Y0(idx); Y1 = Y1(idx);
			% получаем матрицу геометрического преобразования
			A = [X1;Y1;ones(size(X1))]/[X0;Y0;ones(size(X1))];
			% вычисляем примерные координаты центра
			CRDS = (eye(2)-A(1:2,1:2)) \ A(1:2,3);
			Xc = [Xc CRDS(1)]; Yc = [Yc CRDS(2)]; % накапливаем центры
		endfor
	endfor
	Xmed = median(Xc)
	Ymed = median(Yc)
endfunction

(закомментированный участок кода нужно раскомментировать, если изображений маловато для того, чтобы набрать приличную статистику). Функция вычисляет центр вращения исходя из матричного способа, приведенного мной в одной из предыдущих публикаций. Принцип прост: мы вычисляем методом наименьших квадратов матрицу преобразования, A, при воздействии которой на вектор r (координаты объекта на первом изображении из пары) получается вектор r' (координаты объекта на втором изображении). Ну, а координаты центра вычисляются на основе предположения, что изображение лишь вращалось, но не смещалось (если предполагать, что оно и вращалось, и смещалось, задача получится некорректной для одной пары изображений и решить ее будет невозможно), исходя из формулы для отвечающих за смещение членов матрицы A: (I-R)c=d. Здесь I - единичная матрица, R - матрица вращения (левый угол 2х2 матрицы A), c - искомый радиус-вектор центра вращения, d - отвечающие за смещение компоненты матрицы A (правый столбец, над единицей).

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

Вычисленные центры вращения

Проблема выбросов, о которой я говорил в начале, действительно возникала из-за попытки вычислений, когда смещение объектов в результате вращения поля было сопоставимо с ошибкой вычисления их положения. Об этом скажу дальше.

Решив-таки погуглить, а не дурак ли я, я нашел лишь еще один вариант в статье McKein, Abbott & King. Они предложили практически то же самое, что я уже делал, но вместо вычисления матрицы A определять угол поворота всей системы, затем вычислять R и на основе выражения для преобразования вращения вокруг произвольной точки как вращения вокруг начала координат с последующим сдвигом (и что я, дурак, об этом сразу не догаладся?), вычислять вектор v (это самое смещение), который и использовать в приведенной выше формуле вместо d.

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

Итак, сама функция:

function [Xc Yc] = matr_center_notmine(X, Y)
%
%		[Xc Yc] = matr_center(X, Y)
% Вычисляет координаты центра кадра, исходя из матричных операций
% Вход:
%   X, Y - массивы с координатами пятен (из функции get_centers)
% Выход:
%	Xc, Yc - массив центров кадров
%
	Xc = []; Yc = [];
	imax = size(X,2);
	for i = 1:imax-1
		for j = i+1:imax
			% выделяем координаты объектов
			X0 = X(:,i)';
			X1 = X(:,j)';
			Y0 = Y(:,i)';
			Y1 = Y(:,j)';
			% получаем примерный угол поворота изображения (вектора - от центра тяжести)
			vec0X = X0 - mean(X0); vec0Y = Y0 - mean(Y0);
			vec1X = X1 - mean(X1); vec1Y = Y1 - mean(Y1);
			phi0 = atan2(vec0Y, vec0X); phi1 = atan2(vec1Y, vec1X);
			% прибавляем к отрицательным 2пи, чтобы избежать переходов -pi+a -> pi-b
			phi0(find(phi0<0)) += 2*pi; phi1(find(phi1<0)) += 2*pi;
			dphi = median(phi1 - phi0); % медианный угол поворота изображения
			R = [cos(dphi) -sin(dphi); sin(dphi) cos(dphi)]; % матрица поворота
			% Получаем вектор v
			v = median(R * [X0; Y0] - [X1; Y1], 2);
			% вычисляем примерные координаты центра
			CRDS = (eye(2)-R) \ v;
			%printf("Center: (%g, %g)\n", CRDS(1), CRDS(2));
			Xc = [Xc -CRDS(1)]; Yc = [Yc -CRDS(2)]; % накапливаем центры
		endfor
	endfor
	Xmed = median(Xc)
	Ymed = median(Yc)
endfunction

На удивление, результаты получились практически идентичными, однако время выполнения получилось в мою пользу:

octave:828> tic;[Xnotmine Ynotmine] = matr_center_notmine(X,Y);toc
Xmed =  1510.9
Ymed =  1370.2
Elapsed time is 1e+01 seconds.

octave:829> tic; [Xmine Ymine] = matr_center(X,Y); toc
Xmed =  1511.0
Ymed =  1370.2
Elapsed time is 8 seconds.

Я могу объяснить эту странную разность во времени тем, что в Octave криво реализован алгоритм вычисления медианы (по-видимому, посредством быстрой сортировки, а не быстрого выбора), зато быстро реализован механизм вычисления наименьшими квадратами.

Ну и наконец приведу обещанные графики погрешностей:

   
Отклонение (в абсолютных единицах) вычисленного центра от медианного значения (справа - увеличенно).

и корреляции ошибок:

Корреляция ошибок: почти идеальна.

В общем, судя по результатам, пользоваться можно одним из двух методов в зависимости от конечной реализации: в Octave шустрее работает мой способ, а на C наверняка более выгодной будет способ McCane сотоварищи.