Engee 文档
Notebook

对星空进行图像处理，获取星星的坐标

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

Engee的模块化结构和交互式脚本编辑器将使我们能够快速开发原型算法并评估其质量。

算法的输入信息

使用星空的照片，我们可以计算相机在天体坐标系中的方向（例如，экваториальной)。与此同时，我们将不得不创建一个图像处理程序，该程序对在低光照条件下在移动平台上拍摄所固有的各种错误不敏感。 这些因素包括热传感器噪声、模糊图像、过饱和像素、来自外部光源的照明等。

在任何条件下拍摄的星空图像的有用来源是该项目的网站。 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 （来自感兴趣区域，"感兴趣区域"）:

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，而我们没有失去明亮的像素。 通过反复试验，我们找到了一种简单且相当令人满意的预处理方法，该方法一次性处理整个图像，解决了噪声问题并允许我们进行二值化。

然而，正如进一步的实验将向我们展示的那样，由于我们应用了全局操作（对整个图像），一些恒星变得几乎看不见，或者在运动过程中模糊形成的线条中有

如果阈值水平的二值化就足够了呢？ 您可以在[文档]中找到不同的二值化算法（https://engee.com/helpcenter/stable/julia/ImageBinarization/index.html )。

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 它可以返回的不是一个点，而是一个涂抹轨迹，并且这些线中的每一条都将被识别为一个单独的恒星。

此过程的输出可以用作预处理过程质量的另一个度量。

我们将制定一个参考程序，我们将测试简化的解决方案。 如果我们的对象重叠，我们将使用[分割方法]方法。 водораздела](https://ru.wikipedia.org/wiki/Сегментация_(обработка_изображений)#Сегментация_методом_водораздела)。其实现可以在另一个演示示例中找到：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 )

确定恒星的坐标

最后一步是将恒星的坐标转换为球面坐标系。 它仍然与相机的轴相关联，因此我们可以将其称为仪器坐标系统，以免将其与其他人混淆–地心，赤道等。.

使用[本文]中的公式（https://www.researchgate.net/figure/Line-projection-onto-the-unit-sphere_fig2_254987768 ），我们将执行像素从传感器的坐标系统的投影 在一个单位球体上，其半径是来自现实世界的点的矢量 ，这是由相机拍摄，将等于1。 像素坐标将相对于帧的中心设置（w/2, h/2)变量 , 射影模型也使用焦距 .

在这里，我们选择了一个现实的焦距（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画布上用标准方向块实现。