bw=edge(double(gimg),'canny');
bw(1:10,:)=0; %去除上边界影响
bw(m-10:m,:)=0; %去除下边界影响
%% 形态学闭运算，区域连通
se=strel('line',6,0);
bw=imclose(bw,se);
se=strel('disk',8);
bw=bwareaopen(bw,2000); %去除面积小于2000的区域
bw=imclose(bw,se); %换结构元素填充区域
se=strel('disk',8);
bw=imopen(bw,se); %去除甘蔗之外区域干扰，采用开运算
bw=bwareaopen(bw,5000); %去除面积小于5000的区域
[L,num]=bwlabel(bw); %连通区域标记
%% 选择出剩下连通区域中面积最大的区域，去除其他干扰
tL=L.*0;
aera=zeros(num,1);
for i=1:num
tL=L.*0;
tL(L==i)=1;
aera(i)=sum(tL(:));
end
ind=find(aera==max(aera));
lbw=L.*0;
lbw(L==ind)=1;
se=strel('disk',8);
lbw=imdilate(lbw,se); %甘蔗区域边缘，及孔洞填充
lbw=imfill(lbw); %填充孔洞
%% 搜索甘蔗区域左侧边界，右侧边界
pos_l=zeros(1,2);
pos_r=zeros(1,2);
for i=1:m
tim=lbw(i,:);
ind=find(tim>0);
if ~isempty(ind)%% 判别ind是否大于0
pos_l=[pos_l;i,min(ind)];
pos_r=[pos_r;i,max(ind)];
end
end
pos_l(1,:)=[];
pos_r(1,:)=[];
%% 显示边界
figure(3)
imshow(lbw,[]);
hold on
plot(pos_l(:,2),pos_l(:,1),'g');
plot(pos_r(:,2),pos_r(:,1),'g');
%% 设置步长L，依次计算左边界上相邻三点，组成三角形的面积
Ls=100;
L_angle=zeros(size(pos_l,1),1)-1;
for i=Ls+1:size(pos_l,1)-Ls
pos1=pos_l(i-Ls,:)-pos_l(i,:);
pos2=pos_l(i+Ls,:)-pos_l(i,:);
H=[pos_l(i-Ls,:) 1;pos_l(i,:) 1;pos_l(i+Ls,:) 1];
L_angle(i)=det(H);
if abs(pos1(2))<50 && abs(pos1(2))<50
L_angle(i)=0;
end
end
%% 设置步长L，依次计算右边界上相邻三点，组成三角形的面积
R_angle=zeros(size(pos_r,1),1)-1;
for i=Ls+1:size(pos_r,1)-Ls
pos1=pos_r(i-Ls,:)-pos_r(i,:);
pos2=pos_r(i+Ls,:)-pos_r(i,:);
H=[pos_r(i-Ls,:) 1;pos_r(i,:) 1;pos_r(i+Ls,:) 1];
R_angle(i)=det(H);
if abs(pos1(2))<50 && abs(pos1(2))<50
R_angle(i)=0;
end
end
%% 搜索面积峰值点，即为突出点
figure(4)
plot(L_angle)
indexL=(1:length(L_angle));
indexR=(1:length(R_angle));
[pksL,locL]=findpeaks(L_angle,indexL,'MinPeakHeight',Ls*Ls*1/2,'MinPeakDistance',Ls*5);
[pksR,locR]=findpeaks(R_angle,indexR,'MinPeakHeight',Ls*Ls*1/2,'MinPeakDistance',Ls*5);
% 显示结果
figure(3)
plot(pos_l(locL,2),pos_l(locL,1),'r*');
plot(pos_r(locR,2),pos_r(locR,1),'r*');
%% 计算矩形区域顶点坐标
con_posL=[pos_l(1,:),pos_r(1,:);pos_l(locL,:),pos_r(locL,:);pos_l(end,:),pos_r(end,:)];
con_posR=[pos_l(1,:),pos_r(1,:);pos_l(locR,:),pos_r(locR,:);pos_l(end,:),pos_r(end,:)];
con_pos=[con_posL;con_posR];
ucon_pos=unique(con_pos,'rows');
[~,ind]=sort(ucon_pos(:,1));
sucon_pos=ucon_pos(ind,:);
%% 画矩形框
figure(5)
imshow(im)
hold on
plot(pos_l(locL,2),pos_l(locL,1),'r*');
plot(pos_r(locR,2),pos_r(locR,1),'r*');
for i=1:size(sucon_pos,1)-1
X=[sucon_pos(i,1),sucon_pos(i,3),sucon_pos(i+1,3),sucon_pos(i+1,1),sucon_pos(i,1)];
Y=[sucon_pos(i,2),sucon_pos(i,4),sucon_pos(i+1,4),sucon_pos(i+1,2),sucon_pos(i,2)];
plot(Y,X,'-r');
hold on
X1=mean([sucon_pos(i,1),sucon_pos(i+1,1)]);
Y1=mean([sucon_pos(i,2),sucon_pos(i+1,2)]);
X2=mean([sucon_pos(i,3),sucon_pos(i+1,3)]);
Y2=mean([sucon_pos(i,4),sucon_pos(i+1,4)]);
plot([Y1,Y2],[X1,X2],'k','linewidth',2.0);
end
