形状特征提取代码

叶片形状特征提取代码（matlab函数）

Outline：存储叶片的cell数组

1、叶片面积

输入Outline（叶片的cell数组） 输出Area（叶片面积）

function [Area]=Count_area(Outline)
         column=size(Outline,2);         %计算出图像数量
         Area=zeros(1,column);
         for i=1:column
            S=numel(Outline{1,i});       %返回元素数组中的像素数量
            area=sum(sum(Outline{1,i}));
            Area(1,i)=area;
         end
end

2、叶片周长

输入Outline（叶片的cell数组） 输出Length（叶片周长）

function [Length]=Count_length(Outline)
         column=size(Outline,2);      %计算出图像数量
         Length=zeros(1,column);         %创建矩阵存储面积数据
         for i=1:column
             length=sum(sum(bwperim(Outline{1,i},4)==1));  %bwperin只会返回图像边缘的值
             Length(1,i)=length;
         end
end

3、叶片周长最小外接矩阵

输入Outline（叶片的cell数组） 输出Rectangle数组（存放矩阵的长、宽、面积、周长、中心点横坐标、中心点纵坐标）

function [Rectangle]=Count_rectangle(Outline)
   column=size(Outline,2);     %计算图片数量
   Rectangle=zeros(6,column);  
   for i=1:column
       [r c]=find(Outline{1,i}==1);
       [rectx,recty,area,perimeter] = minboundrect(c,r,'a'); 
       for j=2:4
            k(j-1)=sqrt((rectx(1)-rectx(j))^2 + (recty(1)-recty(j))^2);
        end
            k1=sort(k,2);                                   
            rectwide=k1(1);    %最小外接矩阵的款
            rectlong=k1(2);    %最小外接矩阵的长
        
        for m=1:3
             if(k(m)==k1(3))
                 rectcenterx=abs(rectx(m+1)+rectx(1))/2;     %外接矩形的中心点坐标
                 rectcentery=abs(recty(m+1)+recty(1))/2;     
             end
        end           %外接矩阵的中心点坐标
          Rectangle(1,i)=rectlong;  %长
          Rectangle(2,i)=rectwide;  %宽
          Rectangle(3,i)=area;      %面积
          Rectangle(4,i)=perimeter; %周长
          Rectangle(5,i)=rectcenterx; %外接矩阵中心点横坐标
          Rectangle(6,i)=rectcentery; %外接矩阵中心点纵坐标
   end
   
end

上面Count_rectangle函数需要调用下面minboundrect函数

