采用matlab提取目标形状特征组合

发表于 2024-05-08 分类于 MATLAB ，形状特征，特征提取，形选识别，图像处理阅读次数：本文字数： 19k 阅读时长 ≈ 17 分钟

如下表所示，提取此27组形状特征：

特征提取代码如下：

clc 
clear all;
% close all;
I = imread('F:\笔记本电脑备份\E盘\桌面文件备份重要\国际选煤\形选识别\实验验证数据（10张图像）\07-12-38-481-IMG_High_16.tif');
I=I(:,:,1);
I=I(1:800,100:1146);
I_k= imread('F:\笔记本电脑备份\E盘\桌面文件备份重要\国际选煤\形选识别\实验验证数据（10张图像）\1.tif');
I_k=I_k(:,:,1);
% I_k= I_k(50:200,100:500);

figure(1),subplot(1,2,1),imshow(I);
subplot(1,2,2),imshow(I_k);

%采用中值滤波去除条纹状缺陷
I= medfilt2(I,[3,3]); 

%计算分割阈值
min_I_k = min(min(I_k));

%基于最小空皮带灰度值的图像分割
[m,n] = size(I);
img_I = I;
for i = 1 : m
    for j = 1 : n
        if img_I (i,j) > min_I_k 
            img_I (i,j) = 0;
        else 
             img_I (i,j) = 1;
        end
    end
end

%删除小面积区域
img_I_d= bwareaopen(img_I,100);
figure(2),subplot(1,2,1),imshow(img_I,[0,1]);
subplot(1,2,2);imshow(img_I_d);
%分离连通域
[DD,FF] = bwlabel(img_I_d,8);
%连通域分离
IMG_D = cell(FF,1);
for c = 1:FF
    DD = bwlabel(img_I_d,8);
    for i = 1:m
        for j = 1:n
             if DD(i,j) ~= c
                 DD(i,j) = 0;  
             end
        end
    end
    IMG_D{c,1} = DD;
end
% figure,imshow( IMG_D{5,1});
%显示连通域标记数量
%在伪色彩图像上进行标号
[L,num]=bwlabel(img_I_d,8);     
%封装函数获取各种参数信息
stats_1 = regionprops(L,'Area','PixelList','Centroid','BoundingBox');
%绘制质心坐标
Cen_1 = cat(1,stats_1.Centroid);
imshow(I);
hold on ;
for i=1:length(stats_1) 
    plot(Cen_1(i,1), Cen_1(i,2), 'r+'); 
end
%绘制最小外接矩形框
Temp_1 = cat(1,stats_1.BoundingBox);
for i=1:length(stats_1) 
    temp = Temp_1(i,:);
    rectangle('position',temp,'edgecolor','g');
end 
%标记外接矩形框起始点坐标
for i = 1:length(stats_1)
    plot(Temp_1(i,1),Temp_1(i,2),'g');
    text(stats_1(i).BoundingBox(1)-10,stats_1(i).BoundingBox(2),num2str(i),'Color','r','FontSize',12);
end
hold on;


%%
%循环特征提取
    imLabel = bwlabel(IMG_D{27,1});
    stats1 = regionprops(imLabel,'Area');
    area = cat(1,stats1.Area);
    index = find(area == max(area));
    img = ismember(imLabel,index);
    figure,imshow(img);
%      figure,imshow(I);
    by5 = bwperim(img);
    figure,imshow(by5);
%最小凸多边形检测
num = 1;
hull = cell(num,1);  %凸多边形轮廓
full_img = cell(num,1);%凸多边形面积区域
for i = 1:num 
    stats = regionprops(bwlabel(img), 'ALL');
    tn = stats.ConvexHull;    
    %提取轮廓的凸多边形上的点
    hull{i,1}=stats.ConvexHull;
    hold on;
     h = plot(tn(:, 1), tn(:, 2), 'r');
    set(h, 'Marker', '.');
