%基本思想：利用pcnn，每當(dāng)有一批像素對應(yīng)的神經(jīng)元點火，對像素值進(jìn)行一次修正。第n次點火的所有神經(jīng)元用矩陣B（n）表示，
%已經(jīng)點火的像素位置標(biāo)記為‘1’，未點火的標(biāo)記為‘0’。通過一個3*3的模板滑過B（n），判斷若模板內(nèi)的值全為‘1’或全為‘0’，
%則這些像素值不進(jìn)行處理，否則若模板中心的值為‘1’，則增加該位置的像素值的大小，中心值為‘0’，則減小像素值。該功能由xiugai(B,K)函數(shù)實現(xiàn)
%Beta取負(fù)值來抑制周圍的神經(jīng)元點火，因為輸入pcnn（X）的是模糊圖像，抑制之后使處理的像素值更加接近原圖像的像素值（抑制參數(shù)通過試驗確定，并不精確）
%總之運行結(jié)果只是對像素值變化較大的邊緣區(qū)域進(jìn)行了像素值的修正。
%通過不斷調(diào)整參數(shù)使峰值信噪比有所升高，但是升高得不是太多
function filter(I) %主函數(shù)
I = imread('e:\程序\000.jpg');
[a,b] = size(I(:,:,1));
D = zeros(a-2,b-2);
E = zeros(a-2,b-2);
F = zeros(a-2,b-2);
K = filter2(fspecial('average',3),I(:,:,1)); %使用3*3且每個單元都為1的模糊核對原圖像進(jìn)行模糊處理
[J,P] = deconvblind(K,fspecial('average',3),10);
subplot(221), imshow(uint8(I(:,:,1))); title('原圖像');
subplot(222), imshow(uint8(K)); title('模糊圖像');
[Edge,Numberofarea,Yz] = pcnn(K);         %對模糊圖像用pcnn處理
%subplot(223), imshow(uint8(Yz - double(I(:,:,1)) + 128)); title('差值');
subplot(223), imshow(uint8(J)); title('R_L去卷積');
for i0 = 3 : a-2
    for i1 = 3 :b-2
        D(i0-2,i1-2) = K(i0,i1);
        E(i0-2,i1-2) = Yz(i0,i1);   %Yz為pcnn處理后的圖像
        F(i0-2,i1-2) = I(i0,i1,1);
    end
end
PSNR1 = 10*log10(a*b*max(max(double(F)))*max(max(double(F)))/sum(sum((D - double(F)).*(D - double(F)))));
PSNR2 = 10*log10(a*b*max(max(double(F)))*max(max(double(F)))/sum(sum((E - double(F)).*(E - double(F))))); %計算峰值信噪比
disp(PSNR1);
disp(PSNR2);
disp(sum(sum(K - double(I(:,:,1)))));
disp(sum(sum(Yz - double(I(:,:,1)))));
%Yzz = histeq(uint8(Yz),256);                     %直方圖均衡化
subplot(224), imshow(uint8(Yz)); title('修正后結(jié)果');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Edge,Numberofarea,Xz] = pcnn(X)           %pcnn處理主函數(shù)，后面幾個函數(shù)都被該函數(shù)調(diào)用
Xz = double(X);
X1 = jiabian(X);
Beta = -0.23;                     %Beta取負(fù)值來抑制周圍的神經(jīng)元點火,可調(diào)整
Yuzhi = 254;
Decay = 0.4;
[a,b] = size(X);
Threshold = zeros(a,b);
S = zeros(a+2,b+2);
B = zeros(a,b);
Y = zeros(a,b);
T = zeros(a,b);
Yz = zeros(a,b);
Edge = zeros(a,b); Numberofarea = zeros(a,b); Numberofarea_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)];
            W = [X1(i0-1,i1-1) X1(i0-1,i1) X1(i0-1,i1+1);
                 X1(i0,i1-1) X1(i0,i1) X1(i0,i1+1);
                 X1(i0+1,i1-1) X1(i0+1,i1) X1(i0+1,i1+1)];
             
             L = 0.125*sum(sum((2*V-1).*(X1(i0,i1)-W)))/X1(i0,i1);
             
             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) < 30
                 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;
                 Y(i0-1,i1-1) = 0;
             end
        end
    end
