PCNN(4)边缘检测

边缘提取首先提供一段MATLAB代码吧：

I=imread('lena.bmp');%  提取图像
BW1=edge(I,'sobel'); %用SOBEL算子进行边缘检测
BW2=edge(I,'roberts');%用Roberts算子进行边缘检测
BW3=edge(I,'prewitt'); %用prewitt算子进行边缘检测
BW4=edge(I,'log'); %用log算子进行边缘检测
BW5=edge(I,'canny'); %用canny算子进行边缘检测
h=fspecial('gaussian’,5);
BW6=edge(I,’canny’);
subplot(2,3,1), imshow(BW1);
title(‘sobel edge check’);
subplot(2,3,2), imshow(BW2);
title(‘sobel edge check’);
subplot(2,3,3), imshow(BW3);
title(‘prewitt edge check’);
subplot(2,3,4), imshow(BW4);
title(‘log edge check’);
subplot(2,3,5), imshow(BW5);
title(‘canny edge check’);
subplot(2,3,6), imshow(BW6);
title(‘gasussian&canny edge check’);%此为用高斯滤波后Canny算子边缘检测结果

	加法运算P	乘法运算M
Roberts算子	3*N2	0
Sobel算子	11*N2	2*N2
Prewitt算子	11*N2	0
Krisch算子	56*N2	0
Lapalacian算子	4*N2	0

       可以看出，Krisch算子的运算量比较大。其次在边缘检测中边缘定位能力和噪声抑制能力方面，有的算子边缘定位能力强，有的抗噪声能力比较好：Roberts算子利用局部差分算子寻找边缘，边缘定位精度较高，但容易丢失一部分边缘，同时由于没经过图像平滑计算，不能抑制噪声。该算子对具有陡峭的低噪声图像响应最好；Sobel算子和Prewitt算子都是对图像进行差分和滤波运算，差别只是平滑部分的权值有些差异，对噪声具有一定的抑制能力，不能完全排除检测结果中出现伪边缘。这两个算子的边缘定位比较准确和完整，但容易出现边缘多像素宽。对灰度渐变和具有噪声的图像处理的较好；Krisch算子对8个方向边缘信息进行检测，因此有较好的边缘定位能力，并且对噪声有一定的抑制作用，该算子的边缘定位能力和抗噪声能力比较理想；Laplacian算子是二阶微分算子，对图像中的阶跃型边缘点定位准确且具有旋转不变性即无方向性。但该算子容易丢失一部分边缘的方向信息，造成不连续的检测边缘，同时抗噪声能力比较差，比较适用于屋脊型边缘检测。

初始值设置：

%%%%%初始值设置
Weight=[0.5 1 0.5;1 0 1;0.5 1 0.5];% 权值矩阵
Beta=0.4;   %内部活动项的连接因子
Yuzhi=240;
Decay=0.3;   %衰减
[a,b]=size(X);
Threshold=zeros(a,b);  %门限值
S=zeros(a+2,b+2);
B=zeros(a,b);%标记样板，表示该Pixel是否被激活过
Y=zeros(a,b);
Edge=zeros(a,b);
Numberofaera=zeros(a,b);
Numberofaera_1=zeros(a,b);
Num_1=0;
Num=0;
n=1;

点火公式：

%%%点火公式
while(sum(sum(B))~=a*b)
    for i0=2:a+1
         for i1=2:b+1
             V=[S(i0-1,i1-1)  S(i0-1,i1) S(i0-1,i1+1);
                S(i0,i1-1)   S(i0,i1)   S(i0,i1+1);
                S(i0+1,i1-1) S(i0+1,i1) S(i0+1,i1+1)];
              L=sum(sum(V.*Weight));
              F=X(i0-1,i1-1);
              U=double(F)*(1+Beta*double(L));  
             if U>=Threshold(i0-1,i1-1)|Threshold(i0-1,i1-1)<60
                T(i0-1,i1-1)=1;
                Threshold(i0-1,i1-1)=Yuzhi;
                Y(i0-1,i1-1)=1;
                if n==1
                    B(i0-1,i1-1)=0;%避免第一次全部激发造成的影响
                else
                    B(i0-1,i1-1)=1;%已发射过的标记
                    Threshold(i0-1,i1-1)=1000000;%相当于不会第二次激活
                end
             else
                 T(i0-1,i1-1)=0;%no use?
                 Y(i0-1,i1-1)=0;
             end
         end
    end
    Threshold(find(B~=1))=exp(-Decay)*Threshold(find(B~=1));
    %被激活过的象素不再参与迭代过程
    if n~=1
       Edge=Edge+judge_edge(Y,n);
       Y(find(Edge<0))=0;%边界点被置零，Y本来是激发的象素，现在边界被置0，也不能说白激发了，B矩阵有记录！当然下次就休想再被激发了
       [Numberofaera_1,Num_1]=bwlabel(Y,4);
       for i=1:a
           for j=1:b
               if Numberofaera_1(i,j)~=0 
                   Numberofaera_1(i,j)=Numberofaera_1(i,j)+Num;
               end
           end
       end
       Numberofaera=Numberofaera+Numberofaera_1;
       Num=Num_1;
    end
    if n==1
        S=zeros(a+2,b+2);
    else
        S=bianhuan(T);
    end
    n=n+1;
    Numberofaera_1=zeros(a,b); 
    if mod(n,1)==0
        disp(n)
        disp(sum(sum(B)))
    end
    
    figure(1)
    imshow(Numberofaera)
    drawnow
    
    frame=getframe(1);
    im=frame2im(frame);
    [imind,cm]=rgb2ind(im,256);
  
end

具体公式参见PCNN（2）。
参考： http://blog.sina.com.cn/s/blog_6163bdeb0100nihw.html