end

%显示提取轮廓的凸多边形上的点
figure(4),subplot(1,2,1),plot(hull{1,1}(:,2),hull{1,1}(:,1));

%显示提取轮廓的凸多边形区域
full_img{i,1}= stats.ConvexImage;
subplot(1,2,2),imshow(full_img{1,1});

%凸包图像的扩充
figure,imshow(img);
for i=1:size(stats)
     rectangle('Position',[stats(i).BoundingBox],'LineWidth',2,'LineStyle','--','EdgeColor','r'),
end

img1= img(ceil(stats(i).BoundingBox(1,2)): ceil(stats(i).BoundingBox(1,2))-1+stats(i).BoundingBox(1,4),ceil(stats(i).BoundingBox(1,1)):ceil(stats(i).BoundingBox(1,1))+stats(i).BoundingBox(1,3)-1);  %提取二值化图像最小外界矩形区域
img2 = full_img{1,1}; %凸包二值化图像

img3 = img2-img1;
figure,subplot(1,2,1),imshow(img3),hold on;
img3= bwareaopen(img3,100); 
subplot(1,2,2),imshow(img3);


%计算颗粒的最大外接圆直径
pic=~img;
figure;
imshow(pic); title('原图（彩图以及灰度图二值图都可以）');
%%%% %%%% %%%% ！！！调用函数计算最小外接圆圆心与半径，且默认图片背景为白；因为我测试图片背景为白色！
[minCnum,center,radius] = m_minCircumferentialCircle(pic);
%%% %%%% %%%% 显示 pic，并在其上画出所有求得的最小外接圆
imshow(I);
for i = 1:minCnum   % 利用画线函数，画最小外接圆
    %%%% 使用画线函数在pic图上  画圆
    theta = 0:2*pi/3600:2*pi;
    Circle1 = center(i,2) + radius(i)*cos(theta);
    Circle2 = center(i,1) + radius(i)*sin(theta);% 构建圆
    hold on;
    plot(Circle1,Circle2,'m','Linewidth',1);% 画线
    title('最小外接圆');
    plot(center(1,2),center(1,1),'r+','Linewidth',2);
end
%绘制半径
global elem;
x = [center(1,2),elem(1,2)];
y = [center(1,1),elem(1,1)];
plot(x,y,'g','Linewidth',1);
%最大外界圆直径
R = 2*radius;

%凹点检测之寻找凹陷边缘
by = bwperim(img3);
figure,imshow(by);
by2 = bwperim(img2);
figure,imshow(by2);
by3 = by - by2;
[m,n] = size(by3);
for i = 1:m
    for j = 1:n
        if by3(i,j) == -1;
            by3(i,j) =0;
        end
    end
end
by3 =  bwareaopen(by3,15);%删除小段边缘线
% [oo pp]= bwlabel(by3);
figure,imshow(by3);

%2021.11.24
%关于边缘线需要进行修改
%二值化取反图像连通域检测
%GZ元胞数组存储元原二值图像重合的边缘线，EZ元宝数组存储与凸包边缘重合的边缘线，且GS与EZ对应
[G,num] = bwlabel(by,8);
%连通域分离
GZ = cell(num,1);
for c = 1:num
    G = bwlabel(by,8);
    for i = 1:m
        for j = 1:n
             if G(i,j) ~= c
                 G(i,j) = 0;  
             end
        end
    end
    GZ{c,1} = G .* (1/c);   
    GZ{c,1} = GZ{c,1}-by2;
    for i = 1:m
    for j = 1:n
        if GZ{c,1}(i,j) == -1;
            GZ{c,1}(i,j) =0;
        end
    end
    end
    GZ{c,1} = bwareaopen(GZ{c,1},5);
end
% figure,imshow(GZ{1,1});

 by4 = by - by3;
 figure,imshow( by4);