%Threshold(find(B ~= 1)) = exp(-Decay)*Threshold(find(B ~= 1)); %閾值按指數(shù)衰減
Threshold(find(B ~= 1)) = Threshold(find(B ~= 1))-2; %閾值設(shè)置成線性衰減
if n ~= 1
    Edge = Edge+judge_edge(Y,n);
    Y(find(Edge<0)) = 0;
    Numberofarea = Numberofarea+Y*10*(30-n);
    Num = Num_1;
end
if n == 1
    S = zeros(a+2,b+2);
else
    S = bianhuan(T);
end
n = n+1;
Numberofarea_1 = zeros(a,b);
Yz = xiugai(B,X);
Xz = Xz + Yz;
if Xz(i0-1,i1-1) > 255
    Xz(i0-1,i1-1) = 255;
end
if Xz(i0-1,i1-1) < 0
    Xz(i0-1,i1-1) = 0;
end
end
%disp(n);
%disp(min(min(uint16(X))));
%disp(max(max(uint16(X))));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Y = jiabian(X)
[m,n] = size(X);
Y = zeros(m+2,n+2);
for i = 1:m+2
    for j = 1:n+2
        if i==1 & j~=1 & j~=n+2
            Y(i,j) = X(1,j-1);
        elseif j==1 & i~=1 & i~=m+2
            Y(i,j) = X(i-1,1);
        elseif i~=1 & j==n+2 & i~=m+2
            Y(i,j) = X(i-1,n);
        elseif i==m+2 & j~=1 & j~=n+2
            Y(i,j) = X(m,j-1);
        elseif i==1 & j==1
            Y(i,j) = X(i,j);
        elseif i==1 & j==n+2
            Y(i,j) = X(1,n);
        elseif i==(m+2) & j==1
            Y(i,j) = X(m,1);
        elseif i==m+2 & j==n+2
            Y(i,j) = X(m,n);
        else
            Y(i,j) = X(i-1,j-1);
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Y = bianhuan(X)
[m,n] = size(X);
Y = zeros(m+2,n+2);
for i = 1:m+2
    for j = 1:n+2
        if i==1 | j==1 | i==m+2 | j==n+2
            Y(i,j) = 0;
        else
            Y(i,j) = X(i-1,j-1);
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Y = judge_edge(X,n)
[a,b] = size(X);
T = jiabian(X);
Y = zeros(a,b);
W = zeros(a,b);
for i = 2:a+1
    for j = 2:b+1
        if T(i,j) == 1 & ((T(i-1,j)==0 & T(i+1,j) == 0) | (T(i,j-1)==0 & T(i,j+1) == 0) | (T(i-1,j-1)==0 & T(i+1,j+1)==0) | (T(i+1,j+1) == 0 & T(i-1,j+1)==0))
            Y(i-1,j-1) = -n;
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Y = xiugai(B,K)
[a,b] = size(B);
C = bianhuan(B);
L = jiabian(K);
for i0 = 2 : a+1
        for i1 = 2 : b+1
             M = [C(i0-1,i1-1) C(i0-1,i1) C(i0-1,i1+1);
                 C(i0,i1-1) C(i0,i1) C(i0,i1+1);
                 C(i0+1,i1-1) C(i0+1,i1) C(i0+1,i1+1)];
             N = [L(i0-1,i1-1) L(i0-1,i1) L(i0-1,i1+1);
                 L(i0,i1-1) L(i0,i1) L(i0,i1+1);
                 L(i0+1,i1-1) L(i0+1,i1) L(i0+1,i1+1)];
             if sum(sum(M)) > 0 & sum(sum(M)) < 9
                 if C(i0,i1) == 1
                     Y(i0-1,i1-1) = 0.05*abs(sum(sum(0.125*(1-M).*(K(i0-1,i1-1)-N)))); %模板中心為‘1’，增加像素值
                 elseif C(i0,i1) == 0               
                     Y(i0-1,i1-1) = -0.05*abs(sum(sum(0.125*M.*(N-K(i0-1,i1-1)))));  %模板中心為‘0’，減少像素值
             
                 end      
             else
                     Y(i0-1,i1-1) = 0;              
             end
        end
end