砂石圆形度测量（MATLAB，图像处理）

砂石圆形度测量（MATLAB，图像处理）

一、设计要求

采集砂石图样，建立圆形度和砂石图像面积及周长的理论表达式，使用MATLAB分析采集到的光学信号并获得相关尺寸，计算得到砂石圆形度。

二、原理分析

1、圆形度的定义

圆形度，即颗粒投影周长与和颗粒投影面积相 等的圆的周长之比。
$$C=L/2\sqrt{A\pi }$$

2、图像处理方法——数学形态法

数字图像的基本特征获取需要对图像边缘进行 检测和有序识别，针对边缘检测的精度越高，相应 的图像处理性能就越强。现代数学形态学在参与数字图像处理中，主要 完成的工序任务为对区域边缘进行精确定位，同时 对封闭轮廓进行信息提取，从而实现对于图像分割 的有效处理。通常情况下，数学形态学所采用的多 尺度滤波机制，会通过尺度空间进行平滑处理来完 成运算，其中，现行高斯尺度空间中，主要采用高 斯尺度滤波器的方式，来对图像当中的噪声进行控 制，最终达到图像简化目的。不过在实际的高斯尺 度空间处理中，针对平滑图像的滤波处理可能导致 图像边缘遭到破坏，无法达到原本的精度诉求，因此需要重新进行功能设定，针对其尺度空间中的定 位需求来进行技术优化。

三、设计步骤

1.RGB 图像转化为灰度图像

2.去除噪声处理，中值滤波法medlifilt2

3.将灰度图像转化为二值图像，调用 im2bw函数

4.对二值图像进行优化处理，填充二值图像前景区域内出现的黑色孔洞，调用 imclose 函数消除不能完全包
含在该圆盘内的目标对象

5.确定目标对象并给予标记

6.二值图像面积提取

7.目标区域边缘提取

四、工程代码

%在图像处理中使用较多的是 HSV 颜色空间，它比 RGB 更接近人们对彩色的感知经验。
% 非常直观地表达颜色的色调、鲜艳程度和明暗程度，方便进行颜色的对比。
%在 HSV 颜色空间下，比 BGR 更容易跟踪某种颜色的物体，常用于分割指定颜色的物体
%%
clear;close all;%清除所有内存变量、图形窗口
I=imread('D:\桌面\Circularity measure_lsy\\stone5.jpg');
%将图像文件 1.jpg 的图像像素数据读入矩阵 I

figure;                      %绘制2*2图像
subplot(2,2,1),imshow(I),title('原图像');

I1=rgb2hsv(I);               %转化RGB值为HSV颜色空间
subplot(2,2,2),imshow(I1),title('HSV图像');

s=I1(:,:,2);                 %显示 S 分量
subplot(2,2,3),imshow(s),title('S分量图像');

s1=imcomplement(s);          %求补运算
subplot(2,2,4),imshow(s),title('S分量图像求补');
%% 
guided = rgb2gray(I);        %将图像转换为灰度图作为引导图像
s=imguidedfilter(s,guided);  %不提取H分量，直接将三维矩阵进行边缘保留滤波
bw=im2bw(s,graythresh(s));   %graythresh(s)计算将灰度图像转化为二值图像所需的门限
                             %将灰度图像转化为二值图像
figure;
subplot(2,2,1),imshow(bw),title('S分量的灰度图像转化为二值图像');

bw1=imfill(bw,'holes');      %填充二值图像中的空洞区域
se1=strel('disk',12);        %创建一个盘形结构元素,其中r指定半径
bw1=imerode(bw1,se1);        %利用盘形平滑处理图像
subplot(2,2,2),imshow(bw1),title('腐蚀图像');

bw1=imdilate(bw1,strel('disk',5));       %使用imdilate函数膨胀图像
subplot(2,2,3),imshow(bw1),title('膨胀图像');
%%
figure;
%bw2=imclose(bw1,strel('disk',10));%采用半径为 10 的圆作为结构元素,进行闭合运算
%subplot(2,2,1),imshow(bw2),title('闭合运算');
bw2=imcomplement(bw1);       %求补运算
subplot(2,2,1),imshow(bw2),title('膨胀图像求补');

bw2=imdilate(bw2,strel('disk',5));        %使用imdilate函数膨胀图像
subplot(2,2,2),imshow(bw2),title('膨胀图像');

[x,y]=size(bw2);             
bw2=imcomplement(bw2);       %求补运算
BW=bwperim(bw2,8);           %获取边缘图像，检测目标的边缘跟踪,用于计算周长
subplot(2,2,3),imshow(BW),title('边缘图像');
%%
%计算垂直方向连续的边界像素点数目%
M1=0;
M2=0;
Ny=0;                        %记录垂直方向边界连续像素点的个数
for i=1:x
    for j=1:y
        if(BW(i,j)>0)
        M2=j;
            if((M2-M1)==1)   %判断是否为垂直方向连续的边界像素点
            Ny=Ny+1;
            end
        M1=M2;
        end
    end
end
%计算水平方向连续的边界像素点数目%
M1=0;
M2=0;
Nx=0;                         %记录水平方向连续边界像素点的个数
for j=1:y
    for i=1:x
        if(BW(i,j)>0)
        M2=i;
            if((M2-M1)==1)    %判断是否为水平方向连续的边界像素点
            Nx=Nx+1;
            end
        M1=M2;
        end
    end
end
SN=sum(sum(BW));              %计算周长像素点的总数
Nd=SN-Nx-Ny;                  %计算奇数码的链码数目
L=sqrt(2)*Nd+Nx+Ny;           %计算周长
[l,num]=bwlabel(bw2);
A=sum(sum(bw2));
stats=regionprops(l,'Area'); %计算目标区域面积
A=stats.Area,
Circularity=L/(2*sqrt(A*pi)), %计算圆形度
C1=A*4*pi/power(L,2);
C2C1=1/power(Circularity,2);
%%
%bw2=imcomplement(bw2);
bw3=double(bw2);
Num_of_Stones=con4_2(bw3,x,y)-1,




*pi/power(L,2);
C2C1=1/power(Circularity,2);
%%
%bw2=imcomplement(bw2);
bw3=double(bw2);
Num_of_Stones=con4_2(bw3,x,y)-1,