% Inputs:
%   IM1, IM2    Equally sized images to be compared (0-255 doubles)
%   LEVEL		Level or "Height" of the pyramid decomposition to use
%                   Note: currently, this also determines the level that is
%                         used when computing cssim -> can't use the
%                         baseband <TO CHANGE>
%                   Note: higher pyramids add more low-frequency subbands
%   NUMOR       Number of orientations to use in the pyramid decomposition
%   GUARDB      Size of the boundry region to exclude
%   WINSIZE     Size of the sliding window (optional, default 7)
%   WTS         Weight vector [HP BP(1...LEVEL) LP], length LEVEL+2 (opt)
%
% Outputs:
%   CSSIM_WTD   Overall metric value with weighted subbands
%   CSSIM       Overall mean metric value of passband only
%   LEV_CSSIM   Vector of metric values for each level (from HP to BB)
%   BAND_CSSIM  Matrix of metric values for each orientation
%   CSSIM_MAP   RxCxOR Matrix of individual metric values (for RxC images)

% Original Author: Zhou Wang
% Current Changes Author: Alan Brooks
%
% 16-Mar-2005   slight modification to also pass back cssim_map's
% 26-Sep-2005   added winsize parameter
% 13-Oct-2005   added ability to access baseband
% 27-Oct-2005   removed mean from cssim_lp computations (since equ assumes)
%               added weights (need sub-function)
%  6-Jun-2006   fixed function name in main function (cosmetic)

function [cssim_wtd, cssim, lev_cssim, band_cssim, cssim_map, ... 
    cssim_lp, cssim_map_lp, cssim_hp, cssim_map_hp, ...
    ssim_lp, ssim_map_lp] = ...
    cssim_index_multi(im1, im2, level, numOr, guardb, winsize, wts)

% build the steerable pyramids
[pyr1, pind] = buildSCFpyr(im1, level, numOr-1);
[pyr2, pind] = buildSCFpyr(im2, level, numOr-1);

%figure,showSpyr(abs(pyr2),pind);

% initialize
if ~exist('wts','var')
    % weight the subbands evenly if no weights given
    numWts = level + 2;
    wts = ones(numWts,1)/numWts;
else
    wts = wts/sum(wts);
end
if ~exist('winsize','var')
    nw = 7;
else
    nw = winsize;
end
window = ones(nw);
window = window./sum(sum(window));
C = 0.03;
band_cssim = zeros(level,numOr);
[sx,sy] = size(im1);
%gb = guardb;

% bandpass subbands at each level
for lev=1:level %(high to low frequency)
    gb = round(guardb/(2^(lev-1)));
    for i=1:numOr
        bandind = i+(lev-1)*numOr+1; 
        [corr_band, varr_band] = ...
            calcCorrAndVarr(pyr1,pyr2,pind,bandind,gb,window);
        map = (2*abs(corr_band) + C)./(varr_band + C);
        cssim_map{lev,i} = map;
        band_cssim(lev,i) = mean2(cssim_map{lev,i});
    end % loop through orientations
end % loop over levels

% highpass and lowpass subbands
    % LP (baseband)
    gb = round(guardb/(2^(level-1)));
    bandind = size(pind,1);
    [corr_band, varr_band] = ...
        calcCorrAndVarr(pyr1,pyr2,pind,bandind,gb,window,1);
    map = (2*abs(corr_band) + C)./(varr_band + C);
    cssim_map_lp = map;
    
        % === EXTRA TEST CODE CALCULATING SSIM ON LP BAND ===
            % do same as W = fspecial('gaussian',srN,sigma)
        ndx = [1:nw]-nw/2-.5;
        sigma = 1.5;
        wts2 = ( (exp(-ndx.^2./(2*sigma^2)))'* ...
                (exp(-ndx.^2./(2*sigma^2)))   ) /(2*pi*sigma^2);
        wts2 = wts2./sum(wts2(:));
        im1_bb = pyrLow(pyr1, pind);
        im2_bb = pyrLow(pyr2, pind);
        maxLev = max([im1_bb(:); im2_bb(:)]);
        [ssim_lp ssim_map_lp] = ...
            ssim_index(im1_bb, im2_bb, [.01 .03], wts2, maxLev);
        % === END TEST CODE ===
        
    % HP
    bandind = 1;
    gb = guardb;
    [corr_band, varr_band] = ...
        calcCorrAndVarr(pyr1,pyr2,pind,bandind,gb,window);
    map = (2*abs(corr_band) + C)./(varr_band + C);
    cssim_map_hp = map;   

% average over the orientations
lev_cssim = mean(band_cssim,2);

% add in the hp and lp (high-to-low frequnecy order)
cssim_lp = mean(cssim_map_lp(:));
cssim_hp = mean(cssim_map_hp(:));
lev_cssim = [cssim_hp lev_cssim' cssim_lp]';

% compute the weighted cssim including baseband and highpass
cssim_wtd = sum(lev_cssim.*wts(:));

% compute the original definition using only lowest freq passband
cssim = lev_cssim(level+1);




% === utility functions ===
function [corr_band, varr_band] = ...
    calcCorrAndVarr(pyr1,pyr2,pind,bandind,gb,window,remove_mean)
if ~exist('remove_mean','var'), remove_mean = 0; end
band1 = pyrBand(pyr1, pind, bandind);
band2 = pyrBand(pyr2, pind, bandind);
band1 = band1(gb+1:end-gb, gb+1:end-gb);
band2 = band2(gb+1:end-gb, gb+1:end-gb);
if remove_mean
    band1 = band1 - mean(band1(:)); % special for LP, remove mean
    band2 = band2 - mean(band2(:));
end
corr = band1.*conj(band2);
varr = abs(band1).^2 + abs(band2).^2;
corr_band = filter2(window, corr, 'valid');
varr_band = filter2(window, varr, 'valid');