% Run our sequence-based method SeqConv
% 
% S ... pairwise similarity matrix
% n ... approx. seq length (n ~ past+1+future)
%
%
%   =====================================================================
%   Copyright (C) 2021  Stefan Schubert, stefan.schubert@etit.tu-chemnitz.de
%   
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%   
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%   
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%   =====================================================================
% 
function SS = SeqConv(S, n)
    % setup possible velocities
    v = 0.8:0.1:1.2;
    
    % initialize output
    SS = -Inf * ones(size(S,1), size(S,2));
    
    % run over all convoluational kernels
    prevFilter = [];
    for i = 1:numel(v)
        % compute filter / convoluational kernel
        H =  getFilter(n,v(i));
        
        % if filter is not identical to the previous one, run convoluation
        if ~isSame(H, prevFilter)            
            % correction C for border area
            C = imfilter(ones(size(S,1), size(S,2)), H / sum(H(:)), 'same');
            
            % convolve
            SS = max(SS, 1./C .* imfilter(S, H, 'same'));
            
            % store filter
            prevFilter = H;
        end
    end
end

% Generate a [~v*n, n] mask H with a line of 1 with slope v. Dimensions of
% H are odd, the line goes through the center element.
% n ... number of ones
% v ... slope
function H = getFilter(n,v)
    assert(v>0);
    assert(v>0.5);

    nh = floor(n/2);
    
    % get bottom right part
    x = 1:nh;
    y = round(x*v);
    hh = max(y);
    wh = max(x);
    
    idx = sub2ind([hh,wh], y, x);
    Hh = zeros(hh,wh);
    Hh(idx) = 1;
    
    % combine: top-left is Hh, centre is 1, bottom-right is Hh
    H = [Hh, zeros(hh, wh+1);
        zeros(1,wh), 1, zeros(1,wh);
        zeros(hh,wh+1), Hh];
end

% Are matrices X and Y identical?
function r = isSame(X,Y)
    if size(X,1) == size(Y,1) && size(X,2) == size(Y,2) && sum(X(:)==Y(:)) == numel(X)
       r = true;
    else
       r = false;
    end
end
