Итак, несмотря на фильтрацию слишком уж близких линий, от выбросов это не спасло. Вот как выглядит распределение вычисленных центров:
(кружочки - вычисленное "новым" методом, крестики - "старым").Вот увеличенная область вблизи медиан:
(даже почти совпали!)Ну и напоследок файлы. Определение параметров линии теперь делается так:
function [Y k flag] = find_line_params(x0,y0, x1,y1)
%
% [Y k] = find_line_params(x0,y0, x1,y1)
% Поиск параметров линии y = Y + kx, проходящей через середину отрезка
% (x0,y0) - (x1,y1), перпендикулярно этому отрезку
% flag == 1, если линия ближе к оси Y, иначе flag == 0
%
% определить угол наклона: если линия ближе к горизонту,
% использовать что есть, иначе - установить флаг и использовать
% отсчеты относительно x
%
flag = (abs(x1-x0) > abs(y1-y0));
if(flag)
[Y k] = getK(x0,y0, x1,y1);
else
[Y k] = getK(y0,x0, y1,x1); % меняем x и y местами
endif
endfunction
function [Y k] = getK(x0,y0, x1,y1)
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
endfunction
Здесь тангенс угла наклона вычисляется относительно наиболее близкой координатной линии. Введенная переменная flag сигнализирует, к какой оси линия ближе.
Ну и сама функция определения координат:
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;
im1 = 0;
Xc = []; Yc = [];
do
[CC xc yc] = get_tree(treshold, i_start);
im1 = 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: images %d and %d\n", ++n, im1, i);
im1 = i;
indexes = find_cluster_c(CC, CC1, dR, dPhi);
sz = size(indexes,1);
YY = []; kk = []; flag = []; % опустошаем массивы с параметрами прямых
for ii = 1:sz
i1 = indexes(ii,1);
i2 = indexes(ii,2);
[Y k f] = find_line_params(xc(i1),yc(i1), xc1(i2),yc1(i2));
YY = [YY Y]; % добавляем очередные
kk = [kk k];
flag = [flag f];
endfor
for j = 1:sz-1
for k = j+1:sz
Y1 = YY(j); Y2 = YY(k);
k1 = kk(j); k2 = kk(k);
switch(bitshift(flag(j), 1) + flag(k))
case 0 % обе линии относительно X
[YC XC ] = get_cross1(Y1, k1, Y2, k2);
case 1 % k-я относительно Y, j-я - относительно X
[XC YC ] = get_cross2(Y1, k1, Y2, k2);
case 2 % j-я относительно Y, k-я - относительно X
[XC YC ] = get_cross2(Y2, k2, Y1, k1);
case 3 % обе линии относительно Y
[XC YC ] = get_cross1(Y1, k1, Y2, k2);
endswitch
Xc = [Xc XC]; Yc = [Yc YC]; % накапливаем центры
endfor
endfor
CC = CC1; xc = xc1; yc = yc1;% передаем параметры для следующей пары
endfor
endfunction
function [XC YC] = get_cross1(Y1, k1, Y2, k2)
% вычисляет пересечение линий с параметрами Y1, Y2 для flag == 0,3
XC = []; YC = [];
if(abs(k1-k2) < 0.1) return; endif
XC = -(Y2-Y1)/(k2-k1);
YC = k1*XC + Y1;
endfunction
function [XC YC] = get_cross2(Y1, k1, Y2, k2)
% вычисляет пересечение линий, у которых Y1 и k1 заданы относительно X,
% а Y2 и k2 - относительно Y
XC = []; YC = [];
if(abs(1/k1-k2) < 0.1) return; endif % отбрасываем слишком уж параллельные линии
XC = (Y1 + k1*Y2) / (1 - k1*k2);
YC = Y2 + k2*XC;
endfunction
(линии вычисляются в зависимости от того, к какой координатной оси ближе каждая прямая).
В общем, метод себя не оправдал. Буду курочить матричный (итерациями).
Journal information