DctVarFusion(im1, im2)

function [fusedDctVarCv, fusedDctVar] = DctVarFusion(im1, im2)

% "DctVarFusion" is the simulation of following multifocus image fusion methods:
% 
% (1) DCT+Variance
% (2) DCT+Variance+CV
% 
% proposed in:
% 
% M.B.A. Haghighat, A. Aghagolzadeh, H. Seyedarabi:
% "Multi-Focus Image Fusion for Visual Sensor Networks in DCT Domain"
% Computers and Electrical Engineering, vol. 37, no. 5, pp. 789-797, Sep. 2011
% http://dx.doi.org/10.1016/j.compeleceng.2011.04.016
% 
% 
% Inputs:
%       im1:   First source image 
%       im2:   Second source image
%               
% Outputs:
%       fusedDctVar:    Fused image as the result of "DCT+Variance" method
%       fusedDctVarCV:  Fused image as the result of "DCT+Variance+CV" method
% 
% 
% (C) Mohammad Bagher Akbari Haghighat, University of Tabriz


if (nargin ~= 2)    % check correct number of arguments
    error('There should be two input images!')
end

% Get input image size
[m,n] = size(im1);
fusedDctVar = zeros(m,n);
fusedDctVarCv = zeros(m,n);
cvMap = zeros(m/8,n/8);	% Consistency verification index map

% Level shifting
im1 = double(im1)-128;
im2 = double(im2)-128;

% Devide source images into 8*8 blocks and perform the fusion
for i = 1:m/8
    for j = 1:n/8
        
        im1Sub = im1(8*i-7:8*i,8*j-7:8*j);
        im2Sub = im2(8*i-7:8*i,8*j-7:8*j);
        
        % Compute the 2-D DCT of 8*8 blocks
        im1SubDct = dct2(im1Sub);
        im2SubDct = dct2(im2Sub);
        
        % Calculate normalized transform coefficients
        im1Norm = im1SubDct ./ 8;
        im2Norm = im2SubDct ./ 8;
        
        % Mean value of 8*8 block of images (Measure for surrounding lumminance)
        im1Mean = im1Norm(1,1);
        im2Mean = im2Norm(1,1);
        
        % Variance of 8*8 block of images
        im1Var = sum(sum(im1Norm.^2)) - im1Mean.^2;
        im2Var = sum(sum(im2Norm.^2)) - im2Mean.^2;
        
        % Fusion
        if im1Var > im2Var
            dctVarSub = im1SubDct;
            cvMap(i,j) = -1;	% Consistency verification index
        else
            dctVarSub = im2SubDct;
            cvMap(i,j) = +1;    % Consistency verification index
        end
        
        % Compute the 2-D inverse DCT of 8*8 blocks and construct fused image
        fusedDctVar(8*i-7:8*i,8*j-7:8*j) = idct2(dctVarSub); % DCT+Variance method
        
    end
end

% Concistency verification using a Majority Filter
fi = ones(7)/49; % Filter kernel 7*7
cvMapFiltered = imfilter(cvMap, fi, 'symmetric'); % Filtered index map
cvMapFiltered = imfilter(cvMapFiltered, fi, 'symmetric');


for i = 1:m/8
    for j = 1:n/8
        % DCT+Variance+CV method
        if cvMapFiltered(i,j) <= 0
            fusedDctVarCv(8*i-7:8*i,8*j-7:8*j) = im1(8*i-7:8*i,8*j-7:8*j);
        else
            fusedDctVarCv(8*i-7:8*i,8*j-7:8*j) = im2(8*i-7:8*i,8*j-7:8*j);
        end
        
    end
end


% Inverse level shifting: 
im1 = uint8(double(im1)+128);
im2 = uint8(double(im2)+128);
fusedDctVar = uint8(double(fusedDctVar)+128);
fusedDctVarCv = uint8(double(fusedDctVarCv)+128);

% Show table
subplot(2,2,1), imshow(im1), title('Focus on right');
subplot(2,2,2), imshow(im2), title('Focus on left');
subplot(2,2,3), imshow(fusedDctVar), title('DCT+Variance fusion result');
subplot(2,2,4), imshow(fusedDctVarCv), title('DCT+Variance+CV fusion result');