www.gusucode.com > images 案例代码 matlab源码程序 > images/PendulumLengthExample.m

    %% 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(TM) and analyze them with the Image Processing Toolbox(TM).

% Copyright 1993-2015 The MathWorks, Inc.

%% Step 1: Acquire Images
% Load the image frames of a pendulum in motion. The frames in the MAT-file
% |pendulum.mat| were 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-file

load pendulum;

%% Step 2: Explore Sequence with IMPLAY
% Run the following command to explore the image sequence in |implay|.

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 using |imcrop|, first perform |imcrop| on one
% frame and store its output image. Then use the previous output's size to
% create a series of frame regions.  For convenience, use the |rect| that
% was loaded by |pendulum.mat| in |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))

%%
imshow(frame_regions(:,:,:,1));

%% Step 4: Segment the Pendulum in Each Frame
% 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 using |imbinarize|, and removing background structures
% using |imopen| and |imclearborder|.

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

for count = 1:nFrames
    fr = frame_regions(:,:,:,count);
    imshow(fr)
    pause(0.2)
    
    gfr = rgb2gray(fr);
    gfr = imcomplement(gfr);
    imshow(gfr)
    pause(0.2)
    
    bw = imbinarize(gfr,.7);  % threshold is determined experimentally 
    bw = imopen(bw,se_disk);
    bw = imclearborder(bw);
    seg_pend(:,:,count) = bw;
    imshow(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.

%%
% Use |regionprops| to calculate the center of the pendulum.
pend_centers = zeros(nFrames,2);
for count = 1:nFrames
    property = regionprops(seg_pend(:,:,count), 'Centroid');
    pend_centers(count,:) = property.Centroid; 
end

%%
% Display pendulum centers using |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');

%% 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 parameters |a|, |b|, |c| as
%
% |x^2 + y^2 + a*x + b*y + c = 0|
%
% where |a = -2*xc|, |b = -2*yc|, and |c = xc^2 + yc^2 - radius^2|.
%
% You can solve for parameters |a|, |b|, and |c| using 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)

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

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