view toolboxes/MIRtoolbox1.3.2/MIRToolbox/@mirsimatrix/mirsimatrix.m @ 0:e9a9cd732c1e tip

first hg version after svn

function varargout = mirsimatrix(orig,varargin)
%   m = mirsimatrix(x) computes the similarity matrix resulting from the 
%       mutual comparison between each possible frame analysis in x.
%           By default, x is the spectrum of the frame decomposition.
%           But it can be any other frame analysis.
%   Optional argument:
%       mirsimatrix(...,'Distance',f) specifies the name of a distance
%           function, from those proposed in the Statistics Toolbox
%               (help pdist).
%           default value: f = 'cosine'
%       mirsimatrix(...,'Dissimilarity') return the dissimilarity matrix,
%           which is the intermediary result before the computation of the 
%           actual similarity matrix. It shows the distance between each
%           possible frame analysis in x.
%       mirsimatrix(...,'Similarity',f) indicates the function f specifying
%           the way the distance values in the dissimilarity matrix are
%           transformed into similarity values.
%           Possible values:
%               f = 'oneminus' (default value) 
%                   corresponding to f(x) = 1-x
%               f = 'exponential'
%                   corresponding to f(x)= exp(-x)
%       mirsimatrix(...,'Width',w) or more simply dissimatrix(...,w) 
%           specifies the size of the diagonal bandwidth, in samples,
%           outside which the dissimilarity will not be computed.
%           if w = inf (default value), all the matrix will be computed.
%       mirsimatrix(...,'Horizontal') rotates the matrix 45 degrees in order to
%           make the first diagonal horizontal, and to restrict on the
%           diagonal bandwidth only.
%       mirsimatrix(...,'TimeLag') transforms the (non-rotated) matrix into
%           a time-lag matrix, making the first diagonal horizontal as well
%           (corresponding to the zero-lag line).
%
%       mirsimatrix(M,r) creates a mirsimatrix similarity matrix based on
%           the Matlab square matrix M, of frame rate r (in Hz.)
%               By default r = 20 Hz.
%           mirsimatrix(M,r,'Dissimilarity') creates instead a mirsimatrix
%               dissimilarity matrix.
%
%   Foote, J. & Cooper, M. (2003). Media Segmentation using Self-Similarity
% Decomposition,. In Proc. SPIE Storage and Retrieval for Multimedia
% Databases, Vol. 5021, pp. 167-75.

%       mirsimatrix(...,'Filter',10) filter along the diagonal of the matrix,
%           using a uniform moving average filter of size f. The result is
%           represented in a time-lag matrix.


        distance.key = 'Distance';
        distance.type = 'String';
        distance.default = 'cosine';
    option.distance = distance;

        simf.key = 'Similarity';
        simf.type = 'String';
        simf.default = 'oneminus';
        simf.keydefault = 'Similarity';        
        simf.when = 'After';
    option.simf = simf;

        dissim.key = 'Dissimilarity';
        dissim.type = 'Boolean';
        dissim.default = 0;
    option.dissim = dissim;
    
        K.key = 'Width';
        K.type = 'Integer';
        K.default = Inf;
    option.K = K;
    
        filt.key = 'Filter';
        filt.type = 'Integer';
        filt.default = 0;
        filt.when = 'After';
    option.filt = filt;

        view.type = 'String';
        view.default = 'Standard';
        view.choice = {'Standard','Horizontal','TimeLag'};
        view.when = 'After';
    option.view = view;
    
        rate.type = 'Integer';
        rate.position = 2;
        rate.default = 20;
    option.rate = rate;
    
specif.option = option;
specif.nochunk = 1;
varargout = mirfunction(@mirsimatrix,orig,varargin,nargout,specif,@init,@main);


function [x type] = init(x,option)
if isnumeric(x)
    m.diagwidth = Inf;
    m.view = 's';
    m.similarity = NaN;
    m.graph = {};
    m.branch = {};
    m = class(m,'mirsimatrix',mirdata);
    m = set(m,'Title','Dissimilarity matrix');
    fp = repmat(((1:size(x,1))-.5)/option.rate,[2,1]);
    x = set(m,'Data',{x},'Pos',[],...
              'FramePos',{{fp}},'Name',{inputname(1)});
else
    if not(isamir(x,'mirsimatrix'))
        if (isamir(x,'miraudio'))
            if isframed(x)
                x = mirspectrum(x);
            else
                x = mirspectrum(x,'Frame',0.05,1);
            end
        end
    end
end
type = 'mirsimatrix';


function m = main(orig,option,postoption)
if iscell(orig)
    orig = orig{1};
