2022年2022年金字塔变换的图像融合代码 .pdf

1 附录 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 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 1 页，共 11 页 - - - - - - - - - 2 M1T = conv2(conv2(es2(undec2(M1), 2), 2*w, valid), 2*w, valid); M1 = M1T + Ei1; % 选择图像有效区域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) .* (Ei1+1); M1 = M1(1:zl(i1),1:sl(i1); end; Y = M1; end; 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 2 页，共 11 页 - - - - - - - - - 3 % 梯度金字塔融合函数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); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 3 页，共 11 页 - - - - - - - - - 4 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 + Ei1; 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) 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 4 页，共 11 页 - - - - - - - - - 5 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); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 5 页，共 11 页 - - - - - - - - - 6 % 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; 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 6 页，共 11 页 - - - - - - - - - 7 附录 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; 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 7 页，共 11 页 - - - - - - - - - 8 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; 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 8 页，共 11 页 - - - - - - - - - 9 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)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); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 10 页，共 11 页 - - - - - - - - - 11 temp2=sum(sum(temp1); temp3=sum(sum(GA+GB); NWP=temp2/temp3; end 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 11 页，共 11 页 - - - - - - - - -