Документация Engee
Notebook

Обработка изображения звездного неба и получение координат звезд

Открыть пример в Engee

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

Модульная структура Engee и интерактивный редактор скриптов позволит нам быстро разработать прототип алгоритма и оценить его качество.

Входная информация для алгоритма

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

Полезным источником снимков звездного неба, сделанных в любых условиях, является сайт проекта https://nova.astrometry.net/

Типичный кадр со звездного датчика выглядит следующим образом:

photo1713619559.jpeg

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

Прежде всего – вот сокращенный вариант кода, который мы получим после всех испытаний.

In [ ]: 
#using Pkg; Pkg.add( [ "Images", "IndirectArrays", "ImageBinarization", "ImageContrastAdjustment", "ImageDraw", "ImageFiltering", "ImageMorphology", "FFTW",  ], io=devnull )

using Images, ImageDraw, FileIO, FFTW, Statistics, IndirectArrays
using ImageBinarization, ImageContrastAdjustment, ImageFiltering, ImageMorphology
gr( leg=:false )

# Загрузка и предобработка
img = load( "photo1713619559.jpeg" )[1:end, 1:end-2]
img = binarize(img, MinimumError() ) 
#img = adjust_histogram( img, AdaptiveEqualization(nbins = 256, rblocks = 4, cblocks = 4, clip = 0.2))

# Сегментация
bw = Gray.(img) .> 0.7
labels = label_components( bw )
centers = component_centroids( labels )

# Первый график
src_img_temp = load( "photo1713619559.jpeg" )[1:end, 1:end-2]
[ draw!( src_img_temp, Ellipse(CirclePointRadius( Point(Int32.(round.([cy, cx]))...), 10; thickness = 5, fill = false)), RGB{N0f8}(1,0,0) ) for (cx,cy) in centers ]
p1 = plot( src_img_temp, showaxis=:false, grid=:false, ticks=:false, aspect_ratio=:equal )

# Второй график
(w,h), f = floor.(size(img)), 0.05
p2 = plot3d( camera=(60,20) )
for (cx,cy) in centers
    (scx,scy) = (.1 .* [w/h, 1] ./ 2) .* ([cx,cy] .- [w,h]./2) ./ [w,h]./2
    Pi = 1/(sqrt(scx^2 + scy^2 + f^2)) * [f; scy;-scx]
    plot!( p2, [0,Pi[1]], [0,Pi[2]], [0,Pi[3]], leg=:false, xlimits=(-0.1,1.1), ylimits=(-.3,.3), zlimits=(-.3,.3),
        showaxis=:true, grid=:true, ticks=:true, linealpha=.4, markershape=[:none, :star], markersize=[0, 2.5],
        title="Количество звезд: $(length(centers))", titlefontsize=10)
end

plot( p1, p2, size=(800,300) )
Out[0]:

Теперь разберем отдельные шаги разработки этого алгоритма.

Цепочка обработки кадра

Первым делом кадр нужно будет обработать. На каждом шаге обработки мы можем применить несколько разных алгоритмов, но всегда полезно разбить задачу на несколько более простых подзадач, на каждой из которых нужны разные алгоритмы.

star_comp_vision_pipeline.svg

Если на входе корректный кадр (засвеченные или слишком смазанные должен отсеивать другой алгоритм), то нам предстоит выполнить несколько стандартных действий:

  1. избавиться от шума и бинаризовать изображение;
  2. сегментировать объекты чтобы работать с ними по-отдельности;
  3. перевести их координаты в нужную нам систему.

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

Подготовка окружения

