% Temporal hysteresis pooling method
%
% K. Seshadrinathan and A. C. Bovik, “Temporal hysteresis model of time 
% varying subjective video quality,” in 2011 IEEE International Conference 
% on Acoustics, Speech and Signal Processing (ICASSP), 2011, pp. 1153 –1156.
%
% Only temporal pooling from Sec. 4.
%
% thtp = temppool_hysteresis(om,sps)
%
% om - vector of objective metric scores for each frame
% sps - scores per second
% tau - hysteresis coefficient: the number of seconds to consider

function thtp = temppool_hysteresis(om,sps,tau)

    if ~isvector(om)
        error('Input om must be a vector')
    end
    if ~isscalar(sps) || sps <= 0
        error('Input sps must be a positive scalar')
    end
    if ~isscalar(tau) || tau <= 0
        error('Input tau must be a positive scalar')
    end

    om=reshape(om,length(om),1); 
    
    alpha = 0.8; % Weight, recommended by the authors to be 0.8

    tau_values = tau*sps;

    x = zeros(length(om),1);
    x(1) = om(1);

    for i=2:length(om)  % apply hysteresis - eq. (7)
        t = [max(1,i-tau_values) i-1];
        x(i) = min( om( t(1):t(2) ) );
    end



    for i=1:length(om)
        t = [i min(i+tau_values,length(om))]; % Advancing window
        k = t(2) - t(1);                 % Number of elements being processed

        w = gausswin(k * 2);        % Create Gaussian window
        w = w(k+1:end);      % Take the descending half
        w = w/(sum(w));   % Normalize so that the function sums to 1

        v = sort(om(t(1):t(2)-1),'ascend');
        
        y(i) = sum(w .* v);    % weighting with Gaussian function

        if (i==length(om))
            y(i) = om(i);
        end
    end
    y=y'; % Transpose to column vector

    g = alpha*y + (1-alpha)*x;
    thtp = mean(g);