%分离边缘线连通域
%二值化取反图像连通域检测
[G,num] = bwlabel(by4,8);
%连通域分离
EZ = cell(num,1);
for c = 1:num
    G = bwlabel(by,8);
    for i = 1:m
        for j = 1:n
             if G(i,j) ~= c
                 G(i,j) = 0;  
             end
        end
    end
    EZ{c,1} = G .* (1/c);
    EZ{c,1} = EZ{c,1}-GZ{c,1};
    for i = 1:m
    for j = 1:n
        if EZ{c,1}(i,j) == -1;
            EZ{c,1}(i,j) =0;
        end
    end
    end
    EZ{c,1} = bwareaopen(EZ{c,1},5);
end
% figure,imshow(EZ{2,1});



%以凹陷边缘点坐标为引导，并为中心，计算在二值图像img上30*30区域内像素和
xx = cell(num,1);
yy = cell(num,1);
for c = 1:num
    [xx{c,1},yy{c,1}]=find(GZ{c,1}==1);
end
%计算边缘线连通域坐标点在原二值化图像中30*30区域像素和
he = cell(num,1);%原二值化图像中30*30区域像素和
figure,imshow(img1),hold on;

for rr = 1:num
    [qq,ww] = size(xx{rr,1});   
    for i = 1 : qq
          if (yy{rr,1}(i,1) >5 && xx{rr,1}(i,1) > 5)&&(yy{rr,1}(i,1) <n-5 && xx{rr,1}(i,1)<m-5)

           he{rr,1}(i,1) = sum(sum(img1(xx{rr,1}(i,1)-5 : xx{rr,1}(i,1)+5,yy{rr,1}(i,1)-5:yy{rr,1}(i,1)+5)));
%           figure,imshow(img1(xx{rr,1}(i,1)-5 : xx{rr,1}(i,1)+5,yy{rr,1}(i,1)-5:yy{rr,1}(i,1)+5));
%           cc = [yy{rr,1}(i,1)-5,yy{rr,1}(i,1)+5,yy{rr,1}(i,1)+5,yy{rr,1}(i,1)-5,yy{rr,1}(i,1)-5];
%           bb = [xx{rr,1}(i,1)-5,xx{rr,1}(i,1)-5,xx{rr,1}(i,1)+5,xx{rr,1}(i,1)+5,xx{rr,1}(i,1)-5];
%           plot(  cc ,  bb ,'r','Linewidth',1);
%          plot(yy{1,1}(i,1),xx{1,1}(i,1),'r+','Linewidth',2); %在二值图上绘制重合边缘线
          end
    end
end
figure,
for rr = 1:num
    plot(he{rr,1}),hold on;
end
%计算最大凹陷深度
%从he中返回最大值标号，根据标号对应至xx和yy寻找最大凹陷深度的点坐标
%得到此坐标计算此坐标点到凸包边缘的最小距离   即为最大凹陷深度
K = zeros(rr,1);
for rr = 1:num
       [k,K(rr,1)]= max(he{rr,1});%从he中返回最大值标号
end

figure,imshow(img1),hold on;
for rr = 1:num
    plot(yy{rr,1}(K(rr,1),1),xx{rr,1}(K(rr,1),1),'r+','Linewidth',2);%根据标号对应至xx和yy寻找最大凹陷深度的点坐标
end

%最大凹陷点到凸包的最小距离，找最小距离点坐标并绘制线
figure,imshow(by2),hold on;
for rr = 1:num
    plot(yy{rr,1}(K(rr,1),1),xx{rr,1}(K(rr,1),1),'r+','Linewidth',2);%根据标号对应至xx和yy寻找最大凹陷深度的点坐标
end

%找到最小距离并绘制凸起缺陷

