Matlab描绘堆积角图像边界并拟合函数求休止角

最新推荐文章于 2023-09-07 22:50:02 发布

别和李俊儒抢辣条

最新推荐文章于 2023-09-07 22:50:02 发布

阅读量6.3k

点赞数 20

文章标签： matlab 图像处理计算机视觉

本文链接：https://blog.csdn.net/qq_35339319/article/details/123407740

版权

休止角（堆积角）

在工程应用中，堆积角（休止角）是一个经常听到的名词，目前测量堆积角普遍使用到堆积角测量仪，通过读取堆积物料高度和堆积物料底基直径对堆积角进行估算，缺点在于最后得到的数据并不是十分的准确，特别是针对于黏性较大的物料。
因此我们可以采用对堆积物料进行拍照，读取边界的形式进行测量，使测量结果更加精确。目前此方法虽已经被广泛使用，但因为其涉及到Matlab，故更多人怕麻烦不愿采用。因此我将一系列工作内容进行一个细致说明，希望帮助到大家。

图像

拍照

如何进行拍照？当拍摄时镜头角度有一个较小偏移时，也会导致最后结果有很大的误差，这时更多的人想到的是固定调平镜头。我想说这样岂不是太麻烦了，在这里我给大家推荐一个方法，采用软件对图像进行垂直水平调整。
我们现在拍摄背景后面放置一张绘有标准矩形的白色打印纸，当然也可以提前将矩形打印在纸上，然后再拍摄，拍摄完成通过软件调试，将图像上的矩形刚好铺满整个画面。该类软件比较多，我就不做过多的赘述。

图像处理

这里我们用到了PhotoShop，怎么处理的我先放几张图片，网上教程太多了，我也不细说，今天我们的任务在于代码。
第一张图片是原始图片，我这里用仿真图片代理代替，第二张是处理完之后的图片，图像大小是804x804像素（可以在PS里面调节，若不是此像素，MATLAB会报错）。

Matlab处理代码

话不多说，上代码！！！
这一段单纯是处理PS里面不想去干的，对图像进行简单的二值化处理并裁剪一半（原因不用多说，如果不裁剪，最后拟合函数就是个二次函数）。

%读取一张图片，并显示（C:\Users\ASUS\Desktop\Matlab-LIJunru\15.jpg这部分内容改为处理完的照片位置）
original_picture=imread('D:\Matlab-LIJunru\15.jpg');%横竖像素点为804
% %把图像转换成灰度图像
% GrayPic=rgb2gray(original_picture);%把RGB图像转化成灰度图像
% figure(2)
% imshow(GrayPic);
% title('RGB图像转化为灰度图像')
% 
%对图像进行二值化处理
figure(1);
subplot(2,4,1);
imshow(original_picture);
title('原始图像')

%thresh=graythresh(original_picture);
thresh = 0.95;
Pic=im2bw(original_picture,thresh);
figure(1);
subplot(2,4,2);
imshow(Pic);
title('RGB图像转化为二值化图像')
imwrite(Pic,'Pic.jpg');
%%
BW1 = Pic;
BW2 = bwperim(BW1,8);
BW3 = imcomplement(BW2);
BW3(:,1) = 1;
BW3(:,804) = 1;
BW3(1,:) = 1;

figure(1);
subplot(2,4,3);
imshow(BW2);
title('边界')
subplot(2,4,4);
imshow(BW3);
title('反相')
imwrite(BW3,'边界.jpg');

BW3(:,[403:804]) = 1;
figure(1);
subplot(2,4,5);
imshow(BW3);
title('取一半')

%BW4=medfilt2(BW3,[804,804]);

完整代码如下，可直接复制粘贴，若有问题可私信我解决:

clc
cla
clear all
close all
%读取一张图片，并显示（C:\Users\ASUS\Desktop\Matlab-LIJunru\15.jpg这部分内容改为处理完的照片位置）
original_picture=imread('D:\Matlab-LIJunru\15.jpg');%横竖像素点为804
% %把图像转换成灰度图像
% GrayPic=rgb2gray(original_picture);%把RGB图像转化成灰度图像
% figure(2)
% imshow(GrayPic);
% title('RGB图像转化为灰度图像')
% 
%对图像进行二值化处理
figure(1);
subplot(2,4,1);
imshow(original_picture);
title('原始图像')

%thresh=graythresh(original_picture);
thresh = 0.95;
Pic=im2bw(original_picture,thresh);
figure(1);
subplot(2,4,2);
imshow(Pic);
title('RGB图像转化为二值化图像')
imwrite(Pic,'Pic.jpg');
%%
BW1 = Pic;
BW2 = bwperim(BW1,8);
BW3 = imcomplement(BW2);
BW3(:,1) = 1;
BW3(:,804) = 1;
BW3(1,:) = 1;

figure(1);
subplot(2,4,3);
imshow(BW2);
title('边界')
subplot(2,4,4);
imshow(BW3);
title('反相')
imwrite(BW3,'边界.jpg');

BW3(:,[403:804]) = 1;
figure(1);
subplot(2,4,5);
imshow(BW3);
title('取一半')

%BW4=medfilt2(BW3,[804,804]);

