Engee 文档
Notebook

光学系统的图像合成

在这个例子中，我们将通过组合视觉目标的许多合成图像来为光学传感器创建一个测试视频。 运动模型将以Engee动态模型的形式设置，我们将在脚本中执行所有处理。

组织视频合成过程

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

我们将创建一个64x64像素的图像并保存一个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

结论

我们已经成功地合成了一个小视频，我们可以在其上测试光学系统的一些任务–跟踪或物体识别。 我们使用几个图形块设置对象运动的模型，并在交互式脚本中进行图像合成。

示例中使用的块