Image synthesis for an optical system

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

Organizing the video synthesis process

Synthetic images are needed to debug optical system models. In this demo, we synthesize such an image for a motion model defined using graphic blocks in a file. aero_vibration.engee.

We will create a 64-pixel-by-64-pixel image and save an 80-frame video (at 10 frames per second).

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

Pkg.add( ["Interpolations"] )

delt = 0.1;      # Период дискретизации модели
num_frames= 140;  # Количество кадров, которые нужно синтезировать
framesize = 64;  # Размер стороны изображения в пикселях

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

Creating a target and its motion model

First of all, we will set the appearance of the target and the model of its movement. We chose a noticeable plus sign as a target, and the resulting matrix stores the intensity values of its pixels. According to the specified motion model, the target will move from the center of the image to its lower right corner.

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="Изображение мишени" )

Combining the target image and background

Let's create a background of two sinusoids and set a model of its movement. Then we will overlay the image of the target on top of it.

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="Фоновое изображение + шум" )

We simulate the angular vibration of the camera

The vibration of the sensor is expressed in the displacement of the camera image. The time series that defines the camera vibration description is calculated in the model aero_vibration.

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

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)
)

We simulate movement in the frame by combining background, target, and vibration.

Now we will synthesize all the frames of the future video and put them into 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.

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="Совмещенное изображение" )

Optical system model (Gaussian "diaphragm model")

The optical system model here boils down to frequency domain image blurring (FFT). There are simpler ways to solve this problem (apply a blur). But the method proposed here allows you to implement a much more complex "aperture model" – you only need to replace the first 5 lines of the next cell with your function.

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="Размытое изображение" )

Creating and playing a video clip

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

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

Let's draw another graph to compare the brightness levels of the object and the background, speeding up the visualization a little.

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 )

Conclusion

We have successfully synthesized a small video on which we can test some task for the optical system – tracking or object recognition. The model of the object's movement was set by us using several graphic blocks, and image synthesis was performed in an interactive script.