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

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

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