Engee 文档
Notebook

光学系统的图像合成

在本示例中,我们将通过合成视觉目标的大量图像,为光学传感器创建一个测试视频。运动模型将以动态模型Engee的形式指定,我们将在脚本中完成所有处理。

视频合成过程的组织

调试光学系统模型需要合成图像。在本演示中,我们将为使用aero_vibration.engee 文件中的图形块定义的运动模型合成这样的图像。

我们将创建一幅 64 x 64 像素的图像,并保存一个 80 帧的剪辑(每秒 10 帧)。

图像将由移动背景上的一个目标组成,观察点将受到随机振动。观察点的运动模式将由示例所附的aero_vibration 模型给出。此外,还将使用高斯滤波器(作为点的散射函数)对图像进行模糊处理。

In [ ]:
Pkg.add( ["Interpolations"] )
In [ ]:
delt = 0.1;      # Период дискретизации модели
num_frames= 140;  # Количество кадров, которые нужно синтезировать
framesize = 64;  # Размер стороны изображения в пикселях

# Создадим 3D матрицу для хранения стопки 2D изображений по вектору времени
out = zeros( framesize, framesize, num_frames);

创建目标及其运动模型

我们首先要做的是设置目标的外观及其运动模型。我们选择了一个可见的 "正 "号作为目标,由此产生的矩阵存储了其像素的强度值。根据指定的运动模型,目标将从图像中心向右下角移动。

In [ ]:
target = [zeros(3,11)
          zeros(1,5) 6 zeros(1,5)
          zeros(1,5) 6 zeros(1,5)
          zeros(1,3) 6 6 6 6 6 zeros(1,3) # В качестве мишени мы изобразим знак "плюс" размером 5 на 5,
          zeros(1,5) 6 zeros(1,5)         #  интенсивность 6 выбрана ради соотношения сигнал/шум ~4
          zeros(1,5) 6 zeros(1,5)         # Мы делаем изображение мишени немного больше (11x11 пикселей)
          zeros(3,11)];                   #  чтобы работать с этим изображением без ошибок интерполяции

target_velx = 1;                  # скорость движения мишени по оси X (пиксели/сек)
target_vely = 1;                  # скорость движения мишени по оси Y (пиксели/сек)
target_x_initially = framesize÷2; # изначально мишень находится по центру изображения по обеим осям
target_y_initially = framesize÷2; # x и y

gr()
heatmap( target.*32,
         c=:grays, framestyle=:none, aspect_ratio=:equal, leg=:false,
         title="Изображение мишени" )
Out[0]:

对齐目标图像和背景

创建一个由两个正弦波组成的背景,并定义其运动模型。然后在其上叠加目标图像。

In [ ]:
using FFTW

backsize = framesize+36;  # Фоновое изображение значительно больше изображения мишени
                          # поскольку мы не хотим потерять мишень из виду в процессе движения
xygrid = (1:backsize) ./ backsize;
B = (2sin.(2pi .* xygrid).^2) * (cos.(2pi .* xygrid).^2)';

psd = fft( B );
psd = real( psd .* conj(psd) );

background = B .+ 0.5*randn(backsize, backsize);    # Добавим гауссовский шум
					                                # с вариацией 0.25 (СКО 0.5)

xoff = 10;     # Изначальное смещение точки начала координат сенсора
yoff = 10;     #  относительно точки (0,0) фонового изображения

driftx = 2;    # Скорость смещения фонового изображения (пиксели/сек)
drifty = 2;

minout = minimum( background );
maxout = maximum( background );

heatmap( (background.-minout) .* 64 ./ (maxout.-minout),
         c=:grays, framestyle=:none, aspect_ratio=:equal, leg=:false, clim=(0,255),
		 title="Фоновое изображение + шум" )
Out[0]:

摄像机角度振动建模

传感器的振动表现为摄像机图像的位移。说明摄像机振动的时间序列在模型aero_vibration 中进行计算。

image.png

我们将使用程序控制命令运行模型(变量delt 已在本例开始时设置)。

In [ ]:
omega = 2*pi*5;       # Вибрация структуры складывается из колебаний 5, 10 и 15 Hz
zeta  = 0.01;         # Коэффициент демпфирования для всех компонентов (мод)

# Откроем и выполним модель
modelName = "aero_vibration";
model = modelName in [m.name for m in engee.get_all_models()] ? engee.open( modelName ) : engee.load( "$(@__DIR__)/$(modelName).engee");
engee.set_param!( modelName, "FixedStep"=>delt );
engee.set_param!( modelName, "StopTime"=>num_frames*delt*2 );

data = engee.run( modelName, verbose=false );

# Результат запуска модели разделим на два вектора
vibtime = data["vibdat"].time
vibdata = data["vibdat"].value
vibLen = length( vibdata )

vibx = vibdata[ 1:vibLen÷2 ]
viby = vibdata[ (vibLen÷2+1):end ]

levarmx = 10;   # Длина рычага, переводящего вращательную вибрацию в продольное движение 
levarmy = 10;   #  по осям X и Y

