Вычисление длины маятника в движении

В этом примере показано, как получить и анализировать изображения объекта в движении.

Этот пример получает изображения маятника в движении. Маятник состоит из синего мяча, присоединенного к нейлоновой строке. Данные изображения получены, когда маятник качают. После того, как полученный, изображения обрабатываются, чтобы определить длину маятника.

Изображения получаются с помощью Image Acquisition Toolbox™ и анализируются с Image Processing Toolbox™.

Получите изображения

Получите серию изображений, чтобы анализировать.

% Access an image acquisition device.
vid = videoinput('winvideo', 1, 'RGB24_352x288');
vid.Timeout = 12;

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

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

% Access the device's video source object selected for acquisition.
src = getselectedsource(vid);

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

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

% Initiate the acquisition.
start(vid);

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

% Remove video input object from memory.
delete(vid)
clear vid

% Display the first frame in the series.
imshow( frames(:, :, :, 1) );

% Display all acquired images.
imaqmontage(frames);

Выберите Region Of Interest

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

Обрезать серию систем координат с помощью imcrop, сначала выполните imcrop на одной системе координат и хранилище ее выход. Затем создайте серию систем координат, имеющих размер предыдущего выхода.

% Determine the total number of frames acquired.
nFrames = size(frames, 4);

% Crop the first frame.
roi = [50 16 222 68];
firstFrame = frames(:,:,:,1);
frameRegion = imcrop(firstFrame, roi);

% Create a storage for the modified image series.
regions = repmat(uint8(0), [size(frameRegion) nFrames]);
for count = 1:nFrames,
    regions(:,:,:,count) = imcrop(frames(:,:,:,count), roi);
    imshow(regions(:,:,:,count));
end

Сегментируйте маятник в каждой системе координат

Поскольку маятник является намного более темным, чем фон, маятник может быть сегментирован в каждой системе координат путем преобразования системы координат в шкалу полутонов, пороговая обработка это и фоновые структуры удаления.

% Initialize array to contain the segmented pendulum frames.
segPend = false([size(frameRegion, 1) size(frameRegion, 2) nFrames]);
centroids = zeros(nFrames, 2);
structDisk = strel('disk', 3);

for count = 1:nFrames,
    % Convert to grayscale.
    fr = regions(:,:,:,count);
    gfr = rgb2gray(fr);
    gfr = imcomplement(gfr);

    % Experimentally determine the threshold.
    bw = im2bw(gfr, .7);
    bw = imopen(bw, structDisk);
    bw = imclearborder(bw);
    segPend(:,:,count) = bw;
    imshow(bw);
end

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

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

% Calculate the pendulum centers.
for count = 1:nFrames,
    property = regionprops(segPend(:, :, count), 'Centroid');
    pendCenters(count,:) = property.Centroid;
end

% Display the pendulum centers and adjust the plot.
figure;
x = pendCenters(:, 1);
y = pendCenters(:, 2);
plot(x, y, 'm.');
axis ij;
axis equal;
hold on;
xlabel('x');
ylabel('y');
title('Pendulum Centers');

Вычислите длину маятника

Путем подбора кривой кругу через центры маятника может быть вычислена длина маятника. Перепишите основное уравнение круга:

(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

Это уравнение может быть решено в MATLAB® с помощью обратной косой черты (\) оператор.

% Solve the equation.
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;
circleRadius = sqrt((xc^2 + yc^2) - c);

% Circle radius is the length of the pendulum in pixels.
pendulumLength = round(circleRadius)

pendulumLength =

   253

% Superimpose results onto the pendulum centers
circle_theta = pi/3:0.01:pi*2/3;
x_fit = circleRadius*cos(circle_theta) + xc;
y_fit = circleRadius*sin(circle_theta) + yc;
plot(x_fit, y_fit, 'b-');
plot(xc, yc, 'bx', 'LineWidth', 2);
plot([xc x(1)], [yc y(1)], 'b-');
titleStr = sprintf('Pendulum Length = %d pixels', pendulumLength);
text(xc-110, yc+100, titleStr);

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