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

Added dictionary decorrelation for multiple dictionaries

--- a/DL/two-step DL/dico_decorr_symetric.m	Wed Mar 14 14:42:52 2012 +0000
+++ b/DL/two-step DL/dico_decorr_symetric.m	Tue Mar 20 12:25:50 2012 +0000
@@ -1,15 +1,17 @@
-function dico = dico_decorr_symetric(dico, mu, amp)
+function dico = dico_decorr_symetric(dico, mu)
     %DICO_DECORR decorrelate a dictionary
     %   Parameters:
-    %   dico: the dictionary
+    %   dico: the dictionary, either a matrix or a cell array of matrices.
     %   mu: the coherence threshold
-    %   amp: the amplitude coefficients, only used to decide which atom to
-    %   project
     %
     %   Result:
-    %   dico: a dictionary close to the input one with coherence mu.
+    %   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-6; % define tolerance for normalisation term alpha
+    eps = 1e-3; % define tolerance for normalisation term alpha
     
     % convert mu to the to the mean direction
     theta = acos(mu)/2;
@@ -23,39 +25,115 @@
     %         rank = randperm(length(dico));
     %     end
     
-    % 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 + 0.01
-        % find pairs of high correlation atoms
-        colors = dico_color(dico, mu);
+    % 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
         
-        % 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               
+        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(:,index(1))+dico(:,index(2));
+                    v1 = dico{i}(:,m)+dico{j}(:,n);
                     v1 = v1/norm(v1);
-                    v2 = dico(:,index(1))-dico(:,index(2));
+                    v2 = dico{i}(:,m)-dico{j}(:,n);
                     v2 = v2/norm(v2);
                     
-                    dico(:,index(1)) = ctheta*v1+stheta*v2;
-                    dico(:,index(2)) = ctheta*v1-stheta*v2;
+                    dico{i}(:,m) = ctheta*v1+stheta*v2;
+                    dico{j}(:,n) = ctheta*v1-stheta*v2;
                 else
-                    v1 = dico(:,index(1))-dico(:,index(2));
+                    v1 = dico{i}(:,m)-dico{j}(:,n);
                     v1 = v1/norm(v1);
-                    v2 = dico(:,index(1))+dico(:,index(2));
+                    v2 = dico{i}(:,m)+dico{j}(:,n);
                     v2 = v2/norm(v2);
                     
-                    dico(:,index(1)) = ctheta*v1+stheta*v2;
-                    dico(:,index(2)) = -ctheta*v1+stheta*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
-        niter = niter+1;
     end
 end
-