Синтез изображения для оптической системы¶

Открыть пример в Engee

В этом примере мы создадим тестовый видеоролик для оптического датчика, объединив ряд синтезированных нами изображений визуальной мишени. Модель движения будет задана в форме динамической модели 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 )

[ Info: Saved animation to /user/start/examples/aerospace/optical_sensor_image_generation/tmp.gif

Выведем еще один график чтобы сравнить уровни яркости объекта и фона, немного ускорив визуализацию.

a = Plots.Animation()

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

gif( a, "2.gif", fps=15 )

[ Info: Saved animation to /user/start/examples/aerospace/optical_sensor_image_generation/2.gif

Заключение¶

Мы успешно синтезировали небольшой ролик, на котором можем протестировать какую-нибудь задачу для оптической системы – слежение или распознавание объектов. Модель движения объекта была задана нами при помощи нескольких графических блоков, а синтез изображений был выполнен в интерактивном скрипте.