Daniel@0: function varargout = mirsimatrix(orig,varargin)
Daniel@0: %   m = mirsimatrix(x) computes the similarity matrix resulting from the 
Daniel@0: %       mutual comparison between each possible frame analysis in x.
Daniel@0: %           By default, x is the spectrum of the frame decomposition.
Daniel@0: %           But it can be any other frame analysis.
Daniel@0: %   Optional argument:
Daniel@0: %       mirsimatrix(...,'Distance',f) specifies the name of a distance
Daniel@0: %           function, from those proposed in the Statistics Toolbox
Daniel@0: %               (help pdist).
Daniel@0: %           default value: f = 'cosine'
Daniel@0: %       mirsimatrix(...,'Dissimilarity') return the dissimilarity matrix,
Daniel@0: %           which is the intermediary result before the computation of the 
Daniel@0: %           actual similarity matrix. It shows the distance between each
Daniel@0: %           possible frame analysis in x.
Daniel@0: %       mirsimatrix(...,'Similarity',f) indicates the function f specifying
Daniel@0: %           the way the distance values in the dissimilarity matrix are
Daniel@0: %           transformed into similarity values.
Daniel@0: %           Possible values:
Daniel@0: %               f = 'oneminus' (default value) 
Daniel@0: %                   corresponding to f(x) = 1-x
Daniel@0: %               f = 'exponential'
Daniel@0: %                   corresponding to f(x)= exp(-x)
Daniel@0: %       mirsimatrix(...,'Width',w) or more simply dissimatrix(...,w) 
Daniel@0: %           specifies the size of the diagonal bandwidth, in samples,
Daniel@0: %           outside which the dissimilarity will not be computed.
Daniel@0: %           if w = inf (default value), all the matrix will be computed.
Daniel@0: %       mirsimatrix(...,'Horizontal') rotates the matrix 45 degrees in order to
Daniel@0: %           make the first diagonal horizontal, and to restrict on the
Daniel@0: %           diagonal bandwidth only.
Daniel@0: %       mirsimatrix(...,'TimeLag') transforms the (non-rotated) matrix into
Daniel@0: %           a time-lag matrix, making the first diagonal horizontal as well
Daniel@0: %           (corresponding to the zero-lag line).
Daniel@0: %
Daniel@0: %       mirsimatrix(M,r) creates a mirsimatrix similarity matrix based on
Daniel@0: %           the Matlab square matrix M, of frame rate r (in Hz.)
Daniel@0: %               By default r = 20 Hz.
Daniel@0: %           mirsimatrix(M,r,'Dissimilarity') creates instead a mirsimatrix
Daniel@0: %               dissimilarity matrix.
Daniel@0: %
Daniel@0: %   Foote, J. & Cooper, M. (2003). Media Segmentation using Self-Similarity
Daniel@0: % Decomposition,. In Proc. SPIE Storage and Retrieval for Multimedia
Daniel@0: % Databases, Vol. 5021, pp. 167-75.
Daniel@0: 
Daniel@0: %       mirsimatrix(...,'Filter',10) filter along the diagonal of the matrix,
Daniel@0: %           using a uniform moving average filter of size f. The result is
Daniel@0: %           represented in a time-lag matrix.
Daniel@0: 
Daniel@0: 
Daniel@0:         distance.key = 'Distance';
Daniel@0:         distance.type = 'String';
Daniel@0:         distance.default = 'cosine';
Daniel@0:     option.distance = distance;
Daniel@0: 
Daniel@0:         simf.key = 'Similarity';
Daniel@0:         simf.type = 'String';
Daniel@0:         simf.default = 'oneminus';
Daniel@0:         simf.keydefault = 'Similarity';        
Daniel@0:         simf.when = 'After';
Daniel@0:     option.simf = simf;
Daniel@0: 
Daniel@0:         dissim.key = 'Dissimilarity';
Daniel@0:         dissim.type = 'Boolean';
Daniel@0:         dissim.default = 0;
Daniel@0:     option.dissim = dissim;
Daniel@0:     
Daniel@0:         K.key = 'Width';
Daniel@0:         K.type = 'Integer';
Daniel@0:         K.default = Inf;
Daniel@0:     option.K = K;
Daniel@0:     
Daniel@0:         filt.key = 'Filter';
Daniel@0:         filt.type = 'Integer';
Daniel@0:         filt.default = 0;
Daniel@0:         filt.when = 'After';
Daniel@0:     option.filt = filt;
Daniel@0: 
Daniel@0:         view.type = 'String';
Daniel@0:         view.default = 'Standard';
Daniel@0:         view.choice = {'Standard','Horizontal','TimeLag'};
Daniel@0:         view.when = 'After';
Daniel@0:     option.view = view;
Daniel@0:     
Daniel@0:         rate.type = 'Integer';
Daniel@0:         rate.position = 2;
Daniel@0:         rate.default = 20;
Daniel@0:     option.rate = rate;
Daniel@0:     
Daniel@0: specif.option = option;
Daniel@0: specif.nochunk = 1;
Daniel@0: varargout = mirfunction(@mirsimatrix,orig,varargin,nargout,specif,@init,@main);
Daniel@0: 
Daniel@0: 
Daniel@0: function [x type] = init(x,option)
Daniel@0: if isnumeric(x)
Daniel@0:     m.diagwidth = Inf;
Daniel@0:     m.view = 's';
Daniel@0:     m.similarity = NaN;
Daniel@0:     m.graph = {};
Daniel@0:     m.branch = {};
Daniel@0:     m = class(m,'mirsimatrix',mirdata);
Daniel@0:     m = set(m,'Title','Dissimilarity matrix');
Daniel@0:     fp = repmat(((1:size(x,1))-.5)/option.rate,[2,1]);
Daniel@0:     x = set(m,'Data',{x},'Pos',[],...
Daniel@0:               'FramePos',{{fp}},'Name',{inputname(1)});
Daniel@0: else
Daniel@0:     if not(isamir(x,'mirsimatrix'))
Daniel@0:         if (isamir(x,'miraudio'))
Daniel@0:             if isframed(x)
Daniel@0:                 x = mirspectrum(x);
Daniel@0:             else
Daniel@0:                 x = mirspectrum(x,'Frame',0.05,1);
Daniel@0:             end
Daniel@0:         end
Daniel@0:     end
Daniel@0: end
Daniel@0: type = 'mirsimatrix';
Daniel@0: 
Daniel@0: 
Daniel@0: function m = main(orig,option,postoption)
Daniel@0: if iscell(orig)
Daniel@0:     orig = orig{1};
Daniel@0: end
Daniel@0: if isa(orig,'mirsimatrix')
Daniel@0:     d = get(orig,'Data');
Daniel@0:     for k = 1:length(d)
Daniel@0:         if iscell(option) && isfield(option,'K') && option.K < orig.diagwidth
Daniel@0:             nl = size(d{k},1);
Daniel@0:             if strcmp(orig.view,'h')
Daniel@0:                 dl = (nl - option.K)/2;
Daniel@0:                 dk = d{k}(ceil(dl):nl-floor(dl),:);
Daniel@0:             else
Daniel@0:                 [spA,spd] = spdiags(d{k},-ceil(option.K/2):ceil(option.K/2));
Daniel@0:                 dk = full(spdiags(spA,spd,nl,size(d{k},2)));
Daniel@0:             end
Daniel@0:             d{k} = dk;
Daniel@0:             orig.diagwidth = option.K;
Daniel@0:         end
Daniel@0:     end
Daniel@0:     m = set(orig,'Data',d);
Daniel@0: else
Daniel@0:     v = get(orig,'Data');
Daniel@0:     d = cell(1,length(v));
Daniel@0:     for k = 1:length(v)
Daniel@0:         vk = v{k};
Daniel@0:         if mirwaitbar
Daniel@0:             handle = waitbar(0,'Computing dissimilarity matrix...');
Daniel@0:         else
Daniel@0:             handle = 0;
Daniel@0:         end
Daniel@0:         if 0 %iscell(vk)
Daniel@0:             try
Daniel@0:                 vk = cell2mat(vk);
Daniel@0:             end
Daniel@0:         end
Daniel@0:         if 0 %iscell(vk) %&& length(vk) > 1 %%% ATTENTION KL!!<<<<<<<<<<<<
Daniel@0:             l = length(vk);
Daniel@0:             dk = NaN(l,l);
Daniel@0:             hK = floor(option.K/2);
Daniel@0:             for i = 1:l
Daniel@0:                 if handle
Daniel@0:                     waitbar(i/l,handle);
Daniel@0:                 end
Daniel@0:                 for j = max(1,i-hK):min(l,i+hK)
Daniel@0:                     dk(i,j) = KL(vk{i},vk{j});
Daniel@0:                 end
Daniel@0:             end
Daniel@0:         else
Daniel@0:             if not(iscell(vk))
Daniel@0:                 vk = {vk};
Daniel@0:             end    
Daniel@0:             for z = 1:length(vk)
Daniel@0:                 vz = vk{z};
Daniel@0:                 ll = size(vz,1);
Daniel@0:                 l = size(vz,2);
Daniel@0:                 nc = size(vz,3);
Daniel@0:                 if ll==1 && nc>1
Daniel@0:                     vz = squeeze(vz)';
Daniel@0:                     ll = nc;
Daniel@0:                     nc = 1;
Daniel@0:                 end
Daniel@0:                 nd = size(vz,4);
Daniel@0:                 if not(isempty(postoption)) && ...
Daniel@0:                     strcmpi(postoption.view,'TimeLag')
Daniel@0:                     lK = ceil(option.K/2);
Daniel@0:                     if isinf(lK)
Daniel@0:                         lK = l;
Daniel@0:                     end
Daniel@0:                     dk{z} = NaN(lK,l,nc);
Daniel@0:                 else
Daniel@0:                     dk{z} = NaN(l,l,nc);
Daniel@0:                 end
Daniel@0:                 for g = 1:nc
Daniel@0:                     if nd == 1
Daniel@0:                         vv = vz;
Daniel@0:                     else
Daniel@0:                         vv = zeros(ll*nd,l);
Daniel@0:                         for h = 1:nd
Daniel@0:                             if iscell(vz)
Daniel@0:                                 for i = 1:ll
Daniel@0:                                     for j = 1:l
Daniel@0:                                         vj = vz{i,j,g,h};
Daniel@0:                                         if isempty(vj)
Daniel@0:                                             vv((h-1)*ll+i,j) = NaN;
Daniel@0:                                         else
Daniel@0:                                             vv((h-1)*ll+i,j) = vj;
Daniel@0:                                         end
Daniel@0:                                     end
Daniel@0:                                 end
Daniel@0:                             else
Daniel@0:                                 tic
Daniel@0:                                 vv((h-1)*ll+1:h*ll,:) = vz(:,:,g,h);
Daniel@0:                                 toc
Daniel@0:                             end
Daniel@0:                         end
Daniel@0:                     end
Daniel@0:                     if isinf(option.K) && not(strcmpi(postoption.view,'TimeLag'))
Daniel@0:                         try
Daniel@0:                             manually = 0;
Daniel@0:                             dk{z}(:,:,g) = squareform(pdist(vv',option.distance));
Daniel@0:                             if option.K < Inf
Daniel@0:                                 [spA,spd] = spdiags(dk{z},...
Daniel@0:                                     -ceil(option.K/2):ceil(option.K/2));
Daniel@0:                                 dk{z}(:,:,g) = full(spdiags(spA,spd,size(dk,1),size(dk{z},2)));
Daniel@0:                             end
Daniel@0:                         catch
Daniel@0:                             %err = lasterror;
Daniel@0:                             %warning(err.message)
Daniel@0:                             %disp('Statistics Toolbox does not seem to be
Daniel@0:                             %installed. Recompute the distance matrix manually.');
Daniel@0:                             manually = 1;
Daniel@0:                         end
Daniel@0:                     else
Daniel@0:                         manually = 1;
Daniel@0:                     end
Daniel@0:                     if manually
Daniel@0:                         disf = str2func(option.distance);
Daniel@0:                         if strcmpi(option.distance,'cosine')
Daniel@0:                             for i = 1:l
Daniel@0:                                 vv(:,i) = vv(:,i)/norm(vv(:,i));
Daniel@0:                             end
Daniel@0:                         end
Daniel@0:                         if not(isempty(postoption)) && ...
Daniel@0:                                 strcmpi(postoption.view,'TimeLag')
Daniel@0:                             lK = ceil(option.K/2);
Daniel@0:                             for i = 1:l
Daniel@0:                                 if mirwaitbar && (mod(i,100) == 1 || i == l)
Daniel@0:                                     waitbar(i/l,handle);
Daniel@0:                                 end
Daniel@0:                                 ij = min(i+lK-1,l);
Daniel@0:                                 dkij = disf(vv(:,i),vv(:,i:ij));
Daniel@0:                                 for j = 1:ij-i+1
Daniel@0:                                     dk{z}(j,i+j-1,g) = dkij(j);
Daniel@0:                                 end
Daniel@0:                             end
Daniel@0:                         else
Daniel@0:                             for i = 1:l
Daniel@0:                                 if mirwaitbar && (mod(i,100) == 1 || i == l)
Daniel@0:                                     waitbar(i/l,handle);
Daniel@0:                                 end
Daniel@0:                                 j = min(i+option.K-1,l);
Daniel@0:                                 dkij = disf(vv(:,i),vv(:,i:j));
Daniel@0:                                 dk{z}(i,i:j,g) = dkij;
Daniel@0:                                 dk{z}(i:j,i,g) = dkij';
Daniel@0:                             end
Daniel@0:                         end
Daniel@0:                     end
Daniel@0:                 end
Daniel@0:             end
Daniel@0:         end
Daniel@0:         d{k} = dk;
Daniel@0:         if handle
Daniel@0:             delete(handle)
Daniel@0:         end
Daniel@0:     end
Daniel@0:     m.diagwidth = option.K;
Daniel@0:     if not(isempty(postoption)) && strcmpi(postoption.view,'TimeLag')
Daniel@0:         m.view = 'l';
Daniel@0:     else
Daniel@0:         m.view = 's';
Daniel@0:     end
Daniel@0:     m.similarity = 0;
Daniel@0:     m.graph = {};
Daniel@0:     m.branch = {};
Daniel@0:     m = class(m,'mirsimatrix',mirdata(orig));
Daniel@0:     m = purgedata(m);
Daniel@0:     m = set(m,'Title','Dissimilarity matrix');
Daniel@0:     m = set(m,'Data',d,'Pos',[]);
Daniel@0: end
Daniel@0: if not(isempty(postoption))
Daniel@0:     if strcmpi(m.view,'s')
Daniel@0:         if strcmpi(postoption.view,'Horizontal')
Daniel@0:             for k = 1:length(d)
Daniel@0:                 for z = 1:length(d{k})
Daniel@0:                     d{k}{z} = rotatesim(d{k}{z},m.diagwidth);
Daniel@0:                     if option.K < m.diagwidth
Daniel@0:                         W = size(d{k}{z},1);
Daniel@0:                         hW = ceil(W/2);
Daniel@0:                         hK = ceil(option.K/2);
Daniel@0:                         d{k}{z} = d{k}{z}(hW-hK:hW+hK,:);
Daniel@0:                         m.diagwidth = option.K;
Daniel@0:                     end
Daniel@0:                 end
Daniel@0:             end
Daniel@0:             m = set(m,'Data',d);
Daniel@0:             m.view = 'h';
Daniel@0:         elseif strcmpi(postoption.view,'TimeLag') || postoption.filt
Daniel@0:             W = ceil(m.diagwidth/2);
Daniel@0:             for k = 1:length(d)
Daniel@0:                 for z = 1:length(d{k})
Daniel@0:                     if isinf(W)
Daniel@0:                         dz = NaN(size(d{k}{z}));
Daniel@0:                     else
Daniel@0:                         dz = NaN(W,size(d{k}{z},2));
Daniel@0:                     end
Daniel@0:                     for l = 1:size(dz,1)
Daniel@0:                         dz(l,l:end) = diag(d{k}{z},l-1)';
Daniel@0:                     end
Daniel@0:                     if option.K < m.diagwidth
Daniel@0:                         dz = dz(1:ceil(option.K/2),:);
Daniel@0:                         m.diagwidth = option.K;
Daniel@0:                     end
Daniel@0:                     d{k}{z}= dz;
Daniel@0:                 end
Daniel@0:             end
Daniel@0:             m = set(m,'Data',d);
Daniel@0:             m.view = 'l';
Daniel@0:         end
Daniel@0:     end
Daniel@0:     if ischar(postoption.simf)
Daniel@0:         if strcmpi(postoption.simf,'Similarity')
Daniel@0:             if not(isequal(m.similarity,NaN))
Daniel@0:                 option.dissim = 0;
Daniel@0:             end
Daniel@0:             postoption.simf = 'oneminus';
Daniel@0:         end
Daniel@0:         if isequal(m.similarity,0) && isstruct(option) ...
Daniel@0:                 && isfield(option,'dissim') && not(option.dissim)
Daniel@0:             simf = str2func(postoption.simf);
Daniel@0:             for k = 1:length(d)
Daniel@0:                 for z = 1:length(d{k})
Daniel@0:                      d{k}{z} = simf(d{k}{z});
Daniel@0:                 end
Daniel@0:             end
Daniel@0:             m.similarity = postoption.simf;
Daniel@0:             m = set(m,'Title','Similarity matrix','Data',d);
Daniel@0:         elseif length(m.similarity) == 1 && isnan(m.similarity) ...
Daniel@0:                 && option.dissim
Daniel@0:             m.similarity = 0;
Daniel@0:         end
Daniel@0:     end
Daniel@0:     if postoption.filt
Daniel@0:         fp = get(m,'FramePos');
Daniel@0:         for k = 1:length(d)
Daniel@0:             for z = 1:length(d{k})
Daniel@0:                 dz = filter(ones(postoption.filt,1),1,d{k}{z});
Daniel@0:                 d{k}{z} = dz(postoption.filt:end,1:end-postoption.filt+1);
Daniel@0:                 fp{k}{z} = [fp{k}{z}(1,1:end-postoption.filt+1);...
Daniel@0:                             fp{k}{z}(1,postoption.filt:end)];
Daniel@0:             end
Daniel@0:         end
Daniel@0:         m = set(m,'Data',d,'FramePos',fp);
Daniel@0:     end
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: function S = rotatesim(d,K)
Daniel@0: if length(d) == 1;
Daniel@0:     S = d;
Daniel@0: else
Daniel@0:     K = min(K,size(d,1)*2);
Daniel@0:     lK = floor(K/2);
Daniel@0:     S = NaN(K,size(d,2),size(d,3));
Daniel@0:     for k = 1:size(d,3)
Daniel@0:         for j = -lK:lK
Daniel@0:             S(lK+1+j,:,k) = [NaN(1,floor(abs(j)/2)) diag(d(:,:,k),j)' ...
Daniel@0:                                                     NaN(1,ceil(abs(j)/2))];
Daniel@0:         end
Daniel@0:     end
Daniel@0: end
Daniel@0: 
Daniel@0: function d = cosine(r,s)
Daniel@0: d = 1-r'*s;
Daniel@0: %nr = sqrt(r'*r);
Daniel@0: %ns = sqrt(s'*s);
Daniel@0: %if or(nr == 0, ns == 0);
Daniel@0: %    d = 1;
Daniel@0: %else
Daniel@0: %    d = 1 - r'*s/nr/ns;
Daniel@0: %end
Daniel@0: 
Daniel@0: 
Daniel@0: function d = KL(x,y)
Daniel@0: % Kullback-Leibler distance
Daniel@0: if size(x,4)>1
Daniel@0:     x(end+1:2*end,:,:,1) = x(:,:,:,2);
Daniel@0:     x(:,:,:,2) = [];
Daniel@0: end
Daniel@0: if size(y,4)>1
Daniel@0:     y(end+1:2*end,:,:,1) = y(:,:,:,2);
Daniel@0:     y(:,:,:,2) = [];
Daniel@0: end
Daniel@0: m1 = mean(x);
Daniel@0: m2 = mean(y);
Daniel@0: S1 = cov(x);
Daniel@0: S2 = cov(y);
Daniel@0: d = (trace(S1/S2)+trace(S2/S1)+(m1-m2)'*inv(S1+S2)*(m1-m2))/2 - size(S1,1);
Daniel@0:     
Daniel@0: 
Daniel@0: function s = exponential(d)
Daniel@0:     s = exp(-d);
Daniel@0:     
Daniel@0:     
Daniel@0: function s = oneminus(d)
Daniel@0:     s = 1-d;