plot(
    plot( vibtime[1:vibLen÷2], vibx, title="Динамика углового движения опоры на которой\n размещен оптический датчик (рад)",
           ylabel="по оси x", grid=true, titlefont=font(10), leg=:false ),
    plot( vibtime[(vibLen÷2+1):end], viby, xlabel="Время, с", ylabel="по оси y", grid=true, leg=:false ),
    layout=(2,1)
)
Out[0]:

结合背景、目标和振动模拟帧中的运动

现在,我们将合成未来视频的所有帧,并将其放入矩阵out 中。由于图像是由移动的目标和背景以及摄像机振动形成的,因此所有帧都会有所不同。下面我们将展示这段视频片段的第一帧。

In [ ]:
using Interpolations

for t = 1:num_frames
    
    # Смещение фона со скоростью driftx и drifty и добавление вибрации
    xshift = driftx*delt*t + levarmx * vibx[t];
    yshift = drifty*delt*t + levarmy * viby[t];
    
    # Интерполяция изображения фона (субпиксельное перемещение изображения)
    x2,y2 = ( xshift:1:xshift+backsize, yshift:1:yshift+backsize )
    itp = LinearInterpolation((1:backsize, 1:backsize), background, extrapolation_bc=Line())
    outtemp = [ itp(x,y) for y in y2, x in x2 ]
    
    # Обрезка изображения до размера framesize (от размера backsize)
    out[:,:,t] = outtemp[xoff:xoff+framesize-1, yoff:yoff+framesize-1];
    
    # Добавим движение мишени
    tpixinx = floor( target_velx * delt * t );  # Сохраним амплитуду смещения перед интерполяцией
    tpixiny = floor( target_vely * delt * t );  #  в целых пикселях
    
    # Интерполяция значений яркости
    txi = target_velx*delt*t - tpixinx;
    tyi = target_vely*delt*t - tpixiny;
    txi,tyi = ( txi+1:txi+11, tyi+1:tyi+11)     # Интерполяция мишени нужна только относительно ее центра
    
    itp = LinearInterpolation( (1:11, 1:11), target, extrapolation_bc=Line() )
    temp = [ itp(x,y) for y in tyi, x in txi ]
    
    # Поместить изображение мишени в точке, определенной исходным смещением камеры и сохраненной ранее амплитудой смещения
    tx = convert(Int64, tpixinx + target_x_initially - 1);
    ty = convert(Int64, tpixiny + target_y_initially - 1);

    out[tx:tx+6, ty:ty+6, t] = temp[9:-1:3, 9:-1:3] + out[tx:tx+6,ty:ty+6,t];
    
end

minout = minimum( out );
maxout = maximum( out );

heatmap( (out[:,:,1].-minout) .* 64 ./ (maxout.-minout),
         c=:grays, framestyle=:none, aspect_ratio=:equal, leg=:false, clim=(0,255),
		 title="Совмещенное изображение" )
Out[0]:

光学系统模型(高斯 "光圈模型)

这里的光学系统模型简化为在频域(FFT)中模糊图像。有更简单的方法来解决这个问题(应用模糊)。但这里提出的方法可以实现更复杂的 "虹膜模型"--您只需用您的函数替换下一个单元格的前 5 行。

In [ ]:
x = 1:framesize;
y = 1:framesize;
sigma = 120;
apfunction = @. exp(-(x-framesize/2).^2/(2*sigma))' * exp(-(y-framesize/2).^2/(2*sigma));
apfunction = fftshift(apfunction);

for j = 1:num_frames
    out[:,:,j] = real.( ifft( apfunction.*fft( out[:,:,j] )));
end

minout = minimum(out);
maxout = maximum(out);

heatmap( (out[:,:,1].-minout) .* 64 ./ (maxout.-minout),
         c=:grays, framestyle=:none, aspect_ratio=:equal, leg=:false, clim=(0,255),
		 title="Размытое изображение" )
Out[0]:

创建并播放剪辑

最后,我们缩放图像,使像素强度值介于 0 和 64 之间,然后输出剪辑。

In [ ]:
minout = minimum(out);
maxout = maximum(out);

a = Plots.Animation()

for j = 1:num_frames
    plt = heatmap( (out[:,:,j] .- minout)*64/(maxout-minout),
        c=:grays, framestyle=:none, aspect_ratio=:equal, leg=:false, clim=(0,255),
        size=(400,400), yflip=true )
    frame( a, plt )
end

gif( a, "1.gif", fps=15 )
[ Info: Saved animation to /user/start/examples/aerospace/optical_sensor_image_generation/1.gif
Out[0]:
No description has been provided for this image

让我们显示另一张图表,比较物体和背景的亮度水平,加快可视化速度。

In [ ]:
a = Plots.Animation()

for j = 1:num_frames
    plt = wireframe( (out[:,:,j] .- minout)*64/(maxout-minout),
        c=:grays, aspect_ratio=:equal, leg=:false,
        size=(500,500), yflip=true, zlimits=(0,64), grid=true )
    frame( a, plt )
end

gif( a, "2.gif", fps=15 )
Out[0]:
No description has been provided for this image

结论

我们已经成功合成了一个小剪辑,我们可以在上面测试光学系统的某些任务--跟踪或物体识别。我们使用多个图形块定义了物体运动模型,并通过交互式脚本进行了图像合成。

示例中使用的块