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

Открыть сценарий

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

Шаг 1: Получение изображений

Загрузить кадры изображения маятника в движение. Кадры в MAT-файле pendulum.mat были получены с помощью следующих функций в панели инструментов получения изображений.

% Access an image acquisition device (video object).
% vidimage=videoinput('winvideo',1,'RGB24_352x288');

% Configure object to capture every fifth frame.
% vidimage.FrameGrabInterval = 5;

% Configure the number of frames to be logged.
% nFrames=50;
% vidimage.FramesPerTrigger = nFrames;

% Access the device's video source.
% src=getselectedsource(vidimage);

% Configure device to provide thirty frames per second.
% src.FrameRate = 30;

% Open a live preview window. Focus camera onto a moving pendulum.
% preview(vidimage);

% Initiate the acquisition.
% start(vidimage);

% Wait for data logging to finish before retrieving the data.
% wait(vidimage, 10);

% Extract frames from memory.
% frames = getdata(vidimage);

% Clean up. Delete and clear associated variables.
% delete(vidimage)
% clear vidimage

% load MAT-file

load pendulum;

Шаг 2: Исследовать последовательность с помощью IMPLAY

Выполните следующую команду для просмотра последовательности изображений в implay.

implay(frames);

Шаг 3: Выбор области, в которой маятник качается

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

Обрезка серии кадров с помощью imcrop, сначала выполнить imcrop на одном кадре и сохраняют его выходное изображение. Затем используйте размер предыдущего вывода для создания последовательности областей кадра. Для удобства используйте rect которая была загружена pendulum.mat в 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
imshow(frames(:,:,:,1))

Шаг 4: Сегментация маятника в каждом кадре

Обратите внимание, что маятник намного темнее фона. Можно сегментировать маятник в каждом кадре, преобразовав кадр в градации серого и установив пороговые значения с помощью imbinarizeи удаление фоновых структур с помощью 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);

    gfr = rgb2gray(fr);
    gfr = imcomplement(gfr);

    bw = imbinarize(gfr,.7);  % threshold is determined experimentally
    bw = imopen(bw,se_disk);
    bw = imclearborder(bw);
    seg_pend(:,:,count) = bw;

    montage({fr,labeloverlay(gfr,bw)});
    pause(0.2)

end

Шаг 5: Найти центр сегментированного маятника в каждом кадре

Видно, что форма маятника варьировалась в разных рамах. Это не серьезный вопрос, потому что вам просто нужен его центр. Центр маятника используется для определения длины маятника.

Использовать regionprops для расчета центра маятника.

pend_centers = zeros(nFrames,2);
for count = 1:nFrames
    property = regionprops(seg_pend(:,:,count), '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('Pendulum Centers');

Шаг 6: Расчет радиуса путем подгонки окружности через центры маятника

Переписать основное уравнение круга:

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

который также может быть переписан как

[x y 1] * [a;b;c] = -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('Pendulum length = %d pixels', pendulum_length));

См. также

Документация