for rr = 1:num
        pp=1;
        for i = 1:m      
            for j = 1:n
                if EZ{rr,1}(i,j) == 1% 边缘像素点
                    elemn = [i j]; % 该像素位置
                    elemn_x{rr,1}(pp,1) = i;
                    elemn_y{rr,1}(pp,1) = j;
                    distance{rr,1}(pp,1)= norm( elemn -[xx{rr,1}(K(rr,1),1),yy{rr,1}(K(rr,1),1)]); % 该像素与最大凹点距离（范数）
                    pp = pp+1;
                end
            end
        end
    ll(rr,1) = min(distance{rr,1});   
end

%最大缺陷距离    
if sum(sum(img3)) == 0 
        MAX = 0;
else
        MAX = max(ll);
end
% figure,
% for rr = 1:num
%     plot(distance{rr,1}),hold on;
% end

%最大缺陷对应凹陷点坐标以及最短距离点坐标
D = zeros(rr,1);
for rr = 1:num
       [d,D(rr,1)]= min(distance{rr,1});%从distance中返回最小值标号
end

figure,imshow(img2),hold on;
for rr = 1:num
    plot(elemn_y{rr,1}(D(rr,1),1),elemn_x{rr,1}(D(rr,1),1),'r+','Linewidth',2);%根据标号对应至xx和yy寻找最大凹陷深度的点坐标
end
%连接最大缺陷对应凹陷点坐标以及最短距离点坐标

figure,imshow(img3),hold on;
for rr = 1:num
    cc = [elemn_y{rr,1}(D(rr,1),1),yy{rr,1}(K(rr,1),1)];
    bb = [elemn_x{rr,1}(D(rr,1),1),xx{rr,1}(K(rr,1),1)];
    plot(  cc ,  bb ,'r','Linewidth',1),hold on;
end


%最大凹陷点坐标
a = cell(num,1);
for rr = 1:num
    plot(yy{rr,1}(K(rr,1),1),xx{rr,1}(K(rr,1),1),'r+','Linewidth',2);%根据标号对应至xx和yy寻找最大凹陷深度的点坐标
    a{rr,1} = [yy{rr,1}(K(rr,1),1),xx{rr,1}(K(rr,1),1)];
end

%获取边缘线连通域坐标点
%边界跟踪
%以距离凹点最远的点为边缘跟踪的初始点
elemn = zeros();
elemn_x = cell(num,1);
elemn_y =cell(num,1);
% figure,imshow(GZ{1,1});

for rr = 1:num
        pp=1;
        [m,n] = size(GZ{rr,1});
        for i = 1:m      
            for j = 1:n
                if GZ{rr,1}(i,j) == 1% 边缘像素点
                    elemn = [i j]; % 该像素位置
                    elemn_x{rr,1}(pp,1) = i;
                    elemn_y{rr,1}(pp,1) = j;
                    distance{rr,1}(pp,1)= norm( elemn -[xx{rr,1}(K(rr,1),1),yy{rr,1}(K(rr,1),1)]); % 该像素与最大凹点距离（范数）
                    pp = pp+1;
                end
            end
        end
    ll(rr,1) = min(distance{rr,1});   
end

%sum(sum(GZ{2,1}))
D = zeros(num,1);
for rr = 1:num
    if length(distance{rr,1})>sum(sum(GZ{rr,1}))
        distance{rr,1} = distance{rr,1}(1:length(elemn_y{rr,1}),:); 
    end
end
for rr = 1:num
       [d,D(rr,1)]= max(distance{rr,1});%从distance中返回最大值标号
end

figure,imshow(img2),hold on;
for rr = 1:num
    plot(elemn_y{rr,1}(D(rr,1),1),elemn_x{rr,1}(D(rr,1),1),'r+','Linewidth',2);%根据标号对应至xx和yy寻找最大凹陷深度的点坐标
end

figure,imshow(by3),hold on;
for rr = 1:num
    cc = [elemn_y{rr,1}(D(rr,1),1),yy{rr,1}(K(rr,1),1)];
    bb = [elemn_x{rr,1}(D(rr,1),1),xx{rr,1}(K(rr,1),1)];
    plot(  cc ,  bb ,'r','Linewidth',1),hold on;
