view DL/two-step DL/dico_decorr_symetric.m @ 200:69ce11724b1f luisf_dev

Added dictionary decorrelation for multiple dictionaries

function dico = dico_decorr_symetric(dico, mu)
    %DICO_DECORR decorrelate a dictionary
    %   Parameters:
    %   dico: the dictionary, either a matrix or a cell array of matrices.
    %   mu: the coherence threshold
    %
    %   Result:
    %   dico: if the input dico was a matrix, then a matrix close to the 
    %   input one with coherence mu.
    %   If the input was a cell array, a cell array of the same size
    %   containing matrices such that the coherence between different cells
    %   is lower than mu.
    
    eps = 1e-3; % define tolerance for normalisation term alpha
    
    % convert mu to the to the mean direction
    theta = acos(mu)/2;
    ctheta = cos(theta);
    stheta = sin(theta);
    
    % compute atom weights
    %     if nargin > 2
    %         rank = sum(amp.*amp, 2);
    %     else
    %         rank = randperm(length(dico));
    %     end
    
    % if only one dictionary is provided, then decorrelate it
    if ~iscell(dico)
        % several decorrelation iterations might be needed to reach global
        % coherence mu. niter can be adjusted to needs.
        niter = 1;
        while max(max(abs(dico'*dico -eye(length(dico))))) > mu + eps
            % find pairs of high correlation atoms
            colors = dico_color(dico, mu);
            
            % iterate on all pairs
            nbColors = max(colors);
            for c = 1:nbColors
                index = find(colors==c);
                if numel(index) == 2
                    if dico(:,index(1))'*dico(:,index(2)) > 0
                        %build the basis vectors
                        v1 = dico(:,index(1))+dico(:,index(2));
                        v1 = v1/norm(v1);
                        v2 = dico(:,index(1))-dico(:,index(2));
                        v2 = v2/norm(v2);
                        
                        dico(:,index(1)) = ctheta*v1+stheta*v2;
                        dico(:,index(2)) = ctheta*v1-stheta*v2;
                    else
                        v1 = dico(:,index(1))-dico(:,index(2));
                        v1 = v1/norm(v1);
                        v2 = dico(:,index(1))+dico(:,index(2));
                        v2 = v2/norm(v2);
                        
                        dico(:,index(1)) = ctheta*v1+stheta*v2;
                        dico(:,index(2)) = -ctheta*v1+stheta*v2;
                    end
                end
            end
            niter = niter+1;
        end
        %if a cell array of dictionaries is provided, decorrelate among
        %different dictionaries only
    else
        niter = 1;
        numDicos = length(dico);
        G = cell(numDicos);
        maxCorr = 0;
        for i = 1:numDicos
            for j = i+1:numDicos
                G{i,j} = dico{i}'*dico{j};
                maxCorr = max(maxCorr,max(max(abs(G{i,j}))));
            end
        end
        
        while maxCorr > mu + eps
            % find pairs of high correlation atoms
            [colors nbColors] = dico_color_separate(dico, mu);

            % iterate on all pairs
            for c = 1:nbColors
                for tmpI = 1:numDicos
                    index = find(colors{tmpI}==c);
                    if ~isempty(index)
                        i = tmpI;
                        m = index;
                        break;
                    end
                end
                for tmpJ = i+1:numDicos
                    index = find(colors{tmpJ}==c);
                    if ~isempty(index)
                        j = tmpJ;
                        n = index;
                        break;
                    end
                end
                
                if dico{i}(:,m)'*dico{j}(:,n) > 0
                    %build the basis vectors
                    v1 = dico{i}(:,m)+dico{j}(:,n);
                    v1 = v1/norm(v1);
                    v2 = dico{i}(:,m)-dico{j}(:,n);
                    v2 = v2/norm(v2);
                    
                    dico{i}(:,m) = ctheta*v1+stheta*v2;
                    dico{j}(:,n) = ctheta*v1-stheta*v2;
                else
                    v1 = dico{i}(:,m)-dico{j}(:,n);
                    v1 = v1/norm(v1);
                    v2 = dico{i}(:,m)+dico{j}(:,n);
                    v2 = v2/norm(v2);
                    
                    dico{i}(:,m) = ctheta*v1+stheta*v2;
                    dico{j}(:,n) = -ctheta*v1+stheta*v2;
                end
            end
            niter = niter+1;
            
            % Remove noegative components and renormalize
            for i = 1:length(dico)
                dico{i} = max(dico{i},0);
                for m = 1:size(dico{i},2)
                    dico{i}(:,m) = dico{i}(:,m)/norm(dico{i}(:,m));
                end
            end
            
            maxCorr = 0;
            for i = 1:numDicos
                for j = i+1:numDicos
                    G{i,j} = dico{i}'*dico{j};
                    maxCorr = max(maxCorr,max(max(abs(G{i,j}))));
                end
            end
        end
    end
end