function [rectx,recty,area,perimeter] = minboundrect(x,y,metric)
if (nargin<3) || isempty(metric)                          %narginp判断输入函数变量个数 判断是否为空 为空取1
   metric = 'a';
 elseif ~ischar(metric)
   error 'metric must be a character flag if it is supplied.'
 else
                                                          % check for 'a'按面积计算 or 'p'按照变长计算
   metric = lower(metric(:)');
   ind = strmatch(metric,{'area','perimeter'});
   if isempty(ind)
    error 'metric does not match either ''area'' or ''perimeter'''
   end
   
  metric = metric(1);
end

  x=x(:);                                                   %将矩阵排成一列
  y=y(:);
  n = length(x);
   if n~=length(y)                                           %计算并判断x y的个数
     error 'x and y must be the same sizes'
   end

if n>3
   edges=convhull(x,y);                                   %得到转折点
   x = x(edges);
   y = y(edges);
   nedges = length(x) - 1;
elseif n>1                                                %n大于1小于等于3
  nedges = n;
  x(end+1) = x(1);
  y(end+1) = y(1);
else
  nedges = n;
end

switch nedges
  case 0
    rectx = [];
    recty = [];
    area = [];
    perimeter = [];
    return
  case 1
    rectx = repmat(x,1,5);
    recty = repmat(y,1,5);
    area = 0;
    perimeter = 0;
    return
  case 2
    rectx = x([1 2 2 1 1]);
    recty = y([1 2 2 1 1]);
    area = 0;
    perimeter = 2*sqrt(diff(x).^2 + diff(y).^2);
    return
end

Rmat = @(theta) [cos(theta) sin(theta);-sin(theta) cos(theta)];  %

ind = 1:(length(x)-1);
edgeangles = atan2(y(ind+1) - y(ind),x(ind+1) - x(ind));         %点与x轴正向夹角

edgeangles = unique(mod(edgeangles,pi/2));                       %unique除去重复项  mod相除取余

nang = length(edgeangles);

area = inf;                                                      %inf表示无穷大量
perimeter = inf;
met = inf;
xy = [x,y];
for i = 1:nang
  rot = Rmat(-edgeangles(i));                                    %构建2维方阵
  xyr = xy*rot;
  xymin = min(xyr,[],1);
  xymax = max(xyr,[],1);
  A_i = prod(xymax - xymin);
  P_i = 2*sum(xymax-xymin);
  
  if metric=='a'
    M_i = A_i;
  else
    M_i = P_i;
  end
  if M_i<met
    met = M_i;
    area = A_i;
    perimeter = P_i;
    rect = [xymin;[xymax(1),xymin(2)];xymax;[xymin(1),xymax(2)];xymin];
  
    rect = rect*rot';
  
    rectx = rect(:,1);
     
    recty = rect(:,2);
  end
end
end 

4、叶片的重心坐标

输入Outline（叶片的cell数组） 输出Gravity数组（存放重心横坐标、重心纵坐标、重心距离边界的最小值）

function [Gravity]=Count_gravity(Outline)
    column=size(Outline,2);
    Gravity=zeros(3,column);
    for i=1:column
          [rows,cols] = size(Outline{1,i}); 
           x = ones(rows,1)*[1:cols];
           y = [1:rows]'*ones(1,cols);   
           area = sum(sum(Outline{1,i})); 
           meanx = sum(sum(Outline{1,i}.*x))/area;    %重心横坐标
           meany = sum(sum(Outline{1,i}.*y))/area;    %重心纵坐标
           Gravity(1,i)=meanx;
           Gravity(2,i)=meany;
           Out_bwperim=bwperim(Outline{1,i},4);
           [y,x] = ind2sub([rows,cols],find(Out_bwperim==1));
           length=0;
           location=0;
           for j=1:size(x)
               if j==1
                   length=sqrt((x(j)-meanx)^2+(y(j)-meany)^2);
                   location=j;
               else
                   if length>sqrt((x(j)-meanx)^2+(y(j)-meany)^2)
                       length=sqrt((x(j)-meanx)^2+(y(j)-meany)^2);
                       location=j;
                   end
               end
           end
           Gravity(3,i)=length;
%            figure
%            imshow(Outline{1,i});hold on;
%            plot(meanx,meany,'r+');                               %十字标出重心位置
           
    
    end
end

5、叶片的凸包

输入Outline（叶片的cell数组） 输出Convex数组（存放凸包面积、凸包周长）

function [Convex]=Count_convex_hull(Outline)
 column=size(Outline,2);
 Convex=zeros(2,column);
 for i=1:column
     convex_length=0;
     [row,col]=size(Outline{1,i});
     [y,x] = ind2sub([row,col],find(Outline{1,i}== 1));  %数组中元素索引值转换为该元素在数组中对应的下标
     
     dt = DelaunayTri(x,y);                              %三角剖分
     if(~size(dt,1))
        return;
     end
     [k,av] = convexHull(dt);     %调用convexhull函数计算图像的凸包
                                  %返回值k为在x、y函数中所处的位置
                                  %av返回凸包的面积
%    显示图像的凸包
%       figure
%      imshow(Outline{1,i});hold on;
%       hold on;
%       plot(x(k), y(k), 'r-');  
%       hold off;
     num=size(k,1);
     k(num+1,1)=k(1,1);
     for j=1:num
         convex_length=sqrt((x(k(j,1))-x(k(j+1,1)))^2+(y(k(j,1))-y(k(j+1,1)))^2)+convex_length;
     end 
     Convex(1,i)=av;              %存储凸包的面积
     Convex(2,i)=convex_length;   %存储凸包的周长   
 end
end

6、叶片长短轴

输入Outline（叶片的cell数组） 输出Shaft数组（存放长轴、短轴）

function [Shaft]=Count_shaft(Outline)
    column=size(Outline,2);
    Shaft=zeros(2,column);
    for i=1:column
        length_long=0;
        length_short=0;
        X1=0;           %存放长短轴的坐标点
        X2=0;
        X3=0;
        X4=0;
        [row,col]=size(Outline{1,i});
        Out_bwperim=bwperim(Outline{1,i},4);
        [y,x] = ind2sub([row,col],find(Out_bwperim==1));  %数组中元素索引值转换为该元素在数组中对应的下标
        for m=1:size(x,1)
            for n=1:size(y,1)
                if(m==1&&n==1)
                    length_long=sqrt((x(m,1)-x(n,1))^2+(y(m,1)-y(n,1))^2);
                     X1=m;
                     X2=n;
                else
                    length_1=sqrt((x(m,1)-x(n,1))^2+(y(m,1)-y(n,1))^2);
                    if(length_1>length_long)
                      length_long=length_1;
                      X1=m;
                      X2=n;  
                    end
                end
             
            end
        end
        k=(y(X2)-y(X1))/(x(X2)-x(X1));  %长轴的斜率
        for m=1:size(x,1)
            for n=1:size(y,1)
               k1=(y(m)-y(n))/(x(m)-x(n)); %短轴的斜率
               min=abs(-(1/k1)-k);
               if min<0.8
                   if(m==1&&n==1)
                    length_short=sqrt((x(m,1)-x(n,1))^2+(y(m,1)-y(n,1))^2);
                     X3=m;
                     X4=n;
                   else
                    length_1=sqrt((x(m,1)-x(n,1))^2+(y(m,1)-y(n,1))^2);
                    if(length_1>length_short)
                      length_short=length_1;
                      X3=m;
                      X4=n;  
                    end
                  end
               end
            end
        end
        Shaft(1,i)=length_long;
        if length_short==0
            length_short=length_long;
        end
        Shaft(2,i)=length_short;
%         figure
%         imshow(Outline{1,i})
%         hold on;
%         plot(x(X1),y(X1),'rx')
%         plot(x(X2),y(X2),'rx')
%         plot([x(X1) x(X2)],[y(X1) y(X2)],'b-')
%         plot(x(X3),y(X3),'rx')
%         plot(x(X4),y(X4),'rx')
%         plot([x(X3) x(X4)],[y(X3) y(X4)],'b-')
        clear Out_bwperim;
    end
   
end

附：如需完整的特征获取的代码，可查看博客下载处/以上函数可直接调用