%function [Edge, Numberofaera]=pcnn(X)
clear all
I=imread('lena256.bmp');
figure,imshow(I);
X=double(I);
weight=[0.25  1 0.25;1 0 1;0.25 1 0.25];
%Beta=0.4;
Beta=0.80;
Yuzhi=195;
%Decay=0.3;
Decay=0.70;
[a,b]=size(X);
Threshold=zeros(a,b);
S=zeros(a+2,b+2);
B=zeros(a,b);
C=zeros(a,b);
Y=zeros(a,b);
Edge=zeros(a,b);
Numberofaera=zeros(a,b);
Numberofaera_1=zeros(a,b);
Num_1=0;
Num=0;
n=1;
%for tt=1:1,
while(sum(sum(B))~=256*256)
    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;
                    C(i0-1,i1-1)=n;
                end
            else
                T(i0-1,i1-1)=0;
                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;
        [Numberofara_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_=zeros(a,b);
end

%end%%%%




S=S*(-1)+1;
t=0;
for i=1:256,
    for j=1:256,
        t=t+S(i,j);
    end
end
        
t=t/256/256;
%S=S*(-1)+1
figure(2);
imshow(Y,[]);  
%figure
%Y=Y*(-1)+1
%imshow(Y)