光学系统的图像合成¶
在本示例中，我们将通过合成视觉目标的大量图像，为光学传感器创建一个测试视频。运动模型将以动态模型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 中进行计算。
我们将使用程序控制命令运行模型（变量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 )
Out[0]:
让我们显示另一张图表，比较物体和背景的亮度水平，加快可视化速度。
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]:
结论¶
我们已经成功合成了一个小剪辑，我们可以在上面测试光学系统的某些任务--跟踪或物体识别。我们使用多个图形块定义了物体运动模型，并通过交互式脚本进行了图像合成。