end
for rr = 1:num
    a1{rr,1}=[elemn_y{rr,1}(D(rr,1),1),elemn_x{rr,1}(D(rr,1),1)];%根据标号对应至xx和yy寻找最大凹陷深度的点坐标
end
%a1为边缘搜索的初始点
% figure,imshow(GZ{rr,1});
for rr = 1:num
        dim=size(GZ{rr,1});
        col=a1{rr,1}(1,2); %计算起始点列坐标
        row=a1{rr,1}(1,1);%计算起始点行坐标        
        figure,imshow(img3),hold on;
        plot(row,col,'r+','Linewidth',2);
        connectivity=8;
        num_points=sum(sum(GZ{rr,1}));
        contour{rr,1}=bwtraceboundary(GZ{rr,1},[col,row],'N',connectivity, num_points);
        %提取边界
        plot(contour{rr,1}(:,2),contour{rr,1}(:,1), 'g+','LineWidth' ,1);
        title('边界跟踪图像');
end

% plot(contour(52,2),contour(52,1), 'm+','LineWidth' ,4);
xx = cell(num,1);
yy = cell(num,1);
for c = 1:num
    xx{c,1}=contour{c,1}(:,1);
    yy{c,1}=contour{c,1}(:,2);
end

%计算边缘线连通域坐标点在原二值化图像中30*30区域像素和
he = cell(num,1);%原二值化图像中30*30区域像素和
figure,imshow(img1),hold on;

for rr = 1:num
    [qq,ww] = size(xx{rr,1});   
    for i = 1 : qq
          if (yy{rr,1}(i,1) >5 && xx{rr,1}(i,1) > 5)&&(yy{rr,1}(i,1) <n-5 && xx{rr,1}(i,1)<m-5)

           he{rr,1}(i,1) = sum(sum(img1(xx{rr,1}(i,1)-5 : xx{rr,1}(i,1)+5,yy{rr,1}(i,1)-5:yy{rr,1}(i,1)+5)));
%           figure,imshow(img1(xx{rr,1}(i,1)-5 : xx{rr,1}(i,1)+5,yy{rr,1}(i,1)-5:yy{rr,1}(i,1)+5));
          cc = [yy{rr,1}(i,1)-5,yy{rr,1}(i,1)+5,yy{rr,1}(i,1)+5,yy{rr,1}(i,1)-5,yy{rr,1}(i,1)-5];
          bb = [xx{rr,1}(i,1)-5,xx{rr,1}(i,1)-5,xx{rr,1}(i,1)+5,xx{rr,1}(i,1)+5,xx{rr,1}(i,1)-5];
          plot(  cc ,  bb ,'r','Linewidth',1);
%          plot(yy{1,1}(i,1),xx{1,1}(i,1),'r+','Linewidth',2); %在二值图上绘制重合边缘线
          end
    end
end
figure,
for rr = 1:num
    plot(he{rr,1}),hold on;
end

K = zeros(rr,1);
for rr = 1:num
    if sum(sum(he{rr,1})~=0)
       [k,K(rr,1)]= max(he{rr,1});%从he中返回最大值标号
    end
        if sum(sum(he{rr,1})==0)
       K(rr,1)= 10%从he中返回最大值标号
    end
end

figure,imshow(img1),hold on;
for rr = 1:num
    plot(yy{rr,1}(K(rr,1),1),xx{rr,1}(K(rr,1),1),'r+','Linewidth',2);%根据标号对应至xx和yy寻找最大凹陷深度的点坐标
end

%凹缺陷的个数GS
GS = num;
%最大凹陷角度
a = cell(num,1);
for rr = 1:num
    plot(yy{rr,1}(K(rr,1),1),xx{rr,1}(K(rr,1),1),'r+','Linewidth',2);%根据标号对应至xx和yy寻找最大凹陷深度的点坐标
    a{rr,1} = [yy{rr,1}(K(rr,1),1),xx{rr,1}(K(rr,1),1)];
