map1=imread('p.jpg');
[row,col,dep]=size(map1);   %行，列，深度值
map=zeros(row,col);
pixsum=row*col;                %像素总数
%灰度化
for i=1:row
    for j=1:col
        map(i,j)=0.11*map1(i,j,1)+0.59*map1(i,j,2)+0.3*map1(i,j,3);
    end
end
map2=map;                         %灰度图副本
%灰度化完成,输出灰度图
figure(1)
imshow(map,[])
title('原始灰度图像')
B=zeros(row,col);       %存储输出图像
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%根据通道二值化
    for i=1:row
        for j=1:col
           if map1(i,j,2)>map1(i,j,1)&&map1(i,j,2)>map1(i,j,3)
              B(i,j)=1;
           else
              B(i,j)=0;
            end
        end
    end
%figure(2)
%imshow(B,[])
%title('二值化后图像')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%以下为边界提取程序
A=ones(3,3);%腐蚀剂
Y1=zeros(row,col);%腐蚀后图像
for i=2:row-1
    for j=2:col-1
        a=sum(sum(B(i-1:i+1,j-1:j+1).*A));
        if a==9
            Y1(i,j)=1;
        else
            Y1(i,j)=0;
        end
    end
end
Y=B-Y1;
num=sum(sum(Y));
fprintf('轮廓图像中边界长度为%d\n',num)
figure(3)
imshow(Y,[])
title('轮廓图像')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%以下为直线检测程序
N=101;       %theta范围为0-pi/2,将其100等分，在进行正余弦计算时再转换成弧度。
rohmax=round(sqrt(row^2+col^2))+1;   %roh的最大值应为
RT=zeros(N,rohmax);              %roh-theta计数器。
for i=1:row
    for j=1:col
        if Y(i,j)==1              %若为边界上的点，则进行计算。
        for theta=0:pi/200:pi/2
            roh=i*cos(theta)+j*sin(theta);
            m=round(theta*200/pi+1);
            n=round(roh);
            RT(m,n)=RT(m,n)+1;
        end
        end
    end
end
%统计完成
p=90;%包含p个以上像素的直线标记出来
cnt=0;%记录直线条数
q=[];
for i=1:N
    for j=1:rohmax
        if RT(i,j)>=p        
            q=[q,RT(i,j)];
            cnt=cnt+1;
            RT(i,j)=1;
        else 
            RT(i,j)=0;
        end
    end
end
%标记完成
Y2=zeros(row,col);                      %标记图
for m=1:N
    for n=1:rohmax
        if RT(m,n)==1                   %遍历RH图，只要roh和theta满足条件，计算响应的（x,y）
            for i=1:row
                for j=1:col
                    if Y(i,j)==1        %遍历边界点，检测是否在满足条件的直线上。
                        z=i*cos((m-1)*pi/200)+j*sin((m-1)*pi/200);  %z=xcostheta+ysintheta
                        w=round(z);
                        if w==n
                            Y2(i,j)=1;
                        end
                    end
                end
            end
        end
    end
end
%检测完成，标记在Y2中
for i=1:row
    for j=1:col
        if Y2(i,j)==1
            map2(i,j)=255;          %检测到的直线标记在灰度图副本中
        end
    end
end
figure(4)
imshow(map2,[])
title('直线检测结果')
fprintf('\n%d条直线长度分别为：\n',cnt)
lentth=q'