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
%灰度化完成,输出灰度图
figure(1)
imshow(map,[])
title('原始灰度图像')
Y=zeros(row,col);       %存储输出图像
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5根据通道二值化
    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)
              Y(i,j)=1;
           else
              Y(i,j)=0;
            end
        end
    end
figure(2)
imshow(Y,[])
title('二值化后图像')
%以下为骨架提取程序
temp=ones(row,col);
N=0;           %统计邻域内非零像素个数。
T=0;           %统计0-1变换次数
for k=1:80      %大循环次数
    for i=2:row-1
        for j=2:col-1
            a=[Y(i,j),Y(i,j-1),Y(i+1,j-1),Y(i+1,j),Y(i+1,j+1),Y(i,j+1),Y(i-1,j+1),Y(i-1,j),Y(i-1,j-1)];%邻域元素排序。中间为1，左边为2，逆时针分别为2---9
            if a(1)==1
                for m=1:8
                    if a(m)==1
                        N=N+1;
                    end
                end
                for n=1:7
                    if a(n)==0&&a(n+1)==1
                        T=T+1;
                    end
                end
                if a(8)==0&&a(1)==1
                    T=T+1;
                end
            end
            if (N>=2)&&(N<=6)&&(T==1)&&(a(2)*a(4)*a(6)==0)&&(a(4)*a(6)*a(8)==0)
                temp(i,j)=0;            %标记要删除的点，标志为1
            end
            N=0;
            T=0;
        end
    end
    %删除点
    for i=1:row
        for j=1:col
            if temp(i,j)==0
                Y(i,j)=0;
            end
        end
    end
    %第一步完成
    for i=2:row-1
        for j=2:col-1
            a=[Y(i,j),Y(i,j-1),Y(i+1,j-1),Y(i+1,j),Y(i+1,j+1),Y(i,j+1),Y(i-1,j+1),Y(i-1,j),Y(i-1,j-1)];%邻域元素排序。中间为1，左边为2，逆时针分别为2---9
            if a(1)==1
                for m=2:9
                    if a(m)==1
                        N=N+1;
                    end
                end
                for n=1:7
                    if a(n)==0&&a(n+1)==1
                        T=T+1;
                    end
                end
                if a(8)==0
                    T=T+1;
                end
            end
            if (N>=2)&&(N<=6)&&(T==1)&&(a(2)*a(4)*a(8)==0)&&(a(2)*a(6)*a(8)==0)
                temp(i,j)==0;            %标记要删除的点，标志为1
            end
            N=0;
            T=0;
        end
    end
    %删除点
    for i=1:row
        for j=1:col
            if temp(i,j)==0
                Y(i,j)=0;
            end
        end
    end    
    %第二步完成
end
%完成，输出
figure(3)
imshow(Y,[])
title('提取的骨架')
%以下为交叉点检测及标志，交叉点标注在灰度图中
temp1=zeros(row,col);
for i=2:row-1
    for j=2:col-1
        c=sum(sum(Y(i-1:i+1,j-1:j+1)));
        if c>=4
            temp1(i,j)=1;
        end
    end
end
for i=1:row
    for j=1:col
        if temp1(i,j)==1
           map(i,j)=255;
        end
    end
end
%标注完成，显示灰度图
figure(4)
imshow(map,[])
title('标注的交叉点')