end


figure,imshow(img1),hold on;
for rr = 1:num

           if (K(rr,1)-5) >0 && (K(rr,1)+5)<length(xx{rr,1})
                plot(yy{rr,1}(K(rr,1)-5,1),xx{rr,1}(K(rr,1)-5,1),'b+','Linewidth',2),hold on;%根据标号对应至xx和yy寻找最大凹陷深度的点坐标
                plot(yy{rr,1}(K(rr,1)+5,1),xx{rr,1}(K(rr,1)+5,1),'g+','Linewidth',2),hold on;
                plot(yy{rr,1}(K(rr,1),1),xx{rr,1}(K(rr,1),1),'r+','Linewidth',2);
           end

end

b = cell(num,1);
c = cell(num,1);
for rr = 1:num
    if (K(rr,1)-5) >0 && (K(rr,1)+5)<length(xx{rr,1})
    b{rr,1} = [yy{rr,1}(K(rr,1)+5,1),xx{rr,1}(K(rr,1)+5,1)];
    c{rr,1} = [yy{rr,1}(K(rr,1)-5,1),xx{rr,1}(K(rr,1)-5,1)];
    end
end

aa = cell(num,1);
bb = cell(num,1);
cc = cell(num,1);
th = zeros(rr,1);
for rr = 1:num
    if (K(rr,1)-5) >0 && (K(rr,1)+5)<length(xx{rr,1})
     aa{rr,1} =norm([b{rr,1}(1,1)-c{rr,1}(1,1),b{rr,1}(1,2)-c{rr,1}(1,2)]);
     bb{rr,1} =norm([a{rr,1}(1,1)-c{rr,1}(1,1),a{rr,1}(1,2)-c{rr,1}(1,2)]);
     cc{rr,1} =norm([c{rr,1}(1,1)-a{rr,1}(1,1),c{rr,1}(1,2)-a{rr,1}(1,2)]);
     acpha{rr,1} =(acos((bb{rr,1}^2+cc{rr,1} ^2-aa{rr,1}^2)/(2*bb{rr,1}*cc{rr,1})))/(2*pi)*360;
     th(rr,1) = acpha{rr,1};
     end
end

if sum(sum(th))==0
    acpha =180;
elseif sum(sum(img3))~=0 
    [jjj,ttt,acpha] = find(th);
    acpha = min(acpha);
end

%计算凹缺陷处的非零像素面积
figure,imshow(img1),hold on;
%计算矩形与边缘相交点的坐标并在图中绘制
for rr = 1:num
          oo = [a{rr,1}(1,1)-5,    a{rr,1}(1,1)-5,   a{rr,1}(1,1)+5,  a{rr,1}(1,1)+5,   a{rr,1}(1,1)-5];
          uu = [a{rr,1}(1,2)-5,    a{rr,1}(1,2)+5,   a{rr,1}(1,2)+5,  a{rr,1}(1,2)-5,   a{rr,1}(1,2)-5];
          plot(  oo ,  uu ,'r','Linewidth',2);   
end

for rr = 1:num
    img_f{rr,1} = img1(a{rr,1}(1,2)-5: a{rr,1}(1,2)+5, a{rr,1}(1,1)-5 : a{rr,1}(1,1)+5);
end
% figure,imshow(img_f{2,1});      
tt = 0;
FL = 0;
for rr = 1:num
    if sum(sum(img3))~=0
    tt(rr,1) = sum(sum(img_f{rr,1}));
    TT(rr,1) = 11*11 - tt(rr,1);
    FL(rr,1) = tt(rr,1)/TT(rr,1);
    end
