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

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

Этот пример показывает вам, как вычислить длину маятника в движении. Можно получить изображения во временных рядах с 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: Выберите Region, где Маятник Качается

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

Обрезать серию систем координат с помощью 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) центр, в терминах параметров aBC как

x^2 + y^2 + a*x + b*y + c = 0

где a = -2*xc, b = -2*yc, и c = xc^2 + yc^2 - radius^2.

Можно решить для параметров aB, и 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));

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