Емельянов Эдуард Владимирович (eddy_em) wrote,
Емельянов Эдуард Владимирович
eddy_em

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

Итак, продолжая разбираться с 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:
(подогнать уровень для бинаризации я не сильно-то старался, поэтому количество найденных пятен различается).

Tags: fits, octave, поиск связанных областей
Subscribe

  • Дела телескопные

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

  • Robotel-2

    В спешном порядке я нарисовал "демон погоды" для мониторинга облачности, чтобы запустить и второй телескоп. В конце года мы с Тимуром примерно…

  • Масштабы

    Сейчас из дома параллельно провожу технические наблюдения на Z-1000 и полуметровом астросибе. На Цейссе провожу мониторинг (примерно по часу на…

promo eddy_em august 17, 2019 12:33 3
Buy for 10 tokens
Юра намедни напечатал корпус для хронометра. Для первого блина получилось неплохо: И еще немного фотографий:
  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your reply will be screened

  • 0 comments