Установим нужные библиотеки. Этот шаг нужно выполнить один раз, раскомментируйте его если некоторые библиотеки не загружаются (убрав знаки # в начале строки, или выделив текст и нажав Ctrl+/).

In [ ]: 
# using Pkg; Pkg.add( [ "Images", "IndirectArrays", "ImageBinarization", "ImageContrastAdjustment",
#                       "ImageDraw", "ImageFiltering", "ImageMorphology", "FFTW" ], io=devnull )
In [ ]: 
using Images, ImageDraw, ImageBinarization, ImageContrastAdjustment
using ImageSegmentation, ImageFiltering, ImageMorphology, IndirectArrays
using FileIO, FFTW

Мы будем строить графики типа heatmap и histogram, которые не очень быстро работают при применении бэкенда plotly(). Поэтому мы заранее выберем не интерактивный, зато более быстрый бэкенд, а заодно заранее зададим параметры, общие для всех графиков.

In [ ]: 
gr( showaxis=:false, grid=:false, ticks=:false, leg=:false, aspect_ratio=:equal )
Out[0]:
Plots.GRBackend()

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

Загрузим изображение

In [ ]: 
photoName = "photo1713619559.jpeg";
img = load( photoName )[1:end, 1:end-2];

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

Определим две зоны, которые нас будут интересовать и назовем ее RoI (от Region of Interest, "интересующая область"):

In [ ]: 
RoI = CartesianIndices( (550:800,450:850) );
RoI2 = CartesianIndices( (540:590, 590:660) );

Выбор алгоритма бинаризации изображения

Воспользуемся командами из библиотек ImageBinarization и ImageContrastAdjustment чтобы получить черно-белое (или близкое к нему) изображение, на котором будет удобно сегментировать отдельные звезды.

Первым делом изучим гистограмму изображения звездного неба. У каждого проектировщика будет собственный набор фильтров, которые нужно будет применить в зависимости от шумов и условий засветки.

In [ ]: 
plot(
    plot( img[RoI], title="Фрагмент исходного изображения" ),
    histogram( vec(reinterpret(UInt8, img)), xlimits=(0,256), ylimits=(0,10), xticks=:true, yticks=:false, leg=:false, aspect_ratio=:none, title="Гистограмма"),
    plot( img[RoI2], title="Фрагмент (ближе)" ),
    size=(800,200), layout=grid(1,3, widths=(0.4,0.3,0.3)), titlefont=font(10) 
)
Out[0]:

На гистограмме много пикселей серых оттенков (область 50-100). Кроме того мы видим что на изображении есть артефакты сжатия. Попробуем просто изменить контраст изображения.

In [ ]: 
img = adjust_histogram( Gray.(img), ContrastStretching(t = 0.5, slope = 12) )

plot(
    plot( img[RoI], title="Фрагмент изображения после бинаризации и выравнивания гистограммы" ),
    histogram( vec(reinterpret(UInt8, img)), xlimits=(0,256), ylimits=(0,10), xticks=:true, yticks=:false, leg=:false, aspect_ratio=:none, title="Гистограмма"),
    plot( img[RoI2], title="Фрагмент (ближе)" ),
    size=(800,200), layout=grid(1,3, widths=(0.4,0.3,0.3)), titlefont=font(10) 
)
Out[0]:

Плотность пикселей в области 50-100 стала намного ниже, большинство пикселей теперь имеют яркость 0-20, при этом мы не потеряли яркие пиксели. Путем проб и ошибок мы нашли простой и вполне удовлетворительный метод предобработки, который обрабатывает всё изображение сразу, решает проблему шума и позволяет перейти к бинаризации.

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

Что, если бинаризации по пороговому уровню будет достаточно? Подобрать другой алгоритм бинаризации можно в документации.

In [ ]: 
img = load( photoName )[1:end, 1:end-2]   # Загрузим изображение заново

img = binarize(img, MinimumError() )      # Бинаризация

plot(
    plot( img[RoI], title="Фрагмент исходного изображения" ),
    histogram( vec(reinterpret(UInt8, img)), xlimits=(0,256), ylimits=(0,10), xticks=:true, yticks=:false, leg=:false, aspect_ratio=:none, title="Гистограмма"),
    plot( img[RoI2], title="Фрагмент (ближе)" ),
    size=(800,200), layout=grid(1,3, widths=(0.4,0.3,0.3)), titlefont=font(10) 
)
Out[0]:

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

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

Положим, что мы довольны алгоритмом фильтрации изображения. Перейдем к следующему этапу – сегментации.

Сегментация

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

In [ ]: 
img = load( photoName )[1:end, 1:end-2]   # Загрузим изображение заново

# Бинаризация
img = binarize(img, MinimumError() );
# Гауссовский фильтр
img = imfilter(img, Kernel.gaussian(3));
# Найдем центры при помощи морфологической операции
img_map = Gray.( thinning( img .> 0.1 ) );
# Перечислим все пики их в массиве координат типа CartesianCoordinates
img_coords = findall( img_map .== 1);

# Выведем график
plot( img, title="Обнаружено звезд: $(length(img_coords))", titlefont=font(11) )
scatter!( last.(Tuple.(img_coords)), first.(Tuple.(img_coords)), c=:yellow, markershape=:star, markersize=4 )
Out[0]:

Неплохой результат, хотя:

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

Выводом этой процедуры можно пользоваться как еще одной метрокой качества процедуры предобработки.

Разработаем эталонную процедуру, относительно которой мы будем проверять упрощённые решения. Если бы наши объекты перекрывались, нам бы пригодился метод сегментации методом водораздела#%D0%A1%D0%B5%D0%B3%D0%BC%D0%B5%D0%BD%D1%82%D0%B0%D1%86%D0%B8%D1%8F_%D0%BC%D0%B5%D1%82%D0%BE%D0%B4%D0%BE%D0%BC_%D0%B2%D0%BE%D0%B4%D0%BE%D1%80%D0%B0%D0%B7%D0%B4%D0%B5%D0%BB%D0%B0). Его реализацию можно найти в другом демонстрационном примере: https://engee.com/helpcenter/stable/interactive-scripts/image_processing/Indexing_and_segmentation.html.

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

In [ ]: 
img = load( photoName )[1:end, 1:end-2]   # Повторим загрузку и предобработку
img = binarize(img, MinimumError() )

# Сегментация
bw = Gray.(img) .> 0.7;                     # Еще раз бинаризуем изображение (нам нужна маска)
labels = label_components( bw );            # Пронумеруем все связанные области (маркеры – целые числа)

Из объекта segments мы можем извлечь области, размеры и координаты каждого объекта при помощи функций библиотеки ImageMorphology (таких как component_boxes, component_centroids...).

In [ ]: 
# Центр каждого объекта
centers = component_centroids( labels );

Чтобы визуально оценить качество работы алгоритма, обозначим обнаруженные объекты на графике.

In [ ]: 
src_img_temp = load( photoName )[1:end, 1:end-2]

for (cy,cx) in centers
    draw!( src_img_temp, Ellipse(CirclePointRadius(Int64(round(cx)), Int64(round(cy)), 10; thickness = 5, fill = false)), RGB{N0f8}(1,0,0) )
end

plot( src_img_temp, title="Обнаружено звезд: $(length(centers))", titlefont=font(11))
Out[0]:

Чтобы поближе рассмотреть этот график, можно сохранить его в отдельный файл.

In [ ]: 
save( "$(photoName)_detection.png", src_img_temp )

Определение координат звезд

Последним этапом мы переведем координаты звезд в сферическую систему координат. Она по-прежнему привязана к осям камеры, поэтому мы можем называть ее приборной системой координат, чтобы не путать с другими – геоцентрическими, экваториальными и т.д..

Пользуясь формулой из этой статьи, мы осуществим проекцию пикселей из системы координат сенсора $m_i$ на единичную сферу, на которой радиус-вектор точки из реального мира $P_i$, которую сфотографировала камера, будет равен 1. Координаты пикселей будут заданы относительно центр кадра (w/2, h/2) переменными $u_i$, $v_i$, в проективной модели также используется фокусное расстояние $f$.

$$P_i \rightarrow m_i = \frac{P_i}{\lVert P_i \rVert} = \frac{1}{\sqrt{u_i^2 + v_i^2 + f^2}} \begin{bmatrix} u_i\\ v_i\\ f \end{bmatrix}$$

Здесь мы выбрали реалистичное фокусное расстояние (50 мм), но совершенно нереалистичный размер сенсора (10 см по длинной стороне), по крайней мере для современной встраиваемой аппаратуры. Это сделано для более наглядной визуализации.

In [ ]: 
f = .05; # Фокусное расстояние в метрах
w,h = floor.(size(img))
sensor_size = .1 .* [w/h, 1]

plot3d( )

# Обработаем все найденные центроиды звезд
for (cx,cy) in centers
    # Вычислим координаты пикселя в системе координат сенсора
    (scx,scy) = (sensor_size ./ 2) .* ([cx,cy] .- [w,h]./2) ./ [w,h]./2
    # Координаты точки на единичной сфере
    Pi = 1/(sqrt(scx^2 + scy^2 + f^2)) * [f; scy;-scx]
    # Выведем на графике очередной радиус-вектор
    plot!( [0,Pi[1]], [0,Pi[2]], [0,Pi[3]], leg=:false,
        xlimits=(-0.1,1.1), ylimits=(-.5,.5), zlimits=(-.5,.5),
        showaxis=:true, grid=:true, ticks=:true, linealpha=.4,
        markershape=[:none, :star], markersize=[0, 2.5])
end

plot!(size = (700,500), camera=(60,20) )
Out[0]:

Мы получили проекцию центроидов звезд на единичной сферу в приборной системе координат.

Заключение

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

После затратной стадии прототипирования алгоритм легко перевести в более эффективный вид (например, переписать на С/С++) и разместить в модельно-ориентированном окружении в качестве блока, из которого можно будет сгенерировать код. Другим направлением развития этой работы может быть разработка модели погрешностей датчика, которую легко воплотить стандартными направленными блоками на холсте Engee.