多焦点图像融合算法附matlab代码

clear all; close all; for jjj = 1:14 name = num2str(jjj,'%02d'); input1=imread(['Images\source' name '_1.tif']); input2=imread(['Images\source' name '_2.tif']); input1 = double(input1); input2 = double(input2); hei2 = size(input1,1); wid2 = size(input1,2); hei = ceil(size(input1,1) /2); wid = ceil(size(input1,2) /2); lvl = size(input1,3); N = 7; win = ones(N,N); thold = conv2(ones(hei,wid),win,'same') *0.4; thold2 = thold * 2; thold3 = conv2(ones(ceil(hei/N),ceil(wid/N)),win,'same') *0.5; C1 = zeros(hei,wid); C2 = zeros(hei,wid); sum_xh12 = zeros(hei-1,wid-1); sum_xh22 = zeros(hei-1,wid-1); map1 = zeros(hei,wid,lvl); map2 = zeros(hei,wid,lvl); map_yat = zeros(hei,wid,lvl); wav = 'haar'; my_fun = @(x) round(x.^2); tic thetime = cputime; for iii = 1:10 [xa1,xh1,xv1,xd1,xxa1,xxh1,xxv1,xxd1]=my_dwt2(input1,wav); [xa2,xh2,xv2,xd2,xxa2,xxh2,xxv2,xxd2]=my_dwt2(input2,wav); for i = 1:lvl C1 = C1 + my_fun(xh1(:,:,i)) + my_fun(xv1(:,:,i)) + my_fun(xd1(:,:,i)); C2 = C2 + my_fun(xh2(:,:,i)) + my_fun(xv2(:,:,i)) + my_fun(xd2(:,:,i)); C1(1:end-1,1:end-1) = C1(1:end-1,1:end-1) + my_fun(xxh1(:,:,i)) + my_fun(xxv1(:,:,i)) + my_fun(xxd1(:,:,i)); C2(1:end-1,1:end-1) = C2(1:end-1,1:end-1) + my_fun(xxh2(:,:,i)) + my_fun(xxv2(:,:,i)) + my_fun(xxd2(:,:,i)); end delta = sign(conv2(C1-C2, win, 'same')); temp = zeros(ceil(hei/N),ceil(wid/N)); for i = 1:N:hei-N+1 for j = 1:N:wid-N+1 temp((i+N-1)/N,(j+N-1)/N) = sum(sum(delta(i:i+N-1,j:j+N-1))); end end if sum(abs(temp(:,end))) == 0 temp(:,end) = temp(:,end-1); end if sum(abs(temp(end,:))) == 0 temp(end,:) = temp(end-1,:); end temp(temp<thold3&temp>-thold3) = 0; temp(temp>=thold3) = 1; temp(temp<=-thold3) = -1; temp1 = conv2(double(temp==1),ones(3,3),'same') ./ conv2(ones(size(temp)),ones(3,3),'same'); temp2 = conv2(double(temp==-1),ones(3,3),'same') ./ conv2(ones(size(temp)),ones(3,3),'same'); temp = zeros(ceil(hei/N),ceil(wid/N)); temp(temp1==1) = 1; temp(temp2==1) = -1; tempt = temp == 1; tempt = 1-bwareaopen(1-tempt,floor(hei*wid/N^2/100)); temp(tempt==1) = 1; tempt = temp == -1; tempt = 1-bwareaopen(1-tempt,floor(hei*wid/N^2/100)); temp(tempt==1) = -1; temp = imresize(temp,N,'nearest'); temp = temp(1:hei,1:wid,:); W = conv2(delta, win, 'same'); map1 = max(sign(W-thold),0); map2 = max(sign(-thold-W),0); map1 = temp==1 | (map1 & temp==0); map2 = temp==-1 | (map2 & temp==0); map3 = 1 - map1 - map2; map_yat = (thold-W)./thold2; if lvl > 1 map_yat = repmat(map_yat,1,1,lvl); map1 = repmat(map1,1,1,lvl); map2 = repmat(map2,1,1,lvl); map3 = repmat(map3,1,1,lvl); end map_yht = abs(xh1) > abs(xh2); map_yvt = abs(xv1) > abs(xv2); map_ydt = abs(xd1) > abs(xd2); yat = xa2.*map_yat + xa1.*(1-map_yat); yht = xh1.*map_yht + xh2.*(1-map_yht); yvt = xv1.*map_yvt + xv2.*(1-map_yvt); ydt = xd1.*map_ydt + xd2.*(1-map_ydt); ya = xa1.*map1 + xa2.*map2 + yat.*map3; yh = xh1.*map1 + xh2.*map2 + yht.*map3; yv = xv1.*map1 + xv2.*map2 + yvt.*map3; yd = xd1.*map1 + xd2.*map2 + ydt.*map3; output = my_idwt2(ya,yh,yv,yd,wav); output = output(1:hei2,1:wid2,:); end toc % imwrite(output(1:hei2,1:wid2,1:lvl)/255, ['R:\source' name '_SSWT.png']); end