Емельянов Эдуард Владимирович (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

  • M$ teams…

    Начал с сегодняшнего дня студентам ЮФУ удаленно лекции читать. У них все завязано на различные корпорации зла. И базовая работа - через teams. ОК,…

  • Почему systemd — дерьмо

    Уже давно на эту статейку натыкался, но все забывал в "закладки" добавить. Вот, добавляю: "systemd — отстой". Советую эту статейку почитать…

  • Разбираемся с утечками

    Занялся полным рефакторингом кода системы управления оптоволоконным спектрографом. А глаз уже "замыленный": ну вот не могу сходу обнаружить, где я…

promo eddy_em september 3, 12:13 8
Buy for 10 tokens
Уже больше полугода занимаюсь разработкой, вот, наконец-то в мастерских взялись за меня и начали выдавать первые детали. Сегодня сделал тестовую сборку (как обычно, местами пришлось "доработать напильником"): Пока прибор без названия (да и как-то не лезет в голову ничего, у меня нет…
  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic
  • 0 comments