声明，本文参考Matlab官方文档
1.最大卡尺直径仿真
% 一、得到每一个像素格的圆心坐标、像素正方形框顶点坐标
% 1：函数find()非零区域（1），返回二维数组[y,x]，（大小为21x2像素）即为 每个白色像素点"中心"的坐标
% 2：用cat构建每一个像素框的四个顶点坐标，得到一个三维数组（21*2*4）
% 3：因为最终要得到的是 像素框顶点的坐标，是一个二维数组，所以用reshape（）进行矩阵转换
% 21*4*2=82*2的矩阵；
% 4.用unique消除重复行元素的数组，这样，每一个像素框的顶点坐标就能得到
%
bw = [ ...
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 1 0
0 0 0 0 0 1 1 1 1 1 0
0 0 0 0 1 1 1 1 1 0 0
0 0 0 1 1 1 1 1 1 0 0
0 0 1 1 1 0 0 0 0 0 0
0 0 1 1 1 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 ];
% bw=[...
% 0 0 0 0 0 0 0 0 0 0 0
% 0 0 0 0 0 0 0 0 0 1 0
% 0 0 0 0 0 0 1 1 1 1 0
% 0 0 0 0 0 1 1 1 1 0 0
% 0 0 0 1 1 1 1 1 1 0 0
% 0 0 1 1 1 0 0 0 0 0 0
% 0 0 1 1 1 0 0 0 0 0 0
% 0 0 0 0 0 0 0 0 0 0 0
% 0 0 0 0 0 0 0 0 0 0 0
%
% ];
imshow(bw,'InitialMagnification','fit')%显示像素图，填充满视图框否则太小
pixelgrid
axis on
xlabel('x(个像素)','FontSize',24)
ylabel('y(个像素)','FontSize',24)
[y,x] = find(bw);%找找出BW矩阵中非零元素所在行和列，并存在y，x中；
hold on%保留bw图，继续用plot在bw上作图
plot(x,y,'xb','MarkerSize',5)%标记处像素重心点（用蓝色的的x符号，大小为5）
corners = [x y] + cat(3,[.5 .5],[.5 -.5],[-.5 .5],[-.5 -.5]);
%cat函数用百于连接两个矩阵或数组，第一度个参数表示按第几维进行连接，1表示第一维
%，即行向，也即垂直方向；2表示第二维，及列向，也即水平方向
corners1 = permute(corners,[1 3 2]);%数组维度交换，从x*y*z变为 x*z*y
corners2 = reshape(corners1,[],2);
corners3 = unique(corners2,'rows');%返回corners值不同的行组成的矩阵
plot(corners3(:,1),corners3(:,2),'sr','MarkerSize',5) %显示每个像素的中心和边角
hold off
% reshape函数重新调整矩阵的行数、列数、维数。
% reshape函数总是将原矩阵A，重组为新矩阵B，这里A、B元素个数需相同。重组的规则如下
% 数据排列规则：对A逐列扫描，对B逐列填充，也就是先处理完第一列，在处理第二列，再第三列。。。
% 先处理低维的，再处理高维的，比如要把4*6的A变为6*4的B，就要先扫描A的第一列，
% 扫描过程中行数不断发生变化，列数隔一段时间变化一次，这就是前面说的：先处理低维
% 再处理高维（行是低维，列比行高一维）
% 因此，把4*6的A变为4*3*2的C，扫描和赋值所遵循的规则就是：低维坐标先改变，高维坐标后改变
%
hold on
k = convhull(corners3);
%求 像素框顶点坐标构成的二维矩阵，所形成的的不规则区域的
%凸点“行索引”，返回值K为原矩阵“行索引”而不是x值，逆
%行索引对应的值为凸点坐标
hull_corners = corners3(k,:);%依据行索引K，得到凸点坐标构成的矩阵
plot(hull_corners(:,1),hull_corners(:,2),'r','LineWidth',3)
%依据 列索引（先x坐标值，后y坐标值，产生一个点，依次连接），绘制各个凸点
plot(hull_corners(:,1),hull_corners(:,2),'ro','MarkerSize',10,'MarkerFaceColor','y')
%每个凸点 着重显示 ro 红色 o型 大小为10号
hold off
dz=hull_corners(:,1)'
dx = hull_corners(:,1) - hull_corners(:,1)';
%14个凸点横坐标x，与彼此相减，得到差值（x轴直角边），共14*14个值
dy = hull_corners(:,2) - hull_corners(:,2)';
%14个凸点纵坐标y，与彼此相减，得到差值（y轴直角边），共14*14个值
pairwise_dist = hypot(dx,dy);
%平方和的平方根（类似直角三角形求斜边长)，得到14*14个“斜边长”矩阵
[max_dist,j] = max(pairwise_dist(:));
% [Y,U]=max(A)：返回行向量Y和U，Y向量记录A的每列的最大值，U向量记录每列最大值的行号
% 这里pariwise_dist（：）从14*14矩阵转化为 列矩阵 14*1
%这里返回 最长“斜边”的大小（10），以及行号j=（68）
[k1,k2] = ind2sub(size(pairwise_dist),j);
%这里返回值是 斜边长度矩阵 ，索引号为 j=68 的行号K1=12,列号K2=5
%因为求
%%% !!!!! 斜边矩阵 ：1.列化求最大值得行号；
%%% 2，行号当做 原斜边 矩阵的索引，求出对应 行列号，
%sub2ind函数是MATLAB中对矩阵索引号检索的函数，matlab按列扫描矩阵
% eg:
%
% A = [1 2 3; 4 5 6; 7 8 9]
%
% [I, J] = ind2sub(size(A), [1,7,9])
% 4
% 意思就是把A中索引值为1,7,9的元素（本例中这些元素是1,3和9）的下标输出出来
% 。结果为：横坐标 I = 1 1 3 ，纵坐标J = 1 3 3
point1 = hull_corners(k1,:);
point2 = hull_corners(k2,:);
hold on
plot([point1(1) point2(1)],[point1(2) point2(2)],'-db','LineWidth',3,'MarkerFaceColor','b')
hold off2.最小卡尺直径feretMin
% hull = [
% 2.5000 5.5000
% 3.5000 4.5000
% 6.5000 2.5000
% 9.5000 1.5000
% 10.5000 1.5000
% 10.5000 3.5000
% 9.5000 5.5000
% 5.5000 7.5000
% 2.5000 7.5000
% 2.5000 5.5000 ];
hull = [
2.5000 5.5000
3.5000 4.5000
4.5000 3.5000
5.5000 2.5000
8.5000 1.5000
9.5000 1.5000
10.5000 1.5000
10.5000 2.5000
10.5000 3.5000
9.5000 5.5000
5.5000 7.5000
4.5000 7.5000
3.5000 7.5000
2.5000 7.5000
2.5000 6.5000
2.5000 5.5000
];
% [out,LM] = bwferet(hull,'MinFeretProperties');%计算图像的最小feret直径
% axis on
% maxLabel = max(LM(:));%计算图片bw2的最多个数
%
% axis on
% h = imshow(LM,[]),title('最小弗雷特直径测量');
% axis0 = h.Parent;
% for labelvalues = 1:maxLabel
% xmin = [out.MinCoordinates{labelvalues}(1,1) out.MinCoordinates{labelvalues}(2,1)];
% ymin = [out.MinCoordinates{labelvalues}(1,2) out.MinCoordinates{labelvalues}(2,2)];
% imdistline(axis0,xmin,ymin); %在图形中显示标尺测量
% end
% colorbar('Ticks',1:maxLabel)
% axis off
plot(hull(:,1),hull(:,2),'r','LineWidth',2)
axis equal
axis ij
axis([0 15 0 10])
hold on
xlabel('x(个像素)','FontSize',24)
ylabel('y(个像素)','FontSize',24)
plot(hull(:,1),hull(:,2),'r*')
[x1,y1] = fullLine(gca,[9.5 5.5],-30);
[x2,y2] = fullLine(gca,[5.5 2.5],-30);
caliper_lines(1) = plot(x1,y1,'k');
caliper_lines(2) = plot(x2,y2,'k');
hold off
set(caliper_lines,'Color',[0.8 0.8 0.8]);
[x3,y3] = fullLine(gca,[9.5 5.5],-27);
[x4,y4] = fullLine(gca,[5.5 2.5],-27);
hold on
plot(x3,y3,'k')
plot(x4,y4,'k')
hold off
delete(caliper_lines);
hold on
plot([5.5 7.8],[2.5 6.4],'b','LineWidth',2)
plot([5.5 9.5 5.5 5.5],[7.5 5.5 2.5 7.5],'Color','g','LineWidth',1)
hold off
3.完整程序，注意，先下载测试图像，或者用自己的不规则颗粒图像也可以。原图像如下：使用时导入图像，记得改图像名字
这
3.1未加噪声的不规则颗粒测量
%---------------------------------regionpros()参数描述-------------------------------------%
%'Centroid'：是1行ndims(L)列的向量，给出每个区域的质心（重心）。
%注意：Centroid 的第一个元素是重心水平坐标（x坐标）、第二个元素是重心垂直坐标（y坐标）。
% 'MajorAxisLength'：是标量，与区域具有相同标准二阶中心矩的椭圆的长轴长度（像素意义下
%'MinorAxisLength'：是标量，与区域具有相同标准二阶中心矩的椭圆的短轴长度（像素意义下）
%'BoundingBox'：是1行ndims(L)*2列的向量，即包含相应区域的最小矩形。
%'Area'：是标量，计算出在图像各个区域中像素总个数。
%'Orientation'：是标量，与区域具有相同标准二阶中心矩的椭圆的长轴与x轴的交角（度）。
%-----------------------------------------------------------------------------------------%
clc
clear all
close all
I = imread('PicWaitToProcess.png');
axis on
Img=imnoise(I,'gaussian',0.05,0.005);
M = size(Img);
if numel(M)>2
gray = rgb2gray(Img);
else
gray = Img;
end
% 创建滤波器
W = fspecial('gaussian',[5,5],1);
G = imfilter(gray, W, 'replicate');
abw=imbinarize(G,'global');%图像二值化
bw=~abw;%背景取反
figure(1),
subplot(2,2,1),imshow(Img),title('加噪图像,模拟船舶烟囱坏境')
axis on;
xlabel('x(个像素)','FontSize',20);
ylabel('y(个像素)','FontSize',20);
subplot(2,2,2),imshow(G),title('滤波去噪得到的图像')
axis on;
xlabel('x(个像素)','FontSize',20);
ylabel('y(个像素)','FontSize',20);
subplot(2,2,3),imshow(abw),title('图像二值化')
axis on
xlabel('x(个像素)','FontSize',20);
ylabel('y(个像素)','FontSize',20);
subplot(2,2,4),imshow(bw),title('背景取反便于测量')
axis on
xlabel('x(个像素)','FontSize',20);
ylabel('y(个像素)','FontSize',20);
% I1 = rgb2gray(G); %把图像变为灰度图像
% f = im2double(I1); % 灰度图转double类型--数据处理，保证精度
axis on
axis on
xlabel('x(个像素)','FontSize',20);
ylabel('y(个像素)','FontSize',20);
figure(2),imshow(bw),title('中心坐标标记');
%输入大小为（550*550）,数据类型为double的随机颗粒分布二值图，背景取反
cc = bwconncomp(bw,4);
%-----------------------------------------------------------------------------------------%
% 二值图像中查找所有连通分量（对象）。结果的准确度取决于对象的大小、连通性参数（
% 4、8 或任意值），以及是否有相互接触的对象（在这种情况下，它们可能被标记为一个对象）。
% 二值图像 bw 中的一些颗粒相互接触。
%-----------------------------------------------------------------------------------------%
bw1 = bwareafilt(bw,13);%选取图片中颗粒最大面积的前十个颗粒
bw1pp = bw > 0;
img_reg = regionprops(bw1pp,bw,{'Centroid','WeightedCentroid','BoundingBox'});
numbw1=numel(img_reg);
% 获取img_reg结构体中，颗粒的个数多少
hold on;
for k = 1 : numbw1
plot(img_reg(k).Centroid(1), img_reg(k).Centroid(2), 'bo')
plot(img_reg(k).WeightedCentroid(1), img_reg(k).WeightedCentroid(2), 'r*')
% 利用循环，为颗粒的重心标记上 红色的星号 ，为颗粒的圆心标记上蓝圈
end
hold off
axis on
xlabel('x(个像素)','FontSize',20);
ylabel('y(个像素)','FontSize',20);
figure(3),imshow(bw1),title('对不规则颗粒进行编号')
bw2 = imfill(bw1,'holes');%在提取的对象区域中填充空洞--这里因为随机颗粒
%来自软件生成，噪音少，所以效果不明显
img_reg1 = regionprops(bw1,'BoundingBox');
rects = cat(1, img_reg1.BoundingBox);
centroid = regionprops(bw1,'Centroid');
for i = 1:size(rects, 1)
rectangle('position', rects(i, :), 'EdgeColor', 'r');
text(centroid(i,1).Centroid(1,1)-15,centroid(i,1).Centroid(1,2)-15, num2str(i),'Color', 'r')
% text中的num2str()用来在图片上标记文本，这里通过循环迭代，依据颗粒重心，
% 进行15像素左上角偏移来标记颗粒序号（（重心有x,y），类似Z字扫描
%边界盒产生的矩形边界 rectangle（）；
%边界形状 矩形，边缘颜色 红色
end
[out,LM] = bwferet(bw2,'MinFeretProperties');%计算图像的最小feret直径
axis on
xlabel('x(个像素)','FontSize',20);
ylabel('y(个像素)','FontSize',20);
maxLabel = max(LM(:));%计算图片bw2的最多个数
axis on
xlabel('x(个像素)','FontSize',20);
ylabel('y(个像素)','FontSize',20);
figure(4),
h = imshow(LM,[]),title('最小弗雷特直径测量');
axis0 = h.Parent;
for labelvalues = 1:maxLabel
xmin = [out.MinCoordinates{labelvalues}(1,1) out.MinCoordinates{labelvalues}(2,1)];
ymin = [out.MinCoordinates{labelvalues}(1,2) out.MinCoordinates{labelvalues}(2,2)];
imdistline(axis0,xmin,ymin); %在图形中显示标尺测量
end
colorbar('Ticks',1:maxLabel)
axis off
[out1,LM1] = bwferet(bw2,'MaxFeretProperties');%计算图像的最大feret直径
axis on;
xlabel('x(个像素)','FontSize',20);
ylabel('y(个像素)','FontSize',20);
figure(5),h1 = imshow(LM1,[]);title('最大弗雷特直径测量');
axis on;
xlabel('x(个像素)','FontSize',20);
ylabel('y(个像素)','FontSize',20);
axis1 = h1.Parent;
for labelvalues = 1:maxLabel
xmax1 = [out1.MaxCoordinates{labelvalues}(1,1) out1.MaxCoordinates{labelvalues}(2,1)];
ymax1 = [out1.MaxCoordinates{labelvalues}(1,2) out1.MaxCoordinates{labelvalues}(2,2)];
imdistline(axis1,xmax1,ymax1);
end
axis on;
xlabel('x(个像素)','FontSize',20);
ylabel('y(个像素)','FontSize',20);
colorbar('Ticks',1:maxLabel)
axis on
xlabel('x(个像素)','FontSize',20);
ylabel('y(个像素)','FontSize',20);
grain = false(size(bw));
grain(cc.PixelIdxList{10}) = true;
figure(6),imshow(grain),title('选择标号为10的颗粒')
labeled = labelmatrix(cc);
%-----------------------------------------------------------------------------------------%
% 使用 labelmatrix 根据 bwconncomp 的输出创建标签矩阵。请注意，labelmatrix
% 将标签矩阵存储在依对象数量得出的最小数值类中
%-----------------------------------------------------------------------------------------%
RGB_label = label2rgb(labeled,'spring','c','shuffle');
axis on
xlabel('x(个像素)','FontSize',20);
ylabel('y(个像素)','FontSize',20);
% figure(8),imshow(RGB_label),title('label2rgb伪彩色图像--颜色相同的表示（图像连通方面）关联程度大，颗粒轮廓相似度高')
%-----------------------------------------------------------------------------------------%
% 使用 label2rgb 选择颜色图、背景颜色以及标签矩阵中的对象如何映射到颜色图中的颜色。在伪彩色
% 图像中，用于标识标签矩阵中每个对象的标签映射到相关联的颜色图矩阵中的不同颜色。
%-----------------------------------------------------------------------------------------%
%-----------------------------------------------------------------------------------------%
% 连通性是描述区域和边界的重要概念
% 两个像素连通的两个必要条件是：
% 1、两个像素的位置是否相邻
% 2、两个像素的灰度值是否满足特定的相似性准则（或者是否相等）
%-----------------------------------------------------------------------------------------%
% 基于等效粒径过滤颗粒
BW_out = bwpropfilt(bw, 'MajorAxisLength', [100, 3000]);
axis on
xlabel('x(个像素)','FontSize',20);
ylabel('y(个像素)','FontSize',20);
% figure(9),imshow(BW_out)
axis on
xlabel('x(个像素)','FontSize',20);
ylabel('y(个像素)','FontSize',20);
[out3,LM3] = bwferet(BW_out,'MaxFeretProperties');%计算图像的最大feret直径
figure(7),h3 = imshow(LM3,[]),title('筛选最大弗雷特直径范围[65,135]的颗粒');
axis3 = h3.Parent;
for labelvalues = 1:6
xmax3 = [out3.MaxCoordinates{labelvalues}(1,1) out3.MaxCoordinates{labelvalues}(2,1)];
ymax3 = [out3.MaxCoordinates{labelvalues}(1,2) out3.MaxCoordinates{labelvalues}(2,2)];
imdistline(axis3,xmax3,ymax3);
end
axis on
xlabel('x(个像素)','FontSize',20);
ylabel('y(个像素)','FontSize',20);
colorbar('Ticks',1:maxLabel)
% % Get properties.
% properties = regionprops(BW_out, {'EquivDiameter', 'MinorAxisLength'});
3.2加噪声后的不规则颗粒测量。同理，记得导入检测图像（第三节的第一张图）
I = imread('DREAL2.png');
axis on
Img=imnoise(I,'gaussian',0.7,2);
M = size(Img);
if numel(M)>2
gray = rgb2gray(Img);
else
gray = Img;
end
% 创建滤波器
W = fspecial('gaussian',[5,5],1);
G = imfilter(gray, W, 'replicate');
abw=imbinarize(G,'global');%图像二值化
bw=~abw;%背景取反
figure(1),
subplot(2,2,1),imshow(Img),title('加噪图像,模拟船舶烟囱坏境')
axis on;
xlabel('x(个像素)');
ylabel('y(个像素)');
subplot(2,2,2),imshow(G),title('滤波去噪得到的图像')
axis on;
xlabel('x(个像素)');
ylabel('y(个像素)');
subplot(2,2,3),imshow(abw),title('图像二值化')
axis on
xlabel('x(个像素)');
ylabel('y(个像素)');
subplot(2,2,4),imshow(bw),title('背景取反便于测量')
axis on
xlabel('x(个像素)');
ylabel('y(个像素)');
% I1 = rgb2gray(G); %把图像变为灰度图像
% f = im2double(I1); % 灰度图转double类型--数据处理，保证精度
axis on
BW_out = bwpropfilt(bw, 'MajorAxisLength', [65, 135]);
axis on
xlabel('x(个像素)','FontSize',20);
ylabel('y(个像素)','FontSize',20);
% figure(9),imshow(BW_out)
axis on
xlabel('x(个像素)','FontSize',20);
ylabel('y(个像素)','FontSize',20);
[out3,LM3] = bwferet(BW_out,'MaxFeretProperties');%计算图像的最大feret直径
figure(2),h3 = imshow(LM3,[]),
axis3 = h3.Parent;
for labelvalues = 1:6
xmax3 = [out3.MaxCoordinates{labelvalues}(1,1) out3.MaxCoordinates{labelvalues}(2,1)];
ymax3 = [out3.MaxCoordinates{labelvalues}(1,2) out3.MaxCoordinates{labelvalues}(2,2)];
imdistline(axis3,xmax3,ymax3);
end
axis on;
xlabel('x(个像素)','FontSize',20);
ylabel('y(个像素)','FontSize',20);
xlabel('x(个像素)','FontSize',20);
ylabel('y(个像素)','FontSize',20);
img_reg1 = regionprops(BW_out,'BoundingBox');
rects = cat(1, img_reg1.BoundingBox);
centroid = regionprops(BW_out,'Centroid');
for i = 1:6
rectangle('position', rects(i, :), 'EdgeColor', 'r');
text(centroid(i,1).Centroid(1,1)-15,centroid(i,1).Centroid(1,2)-15, num2str(i),'Color', 'r')
% text中的num2str()用来在图片上标记文本，这里通过循环迭代，依据颗粒重心，
% 进行15像素左上角偏移来标记颗粒序号（（重心有x,y），类似Z字扫描
%边界盒产生的矩形边界 rectangle（）；
%边界形状 矩形，边缘颜色 红色
end
zoom on
figure(3),imshow(Img);
axis on;
xlabel('x（灰度级） 亮度范围','FontSize',18);
ylabel('y(个像素) 每一灰度级所对应的像素个数','FontSize',18);
1126
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
