Engee 文档
Notebook

处理星空图像并获取恒星坐标

在本演示中，我们将展示一种算法的开发过程，该算法用于检测图像中恒星的中心点，并将其转换为球面坐标系。

通过Engee的模块化结构和交互式脚本编辑器，我们可以快速开发出算法原型，并对其质量进行评估。

算法的输入信息

通过一张星空照片，我们可以计算出相机在天体坐标系中的方位（如赤道）。在此过程中，我们必须创建一个图像处理程序，该程序对在移动平台上进行弱光摄影时固有的各种误差不敏感。这些误差包括传感器热噪声、图像模糊、像素过饱和、外部光源背光等。
<br

该项目网站 https://nova.astrometry.net/

是在各种条件下拍摄星空图像的有用来源。

星空传感器的典型图像如下：

photo1713619559.jpeg

我们会看到一些噪点、不同的亮度区（传感器或初级处理错误）以及运动模糊。所有这些都不会妨碍我们找到恒星的*中心点，然后在知道相机参数的情况下，确定每颗恒星在相机所附球面坐标系中的坐标。

In [ ]: 
Pkg.add(["Statistics", "ImageContrastAdjustment", "ImageSegmentation", "ImageBinarization", "ImageDraw", "IndirectArrays", "ImageFiltering", "ImageMorphology"])
In [ ]: 
# Установка библиотек: закомментируйте при повторном запуске
using Pkg; Pkg.add( [ "Images", "IndirectArrays", "ImageBinarization", "ImageContrastAdjustment", "ImageDraw", "ImageFiltering", "ImageMorphology", "FFTW",  ], io=devnull )

首先，这里是经过所有测试后我们将得到的代码的缩略版。

In [ ]: 
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，但我们并没有丢失任何明亮的像素。通过反复试验，我们找到了一种简单且相当令人满意的预处理方法，它可以一次性处理整个图像，解决噪点问题，并允许我们继续进行二值化处理。

不过，正如进一步的实验所显示的那样，由于我们采用的是全局操作（对整个图像），一些星星变得几乎看不见了，或者在运动模糊形成的线条中出现了不连续性。

如果通过阈值进行二值化就足够了呢？您可以在 documentation 中找到另一种二值化算法。

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]:

在现阶段，我们只需要避免合并对象。在额外的实验中，一些操作（如 morphological from the libraryImageMorphology ）会明显改变框架上线条的方向，并可能改变其中心点的坐标。

<br

如果某些恒星未被识别或其形态发生了变化，我们可以在后期与恒星目录进行比较时检测到它们，为了使每颗恒星的定位精度达到亚像素级，我们可以对找到的每幅图像分别进行处理。

假设我们对图像过滤算法感到满意。让我们进入下一步--分割。

分割

在这一步中，我们必须从恒星图像的二维矩阵中提取恒星的中心点。最简单的方法是模糊图像并提取局部最大值。然而，为了使其有效，我们必须选择大量参数，并确保其对算法应用的每个特定场景中出现的干扰具有鲁棒性。

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 可能会返回一条涂抹轨迹，而不是一个点，每条涂抹线都会被识别为一颗单独的恒星。

该程序的输出结果可作为衡量预处理程序质量的另一个指标。

让我们制定一个参考程序，用来检查简化后的解决方案。如果对象重叠，我们可以使用[分水岭分割]方法(https://ru.wikipedia.org/wiki/%D0%A1%D0%B5%D0%B3%D0%BC%D0%B5%D0%BD%D1%82%D0%B0%D1%86%D0%B8%D1%8F_(%D0%BE%D0%B1%D1%80%D0%B0%D0%B1%D0%BE%D1%82%D0%BA%D0%B0_%D0%B8%D0%B7%D0%BE%D0%B1%D1%80%D0%B0%D0%B6%D0%B5%D0%BD%D0%B8%D0%B9)#%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。

但是，由于我们的对象不会重叠，而且可以很好地相互区分，因此我们需要使用ImageMorphology 中的函数label_components 。

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。像素的坐标将通过变量$u_i$,$v_i$ 设置为相对于画面中心 (w/2,h/2) ，投影模型还使用焦距$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]:

我们将恒星的中心点投影到仪器坐标系中的单位球面上。

结论

我们成功地处理了来自星空传感器的帧，并将恒星坐标从平面坐标系转换到球面坐标系。为此，我们尝试了不同库中已经实现的几种算法。在交互式脚本的不同代码单元中对算法的不同步骤进行处理的方法，使我们能够更快地进行大量计算实验，输出非常复杂和清晰的图解，从而加快了研究进程，并使我们能够更快地找到令人满意的解决过滤恒星图像问题的方法。

在耗资巨大的原型设计阶段之后，该算法可以很容易地转化为更有效的形式（例如用 C/C++ 重写），并作为一个代码块放置在面向模型的环境中。这项工作的另一个发展方向是开发传感器误差模型，这可以很容易地通过Engee画布中的标准方向块来体现。