【MATLAB图像融合】[15]一个自适应PCNN的DEMO详解

出处: Qu Xiao-Bo <qxb_xmu [at] yahoo.com.cn> Aug.28,2008
先上代码，依然来自屈小波教授：

A=imread('1.jpg');
B=imread('2.jpg');
im=imread('1.jpg');
im1=imread('2.jpg');
im=double(im);
im1=double(im1);
[p,q]=size(im);
link_arrange=5;
np=20;

F_NA=Normalized(im);
F_NB=Normalized(im1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha_L=1;
alpha_Theta=2.0;
alpha_F=0.1;
beta=0; 
 vF=0.5;
vL=0.2;
vTheta=20;
L1=zeros(p,q);
L2=zeros(p,q);
U1=zeros(p,q);
U2=zeros(p,q);
Y1=zeros(p,q);
Y2=zeros(p,q);
F1=Y1;
F2=Y1;
Y0=zeros(p,q);
Y01=zeros(p,q);
R=zeros(p,q);
Theta1=zeros(p,q);
Theta2=zeros(p,q);

% Compute the linking strength.
center_x=round(link_arrange/2);
center_y=round(link_arrange/2);
W=zeros(link_arrange,link_arrange);
for i=1:link_arrange
    for j=1:link_arrange
        if (i==center_x)&&(j==center_y)
            W(i,j)=0;
        else
            W(i,j)=1./sqrt((i-center_x).^2+(j-center_y).^2);
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I1=F_NA;
I2=F_NB;
for n=1:np
    [a,b]=size(W);
      for i=1:p
       for j=1:q
    for c=1:a
        for d=1:b
           beta=beta+[im(c,d)-im(c+1,d)].^2+[im(c,d)-im(c,d+1)].^2;
        end
    end
    K1(i,j)=conv2(Y1(i,j),W,'same');
    K2(i,j)=conv2(Y2(i,j),W,'same');
    F1(i,j)=exp(-alpha_F)*F1(i,j)+vF*K1(i,j)+I1(i,j);
    F2(i,j)=exp(-alpha_F)*F2(i,j)+vF*K2(i,j)+I2(i,j);
    L1(i,j)=exp(-alpha_L)*L1(i,j)+vL*K1(i,j);
    L2(i,j)=exp(-alpha_L)*L2(i,j)+vL*K2(i,j);
    Theta1(i,j)=exp(-alpha_Theta)*Theta1(i,j)+vTheta*Y1(i,j);
    Theta2(i,j)=exp(-alpha_Theta)*Theta2(i,j)+vTheta*Y2(i,j);
    U1(i,j)=F1(i,j).*(1+beta*L1(i,j));
    U2(i,j)=F2(i,j).*(1+beta*L2(i,j));
    Y1(i,j)=im2double(U1(i,j)>Theta1(i,j));
    Y2(i,j)=im2double(U2(i,j)>Theta2(i,j));
    Y0(i,j)=Y0(i,j)+Y1(i,j);
    Y01(i,j)=Y01(i,j)+Y2(i,j);
    end
      end
      Y0=medfilt2(Y0,[3,3]);
      Y01=medfilt2(Y01,[3,3]);
      for i=1:p
          for j=1:q
              if abs(Y0(i,j)-Y01(i,j))<=0.015
                 R(i,j)=(I1(i,j)+I2(i,j))/2;
              end 
              if abs(Y0(i,j)-Y01(i,j))>0.015 &&Y0(i,j)>Y01(i,j)
                 R(i,j)=I1(i,j);
              end
              if abs(Y0(i,j)-Y01(i,j))>0.015 && Y0(i,j)<Y01(i,j)
                 R(i,j)=I2(i,j);
              end 
          end
      end  
end
subplot 131;imshow(A);
subplot 132;imshow(B);
subplot 133;imshow(R);

function [normalized_matrix,cmin,cmax]=Normalized(matrix)
input_matrix=abs(matrix);
Max_input=max(input_matrix(:));
Min_input=min(input_matrix(:));
min_matrix=ones(size(input_matrix)).*Min_input;
normalized_matrix=(input_matrix-min_matrix)./(Max_input-Min_input+eps);
cmin=Min_input;
cmax=Max_input;
end;

对本代码进行以下解释：
①、本代码基于PCNN，而不是S-PCNN。
②、PCNN的基本原理在前章有解释。
③、算法的基本流程是归一化；定义PCNN的如下参数：W,LINK RANGE,TIME,aF,vF,aL,vL,aT,vT;定义EOL（energy of laplacian）作为β的自适应值；通过阈值对点火图进行融合得到融合图像。
④、归一化函数算法：选出一个最大值和最小值；让矩阵中所有的数值都减去最小值；得到的新矩阵再除以（最大值-最小值）就得到归一化结果；过程中取绝对值计算。
⑤、权值（W）和链接域（LINK RANGE）：通常取奇数，经测试，在W=3时算法却莫名其妙无法计算结果，检查无果，W=5,6,7等时都可正常运行。
⑥、需要注意，步骤③最终用点火图完成融合，不是把点火图融合在一起，而是根据点火图的情况去相应地对源图像像素进行融合操作，故PCNN图像融合也属于像素级的图像融合。