Синтез изображения для оптической системы¶
В этом примере мы создадим тестовый видеоролик для оптического датчика, объединив ряд синтезированных нами изображений визуальной мишени. Модель движения будет задана в форме динамической модели Engee, а всю обработку мы выполним в скрипте.
Организация процесса синтеза видеоролика¶
Для отладки моделей оптических систем нужны синтетические изображения. В этой демонстрации мы синтезируем такое изображение для модели движения, заданной при помощи графических блоков в файле aero_vibration.engee
.
Мы создадим изображение размером 64 на 64 пикселя и сохраним ролик длиной в 80 кадров (с частотой 10 кадров в секунду).
Изображение будет состоять из мишени, находящейся на подвижном фоне, а точка наблюдения будет подвержена случайной вибрации. Характер движения точки наблюдения будет задан моделью aero_vibration
, приложенной к примеру. Вдобавок изображение будет размыто при помощи гауссовского фильтра (в виде функции рассеяния точки).
delt = 0.1; # Период дискретизации модели
num_frames= 140; # Количество кадров, которые нужно синтезировать
framesize = 64; # Размер стороны изображения в пикселях
# Создадим 3D матрицу для хранения стопки 2D изображений по вектору времени
out = zeros( framesize, framesize, num_frames);
Создание мишени и ее модели движения¶
Первым делом мы зададим внешний вид мишени и модель ее движения. В качестве мишени мы выбрали заметный знак "плюс", а полученная матрица хранит значения интенсивности его пикселей. Согласно заданной модели движения, мишень будет двигаться от центра изображения к его нижнему правому углу.
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="Изображение мишени" )
Совмещение изображения мишени и фона¶
Создадим фон из двух синусоид и зададим модель его движения. Затем наложим поверх него изображение мишени.
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="Фоновое изображение + шум" )
Моделируем угловую вибрацию камеры¶
Вибрация сенсора выражается в смещении изображения камеры. Временной ряд, который задает описание вибрации камеры, рассчитываются в модели aero_vibration
.
Мы запустим модель при помощи команд программного управления (переменная delt
задается в начале этого примера).
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
. Все кадры будут разными, поскольку изображение формируется из движущейся цели и фона, а камера вибрирует. Ниже мы покажем первый кадр этого ролика.
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="Совмещенное изображение" )
Модель оптической системы (гауссовская "модель диафрагмы")¶
Модель оптической системы здесь сводится к размытию изображения в частотной области (FFT). Есть более простые способы решить такую задачу (применить размытие). Но предложенный здесь метод позволяет реализовать и гораздо более сложную "модель диафрагмы" – нужно только заменить первые 5 строчек следующей ячейки на вашу функцию.
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="Размытое изображение" )
Создание и воспроизведение ролика¶
И наконец мы масштабируем изображение таким образом, чтобы значение интенсивности пикселей находилось между 0 и 64 и выводим ролик.
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, fps=10 )