Продолжение эпопеи с центрами

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

Вот увеличенная область вблизи медиан:

Ну и напоследок файлы. Определение параметров линии теперь делается так:

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

Ну и сама функция определения координат:

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