金字塔变换的图像融合Matlab源码

附录1 金字塔变换图像融合方法程序

% 拉普拉斯金字塔融合函数

function Y = fuse_lap(M1, M2, zt, ap, mp)

% M1、M2为源图像

% zt为融合层数，ap为高频子带图像选择系数，mp为低频子带图像选择系数[z1 s1] = size(M1);

[z2 s2] = size(M2);

if (z1 ~= z2) || (s1 ~= s2)

error('输入源图像大小不一致');

end;

% 高斯窗口函数

w = [1 4 6 4 1] / 16;

E = cell(1,zt);

for i1 = 1:zt

[z s] = size(M1);

zl(i1) = z; sl(i1) = s;

% 图像尺寸为奇数还是偶数

if (floor(z/2) ~= z/2), ew(1) = 1; else ew(1) = 0; end;

if (floor(s/2) ~= s/2), ew(2) = 1; else ew(2) = 0; end;

% 若为奇数，扩展为偶数

if (any(ew))

M1 = adb(M1,ew);

M2 = adb(M2,ew);

end;

% M1与M2低通滤波

G1 = conv2(conv2(es2(M1,2), w, 'valid'),w', 'valid');

G2 = conv2(conv2(es2(M2,2), w, 'valid'),w', 'valid');

% G1与G2下采样、上采样低通滤波后的膨胀序列

M1T = conv2(conv2(es2(undec2(dec2(G1)), 2), 2*w, 'valid'),2*w', 'valid');

M2T = conv2(conv2(es2(undec2(dec2(G2)), 2), 2*w, 'valid'),2*w', 'valid');

% 高频子带图像系数选择

E(i1) = {selg(M1-M1T, M2-M2T, ap)};

% G11与G2下采样

M1 = dec2(G1);

M2 = dec2(G2);

end;

% 低频子带图像系数选择

M1 = selh(M1,M2,mp);

% 图像重构

for i1 = zt:-1:1

M1T = conv2(conv2(es2(undec2(M1), 2), 2*w, 'valid'), 2*w', 'valid');

M1 = M1T + E{i1};

% 选择图像有效区域

M1 = M1(1:zl(i1),1:sl(i1));

end;

Y = M1;

end;

% 对比度金字塔融合函数

function Y = fuse_con(M1, M2, zt, ap, mp)

[z1 s1] = size(M1);

[z2 s2] = size(M2);

if (z1 ~= z2) | (s1 ~= s2)

error('输入图像尺寸大小不一致');

end;

w = [1 4 6 4 1] / 16;

eps = 1e-6;

E = cell(1,zt);

for i1 = 1:zt

[z s] = size(M1);

zl(i1) = z; sl(i1) = s;

if (floor(z/2) ~= z/2), ew(1) = 1; else, ew(1) = 0; end;

if (floor(s/2) ~= s/2), ew(2) = 1; else, ew(2) = 0; end;

if (any(ew))

M1 = adb(M1,ew);

M2 = adb(M2,ew);

end;

G1 = conv2(conv2(es2(M1,2), w, 'valid'),w', 'valid');

G2 = conv2(conv2(es2(M2,2), w, 'valid'),w', 'valid');

M1T = conv2(conv2(es2(undec2(dec2(G1)), 2), 2*w, 'valid'),2*w', 'valid'); M2T = conv2(conv2(es2(undec2(dec2(G2)), 2), 2*w, 'valid'),2*w', 'valid'); E(i1) = {selg(M1./(M1T+eps)-1, M2./(M2T+eps)-1, ap)};

M1 = dec2(G1);

M2 = dec2(G2);

end;

M1 = selh(M1,M2,mp);

for i1 = zt:-1:1

M1T = conv2(conv2(es2(undec2(M1), 2), 2*w, 'valid'), 2*w', 'valid');

M1 = (M1T+eps) .* (E{i1}+1);

M1 = M1(1:zl(i1),1:sl(i1));

end;

Y = M1;

end;

% 梯度金字塔融合函数

function Y = fuse_gra(M1, M2, zt, ap, mp)

[z1 s1] = size(M1);

[z2 s2] = size(M2);

if (z1 ~= z2) | (s1 ~= s2)

error('输入图像大小不一致');

end;

w = [1 4 6 4 1] / 16;

v = [1 2 1] / 4;% 核函数

% 梯度算子

d1 = [1 -1];

d2 = [0 -1; 1 0] / sqrt(2);

d3 = [-1 1];

d4 = [-1 0; 0 1] / sqrt(2);

% 计算导数

d1e = conv2(d1,d1);

d1e = [zeros(1,3); d1e; zeros(1,3)];

d2e = conv2(d2,d2);

d3e = d1e';

d4e = conv2(d4,d4);

E = cell(1,zt);

for i1 = 1:zt

[z s] = size(M1);

zl(i1) = z; sl(i1) = s;

if (floor(z/2) ~= z/2), ew(1) = 1; else, ew(1) = 0; end;

if (floor(s/2) ~= s/2), ew(2) = 1; else, ew(2) = 0; end;

if (any(ew))

M1 = adb(M1,ew);

M2 = adb(M2,ew);

end;

% 梯度金字塔的建立

G1 = conv2(conv2(es2(M1,2), w, 'valid'),w', 'valid');

G2 = conv2(conv2(es2(M2,2), w, 'valid'),w', 'valid');

Z1 = es2(M1+conv2(conv2(es2(M1, 1), v, 'valid'), v', 'valid'), 1); Z2 = es2(M2+conv2(conv2(es2(M2, 1), v, 'valid'), v', 'valid'), 1);

B = zeros(size(M1));

% 方向拉普拉斯金字塔的建立

D1 = conv2(Z1, d1e, 'valid');

D2 = conv2(Z2, d1e, 'valid');

B = B + selg(D1, D2, ap);

D1 = conv2(Z1, d2e, 'valid');

D2 = conv2(Z2, d2e, 'valid');

B = B + selg(D1, D2, ap);

D1 = conv2(Z1, d3e, 'valid');

D2 = conv2(Z2, d3e, 'valid');

B = B + selg(D1, D2, ap);

D1 = conv2(Z1, d4e, 'valid');

D2 = conv2(Z2, d4e, 'valid');

B = B + selg(D1, D2, ap);

E(i1) = {-B/8};

M1 = dec2(G1);

M2 = dec2(G2);

end;

M1 = selh(M1,M2,mp);

for i1 = zt:-1:1

M1T = conv2(conv2(es2(undec2(M1), 2), 2*w, 'valid'), 2*w', 'valid'); M1 = M1T + E{i1};

M1 = M1(1:zl(i1),1:sl(i1));

end;

Y = M1;

end;

% 图像2抽取函数

function Y = dec2(X);

[a b] = size(X);

Y = X(1:2:a, 1:2:b);

end;

% 图像2插值函数

function Y = undec2(X)

[z s] = size(X);

Y = zeros(2*z, 2*s);

Y(1:2:2*z,1:2:2*s) = X;

end;

% 图像扩展函数

function Y = adb(X, bd)

[z s] = size(X);

Y = zeros(z+bd(1),s+bd(2));

Y(1:z,1:s) = X;

if (bd(1) > 0)

Y(z+1:z+bd(1),1:s) = X(z-1:-1:z-bd(1),1:s);

end;

if (bd(2) > 0)

Y(1:z,s+1:s+bd(2)) = X(1:z,s-1:-1:s-bd(2));

end;

if (bd(1) > 0 & bd(2) > 0)

Y(z+1:z+bd(1),s+1:s+bd(2)) = X(z-1:-1:z-bd(1),s-1:-1:s-bd(2));

end;

% 图像边缘像素处理函数

function Y = es2(X, n)

[z s] = size(X);

Y = zeros(z+2*n, s+2*n);

Y(n+1:n+z,n:-1:1)= X(:,2:1:n+1);

Y(n+1:n+z,n+1:1:n+s)= X;

Y(n+1:n+z,n+s+1:1:s+2*n)= X(:,s-1:-1:s-n);

Y(n:-1:1,n+1:s+n)= X(2:1:n+1,:);

Y(n+z+1:1:z+2*n,n+1:s+n)= X(z-1:-1:z-n,:);

% 低频子带图像选择函数

function Y = selh(M1, M2, mp)

switch (mp)

case 1, Y = M1;% 选择图像M1低频子带图像作为融合函数低频成分

case 2, Y = M2;% 选择图像M2低频子带图像作为融合函数低频成分

case 3, Y = (M1 + M2)/2;% 选择图像M1与M2加权平均低频子带图像作为融合函数低频成分

otherwise, error('低频成分选择错误');

end;

% 高频子带图像选择函数

function Y = selg(M1, M2, ap)

% 判断输入图像是否大小一致

[z1 s1] = size(M1);

[z2 s2] = size(M2);

if (z1 ~= z2) | (s1 ~= s2)

error('输入图像大小不一致');

end;

switch(ap(1))

case 1,

% 绝对值最大融合准则

mm = (abs(M1)) > (abs(M2));

Y = (mm.*M1) + ((~mm).*M2);

case 2,

% 基于窗口系数加权平均融合准则

um = ap(2); th = .75;% 设定阈值

% 窗口区域加权平均能量

S1 = conv2(es2(M1.*M1, floor(um/2)), ones(um), 'valid');

S2 = conv2(es2(M2.*M2, floor(um/2)), ones(um), 'valid');

% 归一化相关度

MA = conv2(es2(M1.*M2, floor(um/2)), ones(um), 'valid');

MA = 2 * MA ./ (S1 + S2 + eps);

% s根据设定阈值选择融合系数

m1 = MA > th; m2 = S1 > S2;

w1 = (0.5 - 0.5*(1-MA) / (1-th));

Y = (~m1) .* ((m2.*M1) + ((~m2).*M2));

Y = Y + (m1 .* ((m2.*M1.*(1-w1))+((m2).*M2.*w1) + ((~m2).*M2.*(1-w1))+((~m2).*M1.*w1)));

case 3,

% 窗口系数绝对值选大融合准则

um = ap(2);

A1 = ordfilt2(abs(es2(M1, floor(um/2))), um*um, ones(um));

A2 = ordfilt2(abs(es2(M2, floor(um/2))), um*um, ones(um));

mm = (conv2((A1 > A2), ones(um), 'valid')) > floor(um*um/2);

Y = (mm.*M1) + ((~mm).*M2);

case 4,

% 系数最大融合准则

mm = M1 > M2;

Y = (mm.*M1) + ((~mm).*M2);

otherwise,

error('高频成分选择错误');

end;

附录2 融合图像性能评价程序

% 互信息函数

function mui = mutinf(M1, M2, F)

function mi = mutinf1(a, b)

a=double(a);

b=double(b);

[Ma,Na] = size(a);

[Mb,Nb] = size(b);

M=min(Ma,Mb);

N=min(Na,Nb);

% 初始化直方图数组

hab = zeros(256,256);

ha = zeros(1,256);

hb = zeros(1,256);

% 归一化

if max(max(a))~=min(min(a))

a = (a-min(min(a)))/(max(max(a))-min(min(a)));

else

a = zeros(M,N);

end

if max(max(b))-min(min(b))

b = (b-min(min(b)))/(max(max(b))-min(min(b)));

else

b = zeros(M,N);

end

a = double(int16(a*255))+1;

b = double(int16(b*255))+1;

% 统计直方图

for i=1:M

for j=1:N

indexx = a(i,j);

indexy = b(i,j) ;

hab(indexx,indexy) = hab(indexx,indexy)+1;% 联合直方图

ha(indexx) = ha(indexx)+1;% M1直方图

hb(indexy) = hb(indexy)+1;% M2直方图

end

end

% 联合信息熵

hsum = sum(sum(hab));

index = find(hab~=0);

p = hab/hsum;

Hab = sum(sum(-p(index).*log(p(index))));

% M1信息熵

hsum = sum(sum(ha));

index = find(ha~=0);

p = ha/hsum;

Ha = sum(sum(-p(index).*log(p(index))));

% M2信息熵

hsum = sum(sum(hb));

index = find(hb~=0);

p = hb/hsum;

Hb = sum(sum(-p(index).*log(p(index))));

% M1与M2互信息

mi = Ha+Hb-Hab;

% M1与M2归一化互信息

mi1 = hab/(Ha+Hb);

end

mui=(mutinf1(M1, F)+mutinf1(M2, F));

end

% 均方差函数

function img_var= variance(M1,M2,img)

M1 = double(M1);

M2 = double(M2);

img = double(img);

[r, c] = size(img);

img_var = (sqrt(sum(sum((M1 - img).^2)) / (r * c ))+ sqrt(sum(sum((M2 - img).^2)) / (r * c )))/2;

end

% 峰值信噪比函数

function psnr = psnr(M1,M2,img)

M1 = double(M1);

M2 = double(M2);

img = double(img);

[r, c] = size(img);

psnr = 10.*(log((r * c )/sum(sum((M1 - img).^2)))+log((r * c )/sum(sum((M2 -

img).^2))))/2;

end

% 平均梯度函数

function avg_gra = avg_gradient(img)

if nargin == 1 % 判断输入变量个数

img = double(img);

[r,c,b] = size(img);

dx = 1;

dy = 1;

for k = 1 : b

band = img(:,:,k);

% 沿x与y方向差分

[dzdx,dzdy] = gradient(band,dx,dy);

s = sqrt((dzdx .^ 2 + dzdy .^2) ./ 2);

g(k) = sum(sum(s)) / ((r - 1) * (c - 1));

end

avg_gra = mean(g);

else

error('Wrong number of input!');

end

% 边缘保持度函数

function EDG_pr=edge_preservation(M1, M2, img)

% 源图像与融合图像的梯度幅值与相角

[GA,ALPHA_A,flag1]=SOBEL_EDGE(M1);

[GB,ALPHA_B,flag2]=SOBEL_EDGE(M2);

[GF,ALPHA_F,flag3]=SOBEL_EDGE(img);

% 各像素边缘信息保留程度

QAF=EDGE_PRE_VAL(GA,ALPHA_A,GF,ALPHA_F); QBF=EDGE_PRE_VAL(GB,ALPHA_B,GF,ALPHA_F); % 源图像与融合图像边缘保持度

EDG_pr=NOR_WGT_VAL(QAF,GA,QBF,GB);

end

function [G,ALPHA,flag]=SOBEL_EDGE(img)

% Sobel滤波梯度幅值与相角求取函数

[row,col]=size(img);

plate_H=[-1,-2,-1;0,0,0;-1,-2,-1];

plate_V=[-1,0,-1;-2,0,-2;-1,0,-1];

SX=conv2(img,plate_H,'same');

SY=conv2(img,plate_V,'same');

flag=0;

G=zeros(row,col);

ALPHA=zeros(row,col);

for(i=1:row)

for(j=1:col)

if(SX(i,j)~=0)

temp=SX(i,j)*SX(i,j)+SY(i,j)*SY(i,j);

G(i,j)=sqrt(temp);

ALPHA(i,j)=atan(SY(i,j)/SX(i,j));

if(ALPHA(i,j)<0|G(i,j)==0)

flag=flag+1;

end

end

end

end

end

function QAF=EDGE_PRE_VAL(GA,ALPHA_A,GF,ALPHA_F) % 边缘信息保留程度求取函数

[row,col]=size(GA);

GAF=zeros(row,col);

AAF=zeros(row,col);

%求取源图像与融合图像的相对幅值和相角

for(i=1:row)

for(j=1:col)

if(GA(i,j)>GF(i,j))

GAF(i,j)=GF(i,j)/GA(i,j);

else

GAF(i,j)=GA(i,j)/GF(i,j);

end

end

end

ALPHA_DIF=ALPHA_A-ALPHA_F;

ALPHA_DIF=(abs(ALPHA_DIF)*2)/pi;

AAF=1-ALPHA_DIF;

% 可调节参数取值

GAMA_G=0.9994;

KG=-15;

DELTA_G=0.5;

GAMA_A=0.9879;

KA=-22;

DELTA_A=0.8;

TEMP=exp(KG*(GAF-DELTA_G));

QAF_G=GAMA_G./(1+TEMP);

TEMP=exp(KA*(AAF-DELTA_A));

QAF_ALPHA=GAMA_A./(1+TEMP);

QAF=QAF_G.*QAF_ALPHA;

size(QAF);

sum(sum(QAF));

end

function NWP=NOR_WGT_VAL(QAF,GA,QBF,GB)

% 边缘保持度求取函数

NWP=0;

temp1=(QAF.*GA+QBF.*GB);

temp2=sum(sum(temp1)); temp3=sum(sum(GA+GB)); NWP=temp2/temp3;

end