Adaptive multi-focus image fusion using a wavelet-based statistical sharpness measure

by

 

17 Apr 2012 (Updated )

Adaptive multi-focus image fusion using a wavelet-based statistical sharpness measure

SP_2012_main.m
function SP_2012_main
% This is a demo program of the paper J. Tian and L. Chen Adaptive
% multi-focus image fusion using a wavelet-based statistical sharpness measure,
% Signal Processing, Vol. 92, No. 9, Sep. 2012, pp. 2137-2146.

close all; clear all; clc;

MA = double(imread('image1_A.bmp'))./255;
MB = double(imread('image1_B.bmp'))./255;

% the input image MUST be gray-scale image.
% The image size MUST be power of 2, such as 128*128, 256*256.

MF = func_SP_2012(MA, MB);
MF = round(MF.*255);
MF(MF<0) = 0;
MF(MF>255) = 255;
        
fprintf('Mutural Information: %.2f\n', func_fusion_evaluate_mutural_information(MA.*255, MB.*255, MF, 256));

%--------------------------------------------------------------------------
%----------------------------Inner Function -------------------------------
%--------------------------------------------------------------------------

function MF = func_SP_2012(MA, MB)
% This is a demo program of the paper 

[nRow, nCol] = size(MA);
wbase = 'Daubechies';
%wbase = 'Haar';
mom = 4;
dwt_level = 1;
L = log2(nRow)-dwt_level;
qmf = func_MakeONFilter(wbase,mom);

window_size = 5;

% apply DWT
WA  = func_FWT2_PO(MA, L, qmf);
[t1,t2] = func_dyad2HH(L);
HHA = WA(t1,t2);
HHA_extend = padarray(HHA,[(window_size-1)/2 (window_size-1)/2],'symmetric','both');
[t1,t2] = func_dyad2HL(L);
HLA = WA(t1,t2);
HLA_extend = padarray(HLA,[(window_size-1)/2 (window_size-1)/2],'symmetric','both');
[t1,t2] = func_dyad2LH(L);
LHA = WA(t1,t2);
LHA_extend = padarray(LHA,[(window_size-1)/2 (window_size-1)/2],'symmetric','both');
[t1,t2] = func_dyad2LL(L);
LLA = WA(t1,t2);
LLA_extend = padarray(LLA,[(window_size-1)/2 (window_size-1)/2],'symmetric','both');

WB  = func_FWT2_PO(MB, L, qmf);
[t1,t2] = func_dyad2HH(L);
HHB = WB(t1,t2);
HHB_extend = padarray(HHB,[(window_size-1)/2 (window_size-1)/2],'symmetric','both');
[t1,t2] = func_dyad2HL(L);
HLB = WB(t1,t2);
HLB_extend = padarray(HLB,[(window_size-1)/2 (window_size-1)/2],'symmetric','both');
[t1,t2] = func_dyad2LH(L);
LHB = WB(t1,t2);
LHB_extend = padarray(LHB,[(window_size-1)/2 (window_size-1)/2],'symmetric','both');
[t1,t2] = func_dyad2LL(L);
LLB = WB(t1,t2);
LLB_extend = padarray(LLB,[(window_size-1)/2 (window_size-1)/2],'symmetric','both');

CA = (im2col(HHA_extend,[window_size window_size],'sliding')).^2+ ...
    (im2col(LHA_extend,[window_size window_size],'sliding')).^2+ ...
    (im2col(HLA_extend,[window_size window_size],'sliding')).^2;
CB = (im2col(HHB_extend,[window_size window_size],'sliding')).^2+ ...
    (im2col(LHB_extend,[window_size window_size],'sliding')).^2+ ...
    (im2col(HLB_extend,[window_size window_size],'sliding')).^2;
CA_entropy = sum(abs(CA),1);
CB_entropy = sum(abs(CB),1);
judge1 = (CA_entropy>CB_entropy);
LLF = LLA(:).*judge1' + LLB(:).*(1-judge1)';

clear CA CB CA_mean CB_mean CA_std CB_std judge1 judge2

% high-frequence subband
LHF = func_fusion_high_subband(LHA, LHB);
HLF = func_fusion_high_subband(HLA, HLB);
HHF = func_fusion_high_subband(HHA, HHB);

% inverse DWT
WF = zeros(size(WA));
[t1,t2] = func_dyad2HH(L);
WF(t1,t2) = reshape(HHF, size(HHA));
[t1,t2] = func_dyad2HL(L);
WF(t1,t2) = reshape(HLF, size(HLA));
[t1,t2] = func_dyad2LH(L);
WF(t1,t2) = reshape(LHF, size(LHA));
[t1,t2] = func_dyad2LL(L);
WF(t1,t2) = reshape(LLF, size(LLA));

MF  = func_IWT2_PO(WF, L, qmf);

