Engee documentation
Notebook

Image synthesis for an optical system

In this example, we will create a test video for an optical sensor by combining a number of images of a visual target that we have synthesised. The motion model will be specified in the form of a dynamic model Engee, and we will do all the processing in a script.

Organisation of the video synthesis process

Synthetic images are needed to debug models of optical systems. In this demonstration, we synthesise such an image for a motion model defined using graphical blocks in the file aero_vibration.engee.

We will create an image of 64 by 64 pixels and save a clip of 80 frames (at 10 frames per second).

The image will consist of a target on a moving background, and the observation point will be subject to random vibration. The pattern of motion of the observation point will be given by the model aero_vibration, attached to the example. In addition, the image will be blurred using a Gaussian filter (as a scattering function of the point).

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

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

Creating a target and its motion model

The first thing we will do is to set the appearance of the target and its motion model. As a target we have chosen a visible "plus" sign, and the resulting matrix stores the intensity values of its pixels. According to the specified motion model, the target will move from the centre of the image to its lower right corner.

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]:

Alignment of the target image and the background

Create a background of two sinusoids and define its motion model. Then overlay the target image on top of it.

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]:

Modelling the angular vibration of the camera

The vibration of the sensor is expressed as a displacement of the camera image. The time series that specifies the description of the camera vibration are calculated in the model aero_vibration.

image.png

We will run the model using programme control commands (the variable delt is set at the beginning of this example).

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]:

Simulate motion in the frame by combining background, target and vibration

Now we will synthesise all the frames of the future video and place them in the matrix out. All frames will be different, because the image is formed from a moving target and background, and the camera vibrates. Below we will show the first frame of this video clip.

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]:

Optical system model (Gaussian "aperture model")

The optical system model here is reduced to blurring the image in the frequency domain (FFT). There are simpler ways to solve such a problem (apply blurring). But the method proposed here allows you to realise a much more complex "iris model" - you just need to replace the first 5 lines of the next cell with your function.

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]:

Creating and playing a clip

Finally, we scale the image so that the pixel intensity value is between 0 and 64 and output the clip.

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

Let's display another graph to compare the brightness levels of the object and the background, speeding up the visualisation a bit.

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

Conclusion

We have successfully synthesised a small clip on which we can test some task for an optical system - tracking or object recognition. The model of object movement was defined by us using several graphical blocks, and image synthesis was performed in an interactive script.

Blocks used in example