end
if isa(orig,'mirsimatrix')
    d = get(orig,'Data');
    for k = 1:length(d)
        if iscell(option) && isfield(option,'K') && option.K < orig.diagwidth
            nl = size(d{k},1);
            if strcmp(orig.view,'h')
                dl = (nl - option.K)/2;
                dk = d{k}(ceil(dl):nl-floor(dl),:);
            else
                [spA,spd] = spdiags(d{k},-ceil(option.K/2):ceil(option.K/2));
                dk = full(spdiags(spA,spd,nl,size(d{k},2)));
            end
            d{k} = dk;
            orig.diagwidth = option.K;
        end
    end
    m = set(orig,'Data',d);
else
    v = get(orig,'Data');
    d = cell(1,length(v));
    for k = 1:length(v)
        vk = v{k};
        if mirwaitbar
            handle = waitbar(0,'Computing dissimilarity matrix...');
        else
            handle = 0;
        end
        if 0 %iscell(vk)
            try
                vk = cell2mat(vk);
            end
        end
        if 0 %iscell(vk) %&& length(vk) > 1 %%% ATTENTION KL!!<<<<<<<<<<<<
            l = length(vk);
            dk = NaN(l,l);
            hK = floor(option.K/2);
            for i = 1:l
                if handle
                    waitbar(i/l,handle);
                end
                for j = max(1,i-hK):min(l,i+hK)
                    dk(i,j) = KL(vk{i},vk{j});
                end
            end
        else
            if not(iscell(vk))
                vk = {vk};
            end    
            for z = 1:length(vk)
                vz = vk{z};
                ll = size(vz,1);
                l = size(vz,2);
                nc = size(vz,3);
                if ll==1 && nc>1
                    vz = squeeze(vz)';
                    ll = nc;
                    nc = 1;
                end
                nd = size(vz,4);
                if not(isempty(postoption)) && ...
                    strcmpi(postoption.view,'TimeLag')
                    lK = ceil(option.K/2);
                    if isinf(lK)
                        lK = l;
                    end
                    dk{z} = NaN(lK,l,nc);
                else
                    dk{z} = NaN(l,l,nc);
                end
                for g = 1:nc
                    if nd == 1
                        vv = vz;
                    else
                        vv = zeros(ll*nd,l);
                        for h = 1:nd
                            if iscell(vz)
                                for i = 1:ll
                                    for j = 1:l
                                        vj = vz{i,j,g,h};
                                        if isempty(vj)
                                            vv((h-1)*ll+i,j) = NaN;
                                        else
                                            vv((h-1)*ll+i,j) = vj;
                                        end
                                    end
                                end
                            else
                                tic
                                vv((h-1)*ll+1:h*ll,:) = vz(:,:,g,h);
                                toc
                            end
                        end
                    end
                    if isinf(option.K) && not(strcmpi(postoption.view,'TimeLag'))
                        try
                            manually = 0;
                            dk{z}(:,:,g) = squareform(pdist(vv',option.distance));
                            if option.K < Inf
                                [spA,spd] = spdiags(dk{z},...
                                    -ceil(option.K/2):ceil(option.K/2));
                                dk{z}(:,:,g) = full(spdiags(spA,spd,size(dk,1),size(dk{z},2)));
                            end
                        catch
                            %err = lasterror;
                            %warning(err.message)
                            %disp('Statistics Toolbox does not seem to be
                            %installed. Recompute the distance matrix manually.');
                            manually = 1;
                        end
                    else
                        manually = 1;
                    end
                    if manually
                        disf = str2func(option.distance);
                        if strcmpi(option.distance,'cosine')
                            for i = 1:l
                                vv(:,i) = vv(:,i)/norm(vv(:,i));
                            end
                        end
                        if not(isempty(postoption)) && ...
                                strcmpi(postoption.view,'TimeLag')
                            lK = ceil(option.K/2);
                            for i = 1:l
                                if mirwaitbar && (mod(i,100) == 1 || i == l)
                                    waitbar(i/l,handle);
                                end
                                ij = min(i+lK-1,l);
                                dkij = disf(vv(:,i),vv(:,i:ij));
                                for j = 1:ij-i+1
                                    dk{z}(j,i+j-1,g) = dkij(j);
                                end
                            end
                        else
                            for i = 1:l
                                if mirwaitbar && (mod(i,100) == 1 || i == l)
                                    waitbar(i/l,handle);
                                end
                                j = min(i+option.K-1,l);
                                dkij = disf(vv(:,i),vv(:,i:j));
                                dk{z}(i,i:j,g) = dkij;
                                dk{z}(i:j,i,g) = dkij';
                            end
                        end
                    end
                end
            end
        end
        d{k} = dk;
        if handle
            delete(handle)
        end
    end
    m.diagwidth = option.K;
    if not(isempty(postoption)) && strcmpi(postoption.view,'TimeLag')
        m.view = 'l';
    else
        m.view = 's';
    end
    m.similarity = 0;
    m.graph = {};
    m.branch = {};
    m = class(m,'mirsimatrix',mirdata(orig));
    m = purgedata(m);
    m = set(m,'Title','Dissimilarity matrix');
    m = set(m,'Data',d,'Pos',[]);
end
if not(isempty(postoption))
    if strcmpi(m.view,'s')
        if strcmpi(postoption.view,'Horizontal')
            for k = 1:length(d)
                for z = 1:length(d{k})
                    d{k}{z} = rotatesim(d{k}{z},m.diagwidth);
                    if option.K < m.diagwidth
                        W = size(d{k}{z},1);
                        hW = ceil(W/2);
                        hK = ceil(option.K/2);
                        d{k}{z} = d{k}{z}(hW-hK:hW+hK,:);
                        m.diagwidth = option.K;
                    end
                end
            end
            m = set(m,'Data',d);
            m.view = 'h';
        elseif strcmpi(postoption.view,'TimeLag') || postoption.filt
            W = ceil(m.diagwidth/2);
            for k = 1:length(d)
                for z = 1:length(d{k})
                    if isinf(W)
                        dz = NaN(size(d{k}{z}));
                    else
                        dz = NaN(W,size(d{k}{z},2));
                    end
                    for l = 1:size(dz,1)
                        dz(l,l:end) = diag(d{k}{z},l-1)';
                    end
                    if option.K < m.diagwidth
                        dz = dz(1:ceil(option.K/2),:);
                        m.diagwidth = option.K;
                    end
                    d{k}{z}= dz;
                end
            end
            m = set(m,'Data',d);
            m.view = 'l';
        end
    end
    if ischar(postoption.simf)
        if strcmpi(postoption.simf,'Similarity')
            if not(isequal(m.similarity,NaN))
                option.dissim = 0;
            end
            postoption.simf = 'oneminus';
        end
        if isequal(m.similarity,0) && isstruct(option) ...
                && isfield(option,'dissim') && not(option.dissim)
            simf = str2func(postoption.simf);
            for k = 1:length(d)
                for z = 1:length(d{k})
                     d{k}{z} = simf(d{k}{z});
                end
            end
            m.similarity = postoption.simf;
            m = set(m,'Title','Similarity matrix','Data',d);
        elseif length(m.similarity) == 1 && isnan(m.similarity) ...
                && option.dissim
            m.similarity = 0;
        end
    end
    if postoption.filt
        fp = get(m,'FramePos');
        for k = 1:length(d)
            for z = 1:length(d{k})
                dz = filter(ones(postoption.filt,1),1,d{k}{z});
                d{k}{z} = dz(postoption.filt:end,1:end-postoption.filt+1);
                fp{k}{z} = [fp{k}{z}(1,1:end-postoption.filt+1);...
                            fp{k}{z}(1,postoption.filt:end)];
            end
        end
        m = set(m,'Data',d,'FramePos',fp);
    end
end


function S = rotatesim(d,K)
if length(d) == 1;
    S = d;
else
    K = min(K,size(d,1)*2);
    lK = floor(K/2);
    S = NaN(K,size(d,2),size(d,3));
    for k = 1:size(d,3)
        for j = -lK:lK
            S(lK+1+j,:,k) = [NaN(1,floor(abs(j)/2)) diag(d(:,:,k),j)' ...
                                                    NaN(1,ceil(abs(j)/2))];
        end
    end
end

function d = cosine(r,s)
d = 1-r'*s;
%nr = sqrt(r'*r);
%ns = sqrt(s'*s);
%if or(nr == 0, ns == 0);
%    d = 1;
%else
%    d = 1 - r'*s/nr/ns;
%end


function d = KL(x,y)
% Kullback-Leibler distance
if size(x,4)>1
    x(end+1:2*end,:,:,1) = x(:,:,:,2);
    x(:,:,:,2) = [];
end
if size(y,4)>1
    y(end+1:2*end,:,:,1) = y(:,:,:,2);
    y(:,:,:,2) = [];
end
m1 = mean(x);
m2 = mean(y);
S1 = cov(x);
S2 = cov(y);
d = (trace(S1/S2)+trace(S2/S1)+(m1-m2)'*inv(S1+S2)*(m1-m2))/2 - size(S1,1);
    

function s = exponential(d)
    s = exp(-d);
    
    
function s = oneminus(d)
    s = 1-d;