pcnn.m

在MATLAB上编写的PCNN去除噪声的程序

%基本思想：利用pcnn，每当有一批像素对应的神经元点火，对像素值进行一次修正。第n次点火的所有神经元用矩阵B（n）表示，
%已经点火的像素位置标记为‘1’，未点火的标记为‘0’。通过一个3*3的模板滑过B（n），判断若模板内的值全为‘1’或全为‘0’，
%则这些像素值不进行处理，否则若模板中心的值为‘1’，则增加该位置的像素值的大小，中心值为‘0’，则减小像素值。该功能由xiugai(B,K)函数实现
%Beta取负值来抑制周围的神经元点火，因为输入pcnn（X）的是模糊图像，抑制之后使处理的像素值更加接近原图像的像素值（抑制参数通过试验确定，并不精确）
%总之运行结果只是对像素值变化较大的边缘区域进行了像素值的修正。
%通过不断调整参数使峰值信噪比有所升高，但是升高得不是太多
function filter(I) %主函数
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,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('修正后结果');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Edge,Numberofarea,Xz] = pcnn(X)           %pcnn处理主函数，后面几个函数都被该函数调用
Xz = double(X);
X1 = jiabian(X);
Beta = -0.23;                     %Beta取负值来抑制周围的神经元点火,可调整
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)); %阈值按指数衰减
Threshold(find(B ~= 1)) = Threshold(find(B ~= 1))-2; %阈值设置成线性衰减
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