function Trace=BoundaryDetect(I)
% author:zll
% date:2018_7_7
%代码还有优化空间
J=im2uint16(I);
J=(J>30000)*10000;
[m,n]=size(J);
% bwlabel
% first step:
% 标记类栈矩阵
% 先腐蚀再膨胀,去除毛边(开运算)
J=MyCorruption(J);
J=Mydilate(J);
stackMatrix = [];
index=1;
for index1=1:m
for index2=1:n
if J(index1,index2)==10000
index3=index1; %初始化栈
index4=index2;
stackMatrix(:,end+1)=[index3;index4];
[c v]=size(stackMatrix);
while c~=0 && v~=0
%出栈
%只在出栈的时候才赋值
index3=stackMatrix(1,end);
index4=stackMatrix(2,end);
J(index3,index4)=index;
stackMatrix(:,end)=[];
% 4种情况讨论,找到标记值入栈
if index3>1
if J(index3-1,index4)==10000
%入栈
stackMatrix(1,end+1)=index3-1;
stackMatrix(2,end)=index4;
end
end
if index4>1
if J(index3,index4-1)==10000
stackMatrix(1,end+1)=index3;
stackMatrix(2,end)=index4-1;
end
end
if index3<m
if J(index3+1,index4)==10000
%入栈
stackMatrix(1,end+1)=index3+1;
stackMatrix(2,end)=index4;
end
end
if index4<n
if J(index3,index4+1)==10000
%入栈
stackMatrix(1,end+1)=index3;
stackMatrix(2,end)=index4+1;
end
end
[c v]=size(stackMatrix);
end
index=index+1;
end
end
end
% 去噪声
% 连通域面积太小的舍去
statistics=zeros(1,index);
for index3=1:m
for index4=1:n
if J(index3,index4)~=0
statistics(1,J(index3,index4))=statistics(1,J(index3,index4))+1;
end
end
end
for index3=1:index
if statistics(1,index3)<100
statistics(1,index3)=0;
end
end
for index3=1:m
for index4=1:n
if J(index3,index4)~=0 && statistics(1,J(index3,index4))==0
J(index3,index4)=0;
end
end
end
% 为提高效率,先找粗略找一遍边界像素(本例此操作无意义,单就粗犷的轮廓已经可以描绘)
% 边界时min(index3+1,256)或max(index3-1,1)一定在数组定义范围内,且取值为判定区域内点,不影响判断结果
% F=zeros(m,n);
% for index3=1:m
% for index4=1:n
% if J(index3,index4)~=0 && F(index3,index4)==0
% if J(max(index3-1,1),max(index4-1,1))==0 || J(max(index3-1,1),index4)==0 || J(max(index3-1,1),min(index4+1,256))==0 || J(index3,max(index4-1,1))==0 || J(index3,min(index4+1,256))==0 || J(min(index3+1,256),max(index4-1,1))==0 || J(min(index3+1,256),index4)==0 || J(min(index3+1,256),min(index4+1,256))==0
% F(index3,index4)=J(index3,index4);
% end
% end
% end
% end
% 开始统计边界坐标,此种方法在米粒的最上端没法保障单像素;因为最开始一个元素需要走两遍。
% flag=1;
% 缺点,在遇到右边界只露半只米的情况可能会损失一个像素的信息
% Trace数组只需要起始点和后续每步的方向值即可画出连通域轨迹,即bwboundaries的效果。
% LabelMatrix每次使用完清0。
Q=J;
% direction=[7 0 1;6 8 2;5 4 3];
Trace=[];
LabelMatrix=[];
num1=0;
for index3=1:m
for index4=1:n
if Q(index3,index4)~=0
index1=index3;
index2=index4;
flag=Q(index3,index4);
%初始化
if Q(index1,min(index2+1,256))==flag && index2+1~=257
LabelMatrix(1:3,1)=[index1; index2; 2];
index2=index2+1;
elseif Q(min(index1+1,256),min(index2+1,256))==flag && index2+1~=257 && index1+1~=257
LabelMatrix(1:3,1)=[index1; index2; 3];
index1=index1+1;
index2=index2+1;
elseif Q(min(index1+1,256),index2)==flag && index1+1~=257
LabelMatrix(1:3,1)=[index1; index2; 4];
index1=index1+1;
elseif Q(min(index1+1,256),max(index2-1,1))==flag && index1+1~=257 && index2-1~=0
LabelMatrix(1:3,1)=[index1; index2; 5];
index1=index1+1;
index2=index2-1;
end
%temp=0;
%去过噪声,故此while语句可行
%需要进入while语句,而且第一次执行判断满足条件
%while index1~=index3 || index2~=index4 || LabelMatrix(3,end)~=LabelMatrix(3,1) || length(LabelMatrix)<4
%不能这样写,因为每次操作index1,index2都有变化
while LabelMatrix(1,end)~=index3 || LabelMatrix(2,end)~=index4 || LabelMatrix(3,end)~=LabelMatrix(3,1) || length(LabelMatrix)<4
%第一个方向
if LabelMatrix(3,end)==3
if Q(max(index1-1,1),index2)==flag && index1-1~=0
LabelMatrix(:,end+1)=[index1;index2;0];
index1=index1-1;
elseif Q(max(index1-1,1),min(index2+1,256))==flag && index1-1~=0 && index2+1~=257
LabelMatrix(:,end+1)=[index1;index2;1];
index1=index1-1;
index2=index2+1;
elseif Q(index1,min(index2+1,256))==flag && index2+1~=257
LabelMatrix(:,end+1)=[index1; index2; 2];
index2=index2+1;
elseif Q(min(index1+1,256),min(index2+1,256))==flag && index2+1~=257 && index1+1~=257
LabelMatrix(:,end+1)=[index1; index2; 3];
index1=index1+1;
index2=index2+1;
elseif Q(min(index1+1,256),index2)==flag && index1+1~=257
LabelMatrix(:,end+1)=[index1; index2; 4];
index1=index1+1;
elseif Q(min(index1+1,256),max(index2-1,1))==flag && index1+1~=257 && index2-1~=0
LabelMatrix(:,end+1)=[index1; index2; 5];
index1=index1+1;
index2=index2-1;
elseif Q(index1,max(index2-1,1))==flag && index2-1~=0
LabelMatrix(:,end+1)=[index1; index2; 6];
index2=index2-1;
elseif Q(max(index1-1,1),max(index2-1,1))==flag && index2-1~=0 && index1-1~=0
LabelMatrix(:,end+1)=[index1; index2; 7];
index2=index2-1;
index1=index1-1;
end
%第二个方向
elseif LabelMatrix(3,end)==4
if Q(max(index1-1,1),min(index2+1,256))==flag && index1-1~=0 && index2+1~=257
LabelMatrix(:,end+1)=[index1;index2;1];
index1=index1-1;
index2=index2+1;
elseif Q(index1,min(index2+1,256))==flag && index2+1~=257
LabelMatrix(:,end+1)=[index1; index2; 2];
index2=index2+1;
elseif Q(min(index1+1,256),min(index2+1,256))==flag && index2+1~=257 && index1+1~=257
LabelMatrix(:,end+1)=[index1; index2; 3];
index1=index1+1;
index2=index2+1;
elseif Q(min(index1+1,256),index2)==flag && index1+1~=257
LabelMatrix(:,end+1)=[index1; index2; 4];
index1=index1+1;
elseif Q(min(index1+1,256),max(index2-1,1))==flag && index1+1~=257 && index2-1~=0
LabelMatrix(:,end+1)=[index1; index2; 5];
index1=index1+1;
index2=index2-1;
elseif Q(index1,max(index2-1,1))==flag && index2-1~=0
LabelMatrix(:,end+1)=[index1; index2; 6];
index2=index2-1;
elseif Q(max(index1-1,1),max(index2-1,1))==flag && index2-1~=0 && index1-1~=0
LabelMatrix(:,end+1)=[index1; index2; 7];
index2=index2-1;
index1=index1-1;
elseif Q(max(index1-1,1),index2)==flag && index1-1~=0
LabelMatrix(:,end+1)=[index1;index2;0];
index1=index1-1;
end
%第三个方向
elseif LabelMatrix(3,end)==5
if Q(index1,min(index2+1,256))==flag && index2+1~=257
LabelMatrix(:,end+1)=[index1; index2; 2];
index2=index2+1;
elseif Q(min(index1+1,256),min(index2+1,256))==flag && index2+1~=257 && index1+1~=257
LabelMatrix(:,end+1)=[index1; index2; 3];
index1=index1+1;
index2=index2+1;
elseif Q(min(index1+1,256),index2)==flag && index1+1~=257
LabelMatrix(:,end+1)=[index1; index2; 4];
index1=index1+1;
elseif Q(min(index1+1,256),max(index2-1,1))==flag && index1+1~=257 && index2-1~=0
LabelMatrix(:,end+1)=[index1; index2; 5];
index1=index1+1;
index2=index2-1;
elseif Q(index1,max(index2-1,1))==flag && index2-1~=0
LabelMatrix(:,end+1)=[index1; index2; 6];
index2=index2-1;
elseif Q(max(index1-1,1),max(index2-1,1))==flag && index2-1~=0 && index1-1~=0
LabelMatrix(:,end+1)=[index1; index2; 7];
index2=index2-1;
index1=index1-1;
elseif Q(max(index1-1,1),index2)==flag && index1-1~=0
LabelMatrix(:,end+1)=[index1;index2;0];
index1=index1-1;
elseif Q(max(index1-1,1),min(index2+1,256))==flag && index1-1~=0 && index2+1~=257
LabelMatrix(:,end+1)=[index1;index2;1];
index1=index1-1;
index2=index2+1;
end
%第四个方向
elseif LabelMatrix(3,end)==6
if Q(min(index1+1,256),min(index2+1,256))==flag && index2+1~=257 && index1+1~=257
LabelMatrix(:,end+1)=[index1; index2; 3];
index1=index1+1;
index2=index2+1;
elseif Q(min(index1+1,256),index2)==flag && index1+1~=257
LabelMatrix(:,end+1)=[index1; index2; 4];
index1=index1+1;
elseif Q(min(index1+1,256),max(index2-1,1))==flag && index1+1~=257 && index2-1~=0
LabelMatrix(:,end+1)=[index1; index2; 5];
index1=index1+1;
index2=index2-1;
elseif Q(index1,max(index2-1,1))==flag && index2-1~=0
LabelMatrix(:,end+1)=[index1; index2; 6];
index2=index2-1;
elseif Q(max(index1-1,1),max(index2-1,1))==flag && index2-1~=0 && index1-1~=0
LabelMatrix(:,end+1)=[index1; index2; 7];
index2=index2-1;
index1=index1-1;
elseif Q(max(index1-1,1),index2)==flag && index1-1~=0
LabelMatrix(:,end+1)=[index1;index2;0];
index1=index1-1;
elseif Q(max(index1-1,1),min(index2+1,256))==flag && index1-1~=0 && index2+1~=257
LabelMatrix(:,end+1)=[index1;index2;1];
index1=index1-1;
index2=index2+1;
elseif Q(index1,min(index2+1,256))==flag && index2+1~=257
LabelMatrix(:,end+1)=[index1; index2; 2];
index2=index2+1;
end
%第五个方向
elseif LabelMatrix(3,end)==7
if Q(min(index1+1,256),index2)==flag && index1+1~=257
LabelMatrix(:,end+1)=[index1; index2; 4];
index1=index1+1;
elseif Q(min(index1+1,256),max(index2-1,1))==flag && index1+1~=257 && index2-1~=0
LabelMatrix(:,end+1)=[index1; index2; 5];
index1=index1+1;
index2=index2-1;
elseif Q(index1,max(index2-1,1))==flag && index2-1~=0
LabelMatrix(:,end+1)=[index1; index2; 6];
index2=index2-1;
elseif Q(max(index1-1,1),max(index2-1,1))==flag && index2-1~=0 && index1-1~=0
LabelMatrix(:,end+1)=[index1; index2; 7];
index2=index2-1;
index1=index1-1;
elseif Q(max(index1-1,1),index2)==flag && index1-1~=0
LabelMatrix(:,end+1)=[index1;index2;0];
index1=index1-1;
elseif Q(max(index1-1,1),min(index2+1,256))==flag && index1-1~=0 && index2+1~=257
LabelMatrix(:,end+1)=[index1;index2;1];
index1=index1-1;
index2=index2+1;
elseif Q(index1,min(index2+1,256))==flag && index2+1~=257
LabelMatrix(:,end+1)=[index1; index2; 2];
index2=index2+1;
elseif Q(min(index1+1,256),min(index2+1,256))==flag && index2+1~=257 && index1+1~=257
LabelMatrix(:,end+1)=[index1; index2; 3];
index1=index1+1;
index2=index2+1;
end
%第六个方向
elseif LabelMatrix(3,end)==0
if Q(min(index1+1,256),max(index2-1,1))==flag && index1+1~=257 && index2-1~=0
LabelMatrix(:,end+1)=[index1; index2; 5];
index1=index1+1;
index2=index2-1;
elseif Q(index1,max(index2-1,1))==flag && index2-1~=0
LabelMatrix(:,end+1)=[index1; index2; 6];
index2=index2-1;
elseif Q(max(index1-1,1),max(index2-1,1))==flag && index2-1~=0 && index1-1~=0
LabelMatrix(:,end+1)=[index1; index2; 7];
index2=index2-1;
index1=index1-1;
elseif Q(max(index1-1,1),index2)==flag && index1-1~=0
LabelMatrix(:,end+1)=[index1;index2;0];
index1=index1-1;
elseif Q(max(index1-1,1),min(index2+1,256))==flag && index1-1~=0 && index2+1~=257
LabelMatrix(:,end+1)=[index1;index2;1];
index1=index1-1;
index2=index2+1;
elseif Q(index1,min(index2+1,256))==flag && index2+1~=257
LabelMatrix(:,end+1)=[index1; index2; 2];
index2=index2+1;
elseif Q(min(index1+1,256),min(index2+1,256))==flag && index2+1~=257 && index1+1~=257
LabelMatrix(:,end+1)=[index1; index2; 3];
index1=index1+1;
index2=index2+1;
elseif Q(min(index1+1,256),index2)==flag && index1+1~=257
LabelMatrix(:,end+1)=[index1; index2; 4];
index1=index1+1;
end
%第七个方向
elseif LabelMatrix(3,end)==1
if Q(index1,max(index2-1,1))==flag && index2-1~=0
LabelMatrix(:,end+1)=[index1; index2; 6];
index2=index2-1;
elseif Q(max(index1-1,1),max(index2-1,1))==flag && index2-1~=0 && index1-1~=0
LabelMatrix(:,end+1)=[index1; index2; 7];
index2=index2-1;
index1=index1-1;
elseif Q(max(index1-1,1),index2)==flag && index1-1~=0
LabelMatrix(:,end+1)=[index1;index2;0];
index1=index1-1;
elseif Q(max(index1-1,1),min(index2+1,256))==flag && index1-1~=0 && index2+1~=257
LabelMatrix(:,end+1)=[index1;index2;1];
index1=index1-1;
index2=index2+1;
elseif Q(index1,min(index2+1,256))==flag && index2+1~=257
LabelMatrix(:,end+1)=[index1; index2; 2];
index2=index2+1;
elseif Q(min(index1+1,256),min(index2+1,256))==flag && index2+1~=257 && index1+1~=257
LabelMatrix(:,end+1)=[index1; index2; 3];
index1=index1+1;
index2=index2+1;
elseif Q(min(index1+1,256),index2)==flag && index1+1~=257
LabelMatrix(:,end+1)=[index1; index2; 4];
index1=index1+1;
elseif Q(min(index1+1,256),max(index2-1,1))==flag && index1+1~=257 && index2-1~=0
LabelMatrix(:,end+1)=[index1; index2; 5];
index1=index1+1;
index2=index2-1;
end
%第八个方向
elseif LabelMatrix(3,end)==2
if Q(max(index1-1,1),max(index2-1,1))==flag && index2-1~=0 && index1-1~=0
LabelMatrix(:,end+1)=[index1; index2; 7];
index2=index2-1;
index1=index1-1;
elseif Q(max(index1-1,1),index2)==flag && index1-1~=0
LabelMatrix(:,end+1)=[index1;index2;0];
index1=index1-1;
elseif Q(max(index1-1,1),min(index2+1,256))==flag && index1-1~=0 && index2+1~=257
LabelMatrix(:,end+1)=[index1;index2;1];
index1=index1-1;
index2=index2+1;
elseif Q(index1,min(index2+1,256))==flag && index2+1~=257
LabelMatrix(:,end+1)=[index1; index2; 2];
index2=index2+1;
elseif Q(min(index1+1,256),min(index2+1,256))==flag && index2+1~=257 && index1+1~=257
LabelMatrix(:,end+1)=[index1; index2; 3];
index1=index1+1;
index2=index2+1;
elseif Q(min(index1+1,256),index2)==flag && index1+1~=257
LabelMatrix(:,end+1)=[index1; index2; 4];
index1=index1+1;
elseif Q(min(index1+1,256),max(index2-1,1))==flag && index1+1~=257 && index2-1~=0
LabelMatrix(:,end+1)=[index1; index2; 5];
index1=index1+1;
index2=index2-1;
elseif Q(index1,max(index2-1,1))==flag && index2-1~=0
LabelMatrix(:,end+1)=[index1; index2; 6];
index2=index2-1;
end
end
end
% %轨迹信息录入Trace数组
Trace=[Trace LabelMatrix];
num1=num1+1;
% Trace(end+1)=LabelMatrix(1,1);
% Trace(end+1)=LabelMatrix(2,1);
% %此时end在随时变化
% for index5=1:length(LabelMatrix)
% Trace(end+1)=LabelMatrix(3,index5);
% end
% Trace(end+1)=10000;
LabelMatrix=[];
%该连通域内所值置0;(此处大可优化,简单起见暂用此法)
for index7=1:m
for index8=1:n
if Q(index7,index8)==flag
Q(index7,index8)=0;
end
end
end
end
end
end
H=zeros(m,n);
for index10=1:length(Trace)
H(Trace(1,index10),Trace(2,index10))=1;
end
Trace
H=mat2gray(H);
imtool(H);
end