end
MS = max(tt);
FLB = max(FL);

        %特征统计目前有19个
        S1 = sum(sum(img));
        C1= sum(sum(by5));
        e1= (4*pi*S1)/C1^2;

        S2 = sum(sum(img2));
        C2= sum(sum(by2));
        e2= (4*pi*S2)/C2^2;

        L1 = R;
        S3 = sum(sum(img3));%凹缺陷面积
        C3 = sum(sum(by));
        %获取凹缺陷最大连通域区域图像
        imLabel = bwlabel(img3);
        stats1 = regionprops(imLabel,'Area');
        area = cat(1,stats1.Area);
        index = find(area == max(area));
        img4 = ismember(imLabel,index);
%         figure,imshow(img4);
        by6 = bwperim(img4);
%          figure,imshow(by6);

        S4 = sum(sum(img4));  %最大凹缺陷面积
        C4 = sum(sum(by6));  %最大凹缺陷周长
        L2 = MAX; %最大凹陷深度

        LR = C1/C2;
        SR = S1/S2;
        eR = e1/e2;
        LDR = L2/L1;

%距离变换特征
D =-bwdist(~img);
imshow(D,[]);
mask = imextendedmin(D,0.8);

imshowpair(I,mask,'blend');
[NUM,kkk] = bwlabel(mask);

%特征再次组合
LK1 = S1/C1;
LK2 = S2/C2;
if C3==0
    LK3 =0;
else
    LK3 = S3/C3;
end
if C4==0
    LK4=0;
else
    LK4 = S4/C4;
end

%计算长宽比
stats2 = regionprops(img,'BoundingBox');
% figure,imshow(img_I_d),hold on;
% plot(stats2.BoundingBox(1,1),stats2.BoundingBox(1,2),'r+');
ck = stats2.BoundingBox(1,3)/stats2.BoundingBox(1,4);%最小外接矩形长宽比


%多重灰度阈值分割
[m,n] = size(I);
[G,num] = bwlabel(img,8);
p = [170,160,140,120,100,80,60];
%p = [160,150,140,120,110,90,70,60]; 实验较好
K = img.*double(I);
% figure,imshow(K,[]);
            q = 1;   
            u = 1;
            E = zeros(m,n); 
            %存取多重灰度阈值分割后图像  
                    while u <= 7
                           g = p(1,u)  ;               %分割阈值循环                
                                    for x = 1:m                     %二值化
                                         for y = 1:n
                                            if K(x,y) > g
                                                  E(x,y) =0;
                                            elseif K(x,y) == 0
                                                 E(x,y) =0;
                                            else
                                                E(x,y) =1;
                                           end
                                         end
                                    end  
%                                 [L,H(q,1)] = bwlabel(E,8);           %联通域检测
                                   A{q,1} = E;
                                 
                                  figure(q), imshow(A{q,1});
                                    q = q+1;
                                    u = u+1;
                    end   
                    for i = 1:7
                        A{i,1}=bwareaopen(A{i,1},30);
                    end
                    for i = 1:7
                        [L,H(i,1)] = bwlabel(A{i,1},8); 
                    end
KH = max(H);
AAAA= [S1,C1,e1,LK1, ck, S2 ,C2,e2,LK2,L1,S3,C3,LK3,S4,C4,LK4,L2,LR,SR,eR, LDR,GS,acpha,MS,FLB,kkk,KH];
close all;

调用函数：m_minCircumferentialCircle.m

%An highlighted block
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              图像处理    
%       从图片相对简单、单一照片中获取其中物体最小外接圆,返回圆心、半径

% 1、调用函数：[PIC,center,radius] = m_minCircumferentialCircle(pic,background)
%                       采用遍历查询确定最小外接圆半径。