%%
% figure(2)
% imshow(BW3);
% title('边界反相');
%%
%提取曲线坐标
% 提取图片中的曲线数据
%% 图片与曲线间的定标
im=BW3;
set(0,'defaultfigurecolor','w')
figure(1);
subplot(2,4,6);
 imshow(im)%显示图片
%  title('yuanru1')
[y,x]=find(im==0);%找出图形中的“黑点”的坐标。该坐标是一维数据。
y=max(y)-y;%将屏幕坐标转换为右手系笛卡尔坐标
y=fliplr(y);%fliplr()——左右翻转数组
plot(x,y,'r.','Markersize', 1);
%disp('请在Figrure中先后点击实际坐标框的两个顶点(左上点和右下点)，即A、B两点. ');
%[Xx,Yy]=ginput(2);%Xx,Yy——指实际坐标框的两个顶点
Xx = [0.587467362924286,804.242819843342];
Yy = [801.938502673797,-0.601604278074774];
% min_x=input('最小的x值');%输入x轴最小值
% max_x=input('最大的x值');%输入x轴最大值
% min_y=input('最小的y值');%输入y轴最小值
% max_y=input('最大的y值');%输入y轴最大值
min_x=0;
max_x=100;
min_y=0;
max_y=100;
x=(x-Xx(1))*(max_x-min_x)/(Xx(2)-Xx(1))+min_x;
y=(y-Yy(1))*(min_y-max_y)/(Yy(2)-Yy(1))+max_y;
plot(x,y,'r.','Markersize', 2);
axis([min_x,max_x,min_y,max_y])%根据输入设置坐标范围
title('由原图片得到的未处理散点图')
%% 将散点转换为可用的曲线
%需处理的问题与解决思路
%(1)散点图中可能一个x对应好几个y <---> 保留mean()-std()到mean()+std()之间的y值 并取平均处理
%(2)曲线的最前端和最后段干扰较大 <---> 去掉曲线整体的前(如5%)和后5%
%(3)曲线的最顶端和最底段干扰较大 <---> 去掉曲线整体的上10%和下10%

%参数预设
rate_x=0.08;%曲线的最前端和最后段删除比例
rate_y=0.05;%曲线的最顶端和最底段删除比例

[x_uni,index_x_uni]=unique(x);%找出有多少个不同的x坐标

                x_uni(1:floor(length(x_uni)*rate_x))=[];%除去前rate_x(如5%)的x坐标
                x_uni(floor(length(x_uni)*(1-rate_x)):end)=[];%除去后rate_x的x坐标
                index_x_uni(1:floor(length(index_x_uni)*rate_x))=[];%除去前rate_x的x坐标
                index_x_uni(floor(length(index_x_uni)*(1-rate_x)):end)=[];%除去后rate_x的x坐标

[mxu,~]=size(x_uni);
[mx,~]=size(x);
for ii=1:mxu
    if ii==mxu
        ytemp=y(index_x_uni(ii):mx);
    else
        ytemp=y(index_x_uni(ii):index_x_uni(ii+1));
    end
    %删除方差过大的异常点
                threshold1 = mean(ytemp) - std(ytemp);
                threshold2 = mean(ytemp) + std(ytemp);
                ytemp(find(ytemp<threshold1)) = [];%删除同一个x对应的一段y中的异常点
                ytemp(find(ytemp>threshold2)) = [];
                %删除距顶端和底端较近的点
                thresholdy = (max_y-min_y) * rate_y;%y坐标向阈值
                ytemp(find(ytemp>max_y-thresholdy)) = [];%删除y轴向距离顶端与底端距离小于rate_y的坐标
                ytemp(find(ytemp<min_y+thresholdy)) = [];
    %剩下的y求均值
    y_uni(ii)=mean(ytemp);
end
%此时很多x_uni点处对应的y_uni为空,即NAN,要进一步删去这些空点
            x_uni(find(isnan(y_uni)))=[];
            y_uni(find(isnan(y_uni)))=[];
%画图
figure(1);
subplot(2,4,7);
plot(x_uni,y_uni),title('经处理后得到的扫描曲线')
axis([min_x,max_x,min_y,max_y])%设置坐标范围
Picxx = plot(x_uni,y_uni);
figure(2);
plot(x_uni,y_uni)
print(gcf,'-dpng','扫描曲线.png')
% 将最终提取到的x与y数据保存
            % curve_val(1,:)=x_uni';
            % curve_val(2,:)=y_uni;
%% 对提取出的数据进行拟合(按实际情况进行修改)
[p,s]=polyfit(x_uni,y_uni',1);%多项式拟合(为避免龙格库塔,多项式拟合阶数不宜太高)
[y_fit,DELTA]=polyval(p,x_uni,s);%求拟合后多项式在x_uni对应的y_fit值
figure(1);
subplot(2,4,8);
plot(x_uni,y_fit),title('拟合后的曲线')
axis([min_x,max_x,min_y,max_y])%根据输入设置坐标范围
a = roundn(p(1),-3);
b = roundn(p(2),-3);

figure(3);
plot(x_uni,y_fit)
print(gcf,'-dpng','拟合曲线.png')

disp=(['y = ' num2str(a) 'x + ' num2str(b)])