[问题] 球体角加速度实时测量程序

[875691156]

最终生成的函数图像不对，请问如何让改正

clc%清空命令窗口
clear%清空工作区
close all%把显示的图像都关掉


xyloObj = VideoReader('5.mp4');%读入视频
nFrames = xyloObj.NumberOfFrames;%读取视频并逐帧保存为图片
vidHeight = xyloObj.Height;%视频每一帧图像的高度
vidWidth = xyloObj.Width;%视频每一帧图像的宽度
time=xyloObj.Duration/nFrames;%每一帧度过的时间
%%%%%%%%%%%%%求风扇中心的坐标
im = read(xyloObj, 1);%读入第一帧，来获得中心点坐标
im1=im(:,:,2)-im(:,:,1);%图像转化为灰度图像
im1=imbinarize(im1);%灰度图像图像二值化
im1 = bwareaopen(im1,8000);%去除二值化图中的噪点
figure,imshow(im1),title('二值化图')%显示图像
se = strel('disk',35);%设置图像膨胀的参数
img = imdilate(im1,se);%图像膨胀
[L1,N1]=bwlabel(~img); %读取连通区域，联通区域个数=N1；
figure, imshow(L1),title('连通区域标号')
%%%%%%%%%%%%%%%风扇中心坐标为[core_x,core_y]
if length(find(L1==1))>length(find(L1==2))
%一共两个连通区域，找到连通区域较小的那个区域，就是中间的圆，求圆的中心点
[a,c]=find(L1==2);%a为圆内所有点的行数，c为圆内所有点的列数
core_x=(max(a)+min(a))/2;%行数最小值和最大值的平均数为圆心的横坐标
core_y=(max(c)+min(c))/2;
else
[a,c]=find(L1==1);
core_x=(max(a)+min(a))/2;
core_y=(max(c)+min(c))/2;
end
figure,imshow(im1)
draw_point(core_x,core_y),title('风扇中心点')%画出风扇中心点
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for n=1:nFrames
I{n} = read(xyloObj, n);
end
%从第一帧开始把视频中所有的帧都读取出来，存入I中

for n=1:nFrames-1%循环语句，从1开始，到（帧数-1）结束
disp(num2str(n)) %在命令行窗口显示计算到第几帧
[im1,i1,j1]=locate_fan(I{n});%调用locate_fan函数，作用是定位到风扇的扇叶区域
[im2,i2,j2]=locate_fan(I{n+1});
alpha1(n)=cal_angel(i1,j1,core_x,core_y);%求第一幅图的角度，调用cal_angel函数
alpha2(n)=cal_angel(i2,j2,core_x,core_y);%求第二幅图的角度
speed(n)=abs(alpha2(n)-alpha1(n))/time;%瞬时角速度=角位移（单位是弧度）除以对应时间
%这里面有很多帧，我挑了两帧显示一下效果
if n==41
figure,imagesc(im1),colormap gray%显示图像
draw_point(i1,j1) %画出扇叶中心
draw_point(core_x,core_y)%画出电扇中心
%（我们求的就是扇叶中心和电扇中心之间的夹角）
title('图1的角度')
figure,imagesc(im2),colormap gray
draw_point(i2,j2)
draw_point(core_x,core_y)
title('图2的角度')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

speed_jiasudu=abs(diff(speed))/time;%角加速度=角速度（单位是弧度/秒）除以时间
disp(speed_jiasudu);% 显示文本或数组
figure,plot(speed_jiasudu)