Продолжаем с Octave и FITS

[eddy_em]

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

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

Загрузив пробный файл (звездное поле), я обнаружил, что функция read_fits_image отображает его неправильно: как бы повернутым на 90°. Поэтому была написана обертка:

function [image header] = fits_read(filename)
% считываем FITS-файл и переворачиваем для правильного отображения
% Координата y отсчитывается в результате сверху! (если не делать flipdim, отображаться будет зеркально)
	[image header] = read_fits_image(filename);
	image = flipdim(image',1); % переворачиваем, чтобы отображалось верно
endfunction

Теперь команда

II = fits_read('object_0001.fit');

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

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

Так как изображения сильно зашумлены, нам придется отфильтровать их. Отфильтруем бинарное изображение медианным фильтром 3х3 пикселя. Это избавит нас от "ложных срабатываний" детектора связанных областей на шумах.

Ну, а в конце-концов, выделив и пронумеровав связанные области на изображении, найдем их центр тяжести по самой наитупейшей формуле.

В результате получается следующая функция:

function [ xc, yc ] = find_centers(Img, treshold)
% Нахождение центров пятен изображения Img с интенсивностью выше treshold
% xc, yc - массивы координат
	tmp = zeros(size(Img));
	tmp(find(Img > treshold)) = 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

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

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