Нахождение длины маятника в движении

Открыть скрипт

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

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

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

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

См. также

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