%--------------------------------------------------------------------------
%----------------------------Inner Function -------------------------------
%--------------------------------------------------------------------------
function subbandF = func_fusion_high_subband(subbandA, subbandB)
window_size = 7;
subbandA_extend = padarray(subbandA,[(window_size-1)/2 (window_size-1)/2],'symmetric','both');
subbandB_extend = padarray(subbandB,[(window_size-1)/2 (window_size-1)/2],'symmetric','both');
CA = im2col(subbandA_extend,[window_size window_size],'sliding');
CB = im2col(subbandB_extend,[window_size window_size],'sliding');

% A binary thresholding is applied here for fast estimation of distribution
% parameters.
for ii=1:size(CA,2)
    temp = CA(:,ii);
    if max(max(abs(temp)))==0
        CA_std(ii) = 0;
    else
        T = graythresh(abs(temp)./max(maax(abs(temp)))).*max(max(abs(temp)));
        temp(abs(temp)<T)=[];
        if numel(temp)~=0
            CA_std(ii) = mean2(abs(temp));
        else
            CA_std(ii) = 0;
        end        
    end
    
    temp = CB(:,ii);
    if max(max(abs(temp)))==0
        CB_std(ii) = 0;
    else    
        T = graythresh(abs(temp)./max(max(abs(temp)))).*max(max(abs(temp)));
        temp(abs(temp)<T)=[];
        if numel(temp)~=0
            CB_std(ii) = mean2(abs(temp));
        else
            CB_std(ii) = 0;
        end
    end
end

judge1 = (CA_std>=CB_std);
judge2 = (CA_std<CB_std);

subbandF = subbandA(:).*judge1' + subbandB(:).*judge2';

%--------------------------------------------------------------------------
%----------------------------Inner Function -------------------------------
%--------------------------------------------------------------------------
function mutural_informationR = func_fusion_evaluate_mutural_information(grey_matrixA,grey_matrixB,grey_matrixF,grey_level)
% Performance evaluation using mutural information.
% Qu Guihong, Zhang Dali, Yan Pingfan. Information measure for performance of image fusion. Electronics Letters, 2002,38(7): 313-15.
% compute mutural information of the image
% grey_matrixA, grey_matrixB, grey_matrixF are grey values of imageA,imageB and fusion image
% grey_level is the grayscale degree of image
% please set grey_level=256
% ---------
% Author:  Qu Xiao-Bo    <quxiaobo [at] xmu.edu.cn>    June 26, 2009
%          Postal address:
% Rom 509, Scientific Research Building # 2,Haiyun Campus, Xiamen University,Xiamen,Fujian, P. R. China, 361005
% Website: http://quxiaobo.go.8866.org

HA=entropy_fusion(grey_matrixA,grey_level);
HB=entropy_fusion(grey_matrixB,grey_level);
HF=entropy_fusion(grey_matrixF,grey_level);
HFA=Hab(grey_matrixF,grey_matrixA,grey_level);
HFB=Hab(grey_matrixF,grey_matrixB,grey_level);
MIFA=HA+HF-HFA;
MIFB=HB+HF-HFB;
mutural_informationR=MIFA+MIFB;

%--------------------------------------------------------------------------
%----------------------------Inner Function -------------------------------
%--------------------------------------------------------------------------
function entropyR=entropy_fusion(grey_matrix,grey_level)
% Author:  Qu Xiao-Bo    <quxiaobo [at] xmu.edu.cn>    June 26, 2009
%          Postal address:
% Rom 509, Scientific Research Building # 2,Haiyun Campus, Xiamen University,Xiamen,Fujian, P. R. China, 361005
% Website: http://quxiaobo.go.8866.org
[row,column]=size(grey_matrix);
total=row*column;
% grey_level=256 

counter=zeros(1,grey_level);
grey_matrix=grey_matrix+1;
for i=1:row
    for j=1:column
        indexx= grey_matrix(i,j);
        counter(indexx)=counter(indexx)+1;
    end
end
total= sum(counter(:));
index = find(counter~=0);
 p = counter/total;
entropyR= sum(sum(-p(index).*log2(p(index))));

%--------------------------------------------------------------------------
%----------------------------Inner Function -------------------------------
%--------------------------------------------------------------------------
function HabR=Hab(grey_matrixA,grey_matrixB,grey_level)
% HabR=Hab(grey_matrixA,grey_matrixB,grey_level)
% compute mutural information of the image
% grey_matrixA , grey_matrixB,grey_matrixF are grey values of imageA,imageB and fusion image
% grey_level is the grayscale degree of image
% --------- 
% ---------
% Author:  Qu Xiao-bo    <quxiaobo429@163.com>    May 7,2006
%          Postal address:
%          Xiamen University, Department of Communication Engineering
%          Xiamen, Fujian, P. R. China, 361005  
[row,column]=size(grey_matrixA);
counter = zeros(256,256);
grey_matrixA=grey_matrixA+1;
grey_matrixB=grey_matrixB+1;
for i=1:row
    for j=1:column
        indexx = grey_matrixA(i,j);
        indexy = grey_matrixB(i,j);
        counter(indexx,indexy) = counter(indexx,indexy)+1;%ֱͼ
    end
end
total= sum(counter(:));
index = find(counter~=0);
p = counter/total;
HabR = sum(sum(-p(index).*log2(p(index))));        

Contact us