tiii1.m

其他书记方面 其他书记方面 其他书记方面 其他书记方面 其他书记方面

load pendulum;
immovie(frames);
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
immovie(frame_regions);
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 = im2bw(gfr,.7);  % threshold is determined experimentally 
    bw = imopen(bw,se_disk);
    bw = imclearborder(bw);
    seg_pend(:,:,count) = bw;
    imshow(bw)
    pause(0.2)
end
for count = 1:nFrames
    lab = bwlabel(seg_pend(:,:,count));
    property = regionprops(lab,'Centroid');
    pend_centers(count,:) = property.Centroid; 
end
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');
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)
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));