%    pic —— 原图（彩图/灰度图/二值图）
%    background —— 原图背景色，0则为黑、1为白；
%                  默认背景为白色，pic的二值图中黑区表示的物体
%    minCnum —— 外接圆个数
%    center —— 圆中心位置，数据格式为： [height1,width1]
%                                     [height2,width2]
%    radius —— 圆半径，数据格式为：[radius1]
%                                [radius2]
% 注：该函数没有边界检测，所画的圆可能超出原图大小。
%     该函数的思想是：1）、对输入图片进行二值化；
%                    2）、利用bwlabel()函数判断其二值图中的连通区域（这里前提是以每一个连通区域作为一个物体）；
 %                    3）、对每个连通区域求其质心，以质心作为其最小外接圆的圆心；
%                     4）、遍历各个像素点，分别计算每个物体的像素点是否在其物体当前的半径内；如果不在则用当前像素点与圆心的距离作为新的半径；从而类似迭代更新得到最小外接圆半径。
%   笔者采用的是2018a版的MatLAB,其中例如二值化函数可能与之前版本不一致                  
% https://blog.csdn.net/BinHeon
% 2019/5/24
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [minCnum,center,radius,elem] = m_minCircumferentialCircle(pic,background)
%%%%  %%%%  %%%%  %%%%  输入参数处理
% if nargin < 2      % 如果只输入pic图片，默认背景为白
%     background = 1;
% end
%%%%  %%%%  %%%%  %%%%  处理图片
% if (ndims(pic) == 3) % 获取输入图片维度，如果是三维图像需要获得其二值图
%     bw1 = rgb2gray(pic);
%     spic = imbinarize(bw1); % 2值化，之后的便是对二值图进行开闭运算
% else
%     spic = imbinarize(pic); % 如果输入图片为灰度图，则直接二值化
% end
% figure;
% imshow(spic);
spic = pic;
%%%%  %%%%  %%%%  %%% 获取二值图像中的连通区域
%%%%  注意：此bwlabel函数标注的是二值图中值为“1”的点，因而需要确定其二值图背景与物体
%%%%        如果背景为白色，则需要对二值图像取反操作。
% if background == 1
    [L,num]=bwlabel(~spic,8);    % 标注二值图中已连接的黑色（0），即背景为白，物体为黑
% else 
%     [L,num]=bwlabel(spic,8);     % 标注二值图中已连接的白色（1），即背景为黑，物体为白
% end

%%%%  %%%%  %%%%  %%% 获取最小外接圆个数
minCnum = num;

%%%%  %%%%  %%%%  %%%%  利用二值图获取物体质心位置
[height,width]=size(spic); % 获取图片大小
%%%%  %%%% 初始化圆心
global center;
center =zeros(minCnum,2);   % 用于记录质心位置的坐标
%%%%%%%% 计算每个物体中心
for cn = 1:minCnum
    % 初始化参数
    sum_x=0;sum_y=0;area=0;
    for i=1:height
        for j=1:width
            if L(i,j) == cn
                sum_x=sum_x+i;
                sum_y=sum_y+j;
                area=area+1;
            end
        end
    end
    %%质心坐标
    center(cn,1) = fix(sum_x/area);
    center(cn,2) = fix(sum_y/area);   
end

%%%%  %%%% %%%% %%%% 遍历每一个像素点求每个物体的半径
global radius;
global elem;
for rc = 1:minCnum
	radius = zeros(minCnum,1);% 初始化一个minCnum*1的矩阵，用于存每个半径值
    radius(:,1) = 2; % 每个最小外接圆，初始半径2个像素点
    for i = 1:height       % 遍历图片，查找物体像素点（二值化之后，物体像素点“0”）
        for j = 1:width
            if L(i,j) == rc % 物体像素点
                elem = [i j]; % 该像素位置
                distance = norm( elem - center(rc,:) ); % 该像素与质心距离（范数）
                if distance > radius(rc,1)
                    radius(rc,1) = round(distance);
                    elem ;
                end
            end
        end
    end
end

最终输出的特征组合：

上述运行代码输出结果为独立目标的，粘连目标输出结果为：

特征提取可视化过程描述详见博士论文。

最终组建的部分训练集：