И.М.Журавель "Краткий курс теории обработки изображений"

В оглавление \ К следующему разделу \ К предыдущему разделу

Поиск длины маятника в движении

Проведем захват серии изображений с помощью пакета Image Acquisition Toolbox и затем проанализируем их с помощью Image Processing Toolbox. Рассмотрим это на примере задачи по вычислению длины маятника в движении.

Содержание

• Шаг 1: Захват изображений.
• Шаг 2: Выбор области колебания маятника.
• Шаг 3: Сегментация маятника в каждом фрейме.
• Шаг 4: Поиск центра сегментированного маятника в каждом фрейме.
• Шаг 5: Вычисление радиуса колебания маятника.

Шаг 1: Захват изображений.

Считываем кадры изображений маятника в движении. Кадры изображений, после их захвата с помощью 
Image Acquisition Toolbox, хранятся в виде MAT-файла pendulum.mat.

% Доступ к захваченным изображениям (видеообъектам).
% vidobj=videoinput('winvideo', 1, 'RGB24_352x288');
% Формирование объекта на основе пяти фреймов.
% set(vidobj, 'FrameGrabInterval', 5);
% Формирование большего числа фреймов.
% nFrames=50;
% set(vidobj, 'FramesPerTrigger', nFrames);
% Доступ к видеоустройству.
% src=getselectedsource(vidobj);
% Настройка исходного устройства на формирование 30 кадров в секунду.
% set(src, 'FrameRate', '30');
% Открытие окна просмотра.
% preview(vidobj);
% Начало захвата изображений.
% start(vidobj);
% Определение окончания захвата данных.
% wait(vidobj, 10);
% Извлечение фреймов из памяти.
% frames=getdata(vidobj);
% Удаление рабочих переменных.
% delete(vidobj)
% clear vidobj
% Считывание MAT-файла.
load pendulum;
immovie(frames);

Шаг 2: Выбор области колебания маятника.

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

При создании новой серии изображений используем функцию imcrop.

nFrames=size(frames, 4);
first_frame=frames(:, :, :, 1);
first_region=imcrop(first_frame, rect);
frame_regions=repmat(uint8(0), [size(first_region) nFrames]);
for count=1:nFrames
  frame_regions(:, :, :, count)=imcrop(frames(:, :, :, count), rect);
end
immovie(frame_regions);

Шаг 3: Сегментация маятника в каждом фрейме.

Отметим, что маятник представлен в виде потемнения на фоне. Операцию сегментирования можно провести, преобразив фрейм в полутоновое изображение и далее использовать функцию im2bw с некоторым порогом. Далее нужно удалить фон с помощью функций imopen и imclearborder.

% Получение массива, содержащего сегментированное изображение маятника.

seg_pend=false([size(first_region, 1) size(first_region, 2) nFrames]);
centroids=zeros(nFrames, 2);
se_disk=strel('disk', 3);
for count=1:nFrames
    fr=frame_regions(:, :, :, count);
    imshow(fr)
    pause(0.2)
    gfr=rgb2gray(fr);
    gfr=imcomplement(gfr);
    imshow(gfr)
    pause(0.2)
    bw=im2bw(gfr, .7);  % пороговая обработка
    bw=imopen(bw, se_disk);
    bw=imclearborder(bw);
    seg_pend(:, :, count)=bw;
    imshow(bw)
    pause(0.2)
end

Шаг 4: Поиск центра сегментированного маятника в каждом фрейме.

Мы можем наблюдать как изменяется форма и расположение маятника в различных фреймах. Теперь необходимо решить непростую задачу поиска центра маятника. Информацию, найденную при решении этой задачи, используем для определения длины маятника. Для вычисления центра маятника будем использовать функции bwlabel и regionprops.

for count=1:nFrames
    lab=bwlabel(seg_pend(:, :, count));
    property=regionprops(lab, 'Centroid');
    pend_centers(count, :)=property.Centroid; 
end
Отображение центра маятника с помощью функции plot.
x=pend_centers(:, 1);
y=pend_centers(:, 2);
figure
plot(x, y, 'm.'), axis ij, axis equal, hold on;
xlabel('x');
ylabel('y');
title('центр маятника');

Шаг 5: Вычисление радиуса колебания маятника.

Запишем основное выражение окружности: (x-xc)^2+(y-yc)^2=radius^2, где (xc, yc) является центром, параметры a, b, c в уравнении x^2+y^2+a*x+b*y+c=0 представлены в виде a=-2*xc, b=-2*yc и c=xc^2+yc^2-radius^2.

Параметры a, b и c можно найти с помощью метода наименьших квадратов. Тогда выражение можно представить как a*x+b*y+c=-(x^2+y^2) или в виде [a; b; c]*[x y 1]=-x^2-y^2. Представим это выражение в виде операции обратного деления.

Радиус окружности является длиной маятника в пикселях.

abc=[x y ones(length(x), 1)]\[-(x.^2+y.^2)];
a=abc(1); b=abc(2); c=abc(3);
xc=-a/2;
yc=-b/2;
circle_radius=sqrt((xc^2+yc^2)-c);
pendulum_length=round(circle_radius)
pendulum_length=
   253
   
Отобразим окружность и ее центр, совмещая с центром маятника.

circle_theta=pi/3:0.01:pi*2/3;
x_fit=circle_radius*cos(circle_theta)+xc;
y_fit=circle_radius*sin(circle_theta)+yc;
plot(x_fit, y_fit, 'b-');
plot(xc, yc, 'bx', 'LineWidth', 2);
plot([xc x(1)], [yc y(1)], 'b-');
text(xc-110,yc+100,sprintf('длина маятника=%d пикселя', pendulum_length));

В оглавление \ К следующему разделу \ К предыдущему разделу