Main Content

Finding the Length of a Pendulum in Motion

This example shows you how to calculate the length of a pendulum in motion. You can capture images in a time series with the Image Acquisition Toolbox™ and analyze them with the Image Processing Toolbox™.

Step 1: Acquire Images

Load the image frames of a pendulum in motion. The frames in the MAT-filependulum.matwere acquired using the following functions in the 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-fileloadpendulum;

Step 2: Explore Sequence with IMPLAY

Run the following command to explore the image sequence inimplay.

implay(frames);

Step 3: Select Region where Pendulum is Swinging

You can see that the pendulum is swinging in the upper half of each frame in the image series. Create a new series of frames that contains only the region where the pendulum is swinging.

To crop a series of frames usingimcrop, first performimcropon one frame and store its output image. Then use the previous output's size to create a series of frame regions. For convenience, use therectthat was loaded bypendulum.matinimcrop.

nFrames = size(frames,4); first_frame = frames(:,:,:,1); first_region = imcrop(first_frame,rect); frame_regions = repmat(uint8(0), [size(first_region) nFrames]);forcount = 1:nFrames frame_regions(:,:,:,count) = imcrop(frames(:,:,:,count),rect);endimshow(frames(:,:,:,1))

步骤4:段每一帧的钟摆

Notice that the pendulum is much darker than the background. You can segment the pendulum in each frame by converting the frame to grayscale, thresholding it usingimbinarize, and removing background structures usingimopenandimclearborder.

initialize array to contain the segmented pendulum frames.

seg_pend = false([size(first_region,1) size(first_region,2) nFrames]); centroids = zeros(nFrames,2); se_disk = strel('disk',3);forcount = 1:nFrames fr = frame_regions(:,:,:,count); gfr = rgb2gray(fr); gfr = imcomplement(gfr); bw = imbinarize(gfr,.7);% threshold is determined experimentallybw = imopen(bw,se_disk); bw = imclearborder(bw); seg_pend(:,:,count) = bw; montage({fr,labeloverlay(gfr,bw)}); pause(0.2)end

Step 5: Find the Center of the Segmented Pendulum in Each Frame

You can see that the shape of the pendulum varied in different frames. This is not a serious issue because you just need its center. You will use the pendulum centers to find the length of the pendulum.

Useregionpropsto calculate the center of the pendulum.

pend_centers = zeros(nFrames,2);forcount = 1:nFrames property = regionprops(seg_pend(:,:,count),'Centroid'); pend_centers(count,:) = property.Centroid;end

Display pendulum centers usingplot.

x = pend_centers(:,1); y = pend_centers(:,2); figure plot(x,y,'m.') axisijaxisequalholdon; xlabel('x'); ylabel('y'); title('Pendulum Centers');

Step 6: Calculate Radius by Fitting a Circle Through Pendulum Centers

Rewrite the basic equation of a circle:

(x-xc)^2 + (y-yc)^2 = radius^2

where(xc、yc)is the center, in terms of parametersa,b,cas

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

wherea = -2*xc,b = -2*yc, andc = xc^2 + yc^2 - radius^2.

You can solve for parametersa,b, andcusing the least squares method. Rewrite the above equation as

a*x + b*y + c = -(x^2 + y^2)

which can also be rewritten as

[x y 1] * [a;b;c] = -x^2 - y^2.

Solve this equation using the backslash(\) operator.

The circle radius is the length of the pendulum in pixels.

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

Superimpose circle and circle center on the plot of pendulum centers.

circle_theta =π/ 3:0.01:π* 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));

See Also

|||||