% Copyright (c) 2017 Indian Institute of Science, Bangalore, India.
% All rights reserved.
% 
% Permission to use, copy, or modify this software and its documentation
% for educational and research purposes only and without fee is
% hereby granted, provided that this copyright notice and the original 
% authors' names appear on all copies and supporting documentation.
% This program shall not be used, rewritten, or adapted as the basis
% of a commercial software or hardware product without first obtaining
% permission of the authors. The authors make no representations about the
% suitability of this software for any purpose. It is provided "as is"
% without express or implied warranty.
% 
% The following papers are to be cited in the bibliography whenever
% the software is used as:
% Hisham Rahman, Rajiv Soundararajan, and R. Venkatesh Babu,
% "Evaluating Multiexposure Fusion Using Image Information",
% Signal Processing Letters, in print, DOI: 10.1109/LSP.2017.2752233.
% 
% The attached code computes a single quality score given a test image
% and set of over/under-exposed reference images.
% 
% For questions and/or comments please email
% Hisham Rahman (hishamrahman@gmail.com).

% Function computes the mssim q-map of fused image with all the distorted
% image and mix their values using pyrW to get one single ms-ssim qmap.
% Then compute overall ssim score by weighting with meanWghtArray.
% Fused and distorted images are in 0-255 double precision.
function [mixed_ssim_val] = ssim_mixing_StPyr(imgTest, imgRef, pyrW, meanWghtArray)
    [H, W, totRef] = size(imgRef);

    % Get the ms-ssim qmap of test image with respect to each ref image.  
    for j = 1:totRef
        [test_cs_array] = msssim_qmap_StPyr(imgTest, imgRef(:, :, j));
        if(j == 1)
            nLev = length(test_cs_array);
            cs_qmap_array = cell(1, nLev);
            final_qmap_array = cell(1, nLev); % To save the final fused qmap.
        end
        for i= 1:nLev
            if(j == 1)
                cs_qmap_array{i} = zeros(size(test_cs_array{i}, 1), size(test_cs_array{i}, 2), totRef);
                final_qmap_array{i} = zeros(size(test_cs_array{i}));
            end
            cs_qmap_array{i}(:, :, j) = test_cs_array{i};
        end  
    end
    
    % For each level, mix the ssim qmap according to the weights.
    for i = 1:nLev
        for j = 1:totRef
            final_qmap_array{i} = final_qmap_array{i} + (pyrW{i}(:,:,j).*cs_qmap_array{i}(:,:,j));
        end
    end

    % Compute the final overall ms-ssim score.
    m_score = zeros(1, nLev);
    for i = 1:nLev
        m_score(1, i) = sum(sum(final_qmap_array{i}.*meanWghtArray{i}))/sum(sum(meanWghtArray{i}));  
    end
    % Original MS-SSIM weights divided into horizontal and vertical
    % subbands.
    weight = [0.1333/2 0.1333/2, 0.2363/2 0.2363/2, 0.3001/2 0.3001/2, 0.2856/2 0.2856/2, 0.0448/2 0.0448/2];
    m_score = m_score.^weight;
    mixed_ssim_val = prod(m_score);    
end





