allow seeding with binary masks and arbitrary background comps

Author: epnev

update_spatial_components.m

@@ -60,15 +60,15 @@
 end
 
 if nargin < 5 || isempty(P); P = preprocess_data(Y,1); end  % etsimate noise values if not present
-if nargin < 4 || isempty(A_); 
+if nargin < 4 || isempty(A_)
     IND = ones(d,size(C,1)); 
 else
     if islogical(A_)     % check if search locations have been provided, otherwise estimate them        
         IND = A_(:,1:K);
         if isempty(C)    
             INDav = double(IND)/diag(sum(double(IND)));          
             px = (sum(IND,2)>0);
-            f = mean(Y(~px,:));
+            [b,f] = fast_nmf(double(Y(~px,:)),[],options.nb,50);
             b = max(Y*f',0)/norm(f)^2;
             C = max(INDav'*Y - (INDav'*b)*f,0);
         end
@@ -116,7 +116,7 @@
 end
 f_ds = Cf_ds(end-size(f,1)+1:end,:);
 
-if strcmpi(options.spatial_method,'constrained');
+if strcmpi(options.spatial_method,'constrained')
     if spatial_parallel         % solve BPDN problem for each pixel
         Nthr = max(20*maxNumCompThreads,round(d*T/2^24));
         Nthr = min(Nthr,round(d/1e3));
@@ -177,7 +177,7 @@
 end
 
 A(isnan(A))=0;
-A = sparse(A);
+A = sparse(double(A));
 A = threshold_components(A,options);  % post-processing of components
 
 fprintf('Updated spatial components \n');

update_temporal_components.m

@@ -112,9 +112,10 @@
 if isempty(fin) || nargin < 5   % temporal background missing
     bk_pix = (sum(A,2)==0);     % pixels with no active neurons
     if isempty(b) || nargin < 3
-        fin = mean(Y(bk_pix,:));
-        fin = fin/norm(fin);
-        b = max(Y*fin',0);
+        [b,fin] = fast_nmf(double(Y(bk_pix,:)),[],options.nb,50);
+%         fin = mean(Y(bk_pix,:));
+%         fin = fin/norm(fin);
+         b = max(Y*fin',0);
     else
         fin = max(b(bk_pix,:)'*Y(bk_pix,:),0)/(b(bk_pix,:)'*b(bk_pix,:));
     end
@@ -124,9 +125,10 @@
 AY = mm_fun([A,double(b)],Y);
 bY = AY(size(A,2)+1:end,:);
 AY = AY(1:size(A,2),:);
+AA = sparse(double(A))'*sparse(double(A));
 
 if isempty(Cin) || nargin < 4    % estimate temporal components if missing    
-    Cin = max((A'*A)\double(AY - (A'*b)*fin),0);  
+    Cin = max(AA\double(AY - (A'*b)*fin),0);  
     ITER = max(ITER,3);
 end
 

utilities/fast_nmf.m

@@ -10,6 +10,9 @@
 H = H';
 
 for iter = 1:max_iter
+    if isempty(B)
+        B = 1;
+    end
     W = max((A*(B*H'))/(H*H'),0);
     H = max((W'*W)\((W'*A)*B),0);
     %disp(iter)
@@ -18,10 +21,14 @@
 function [W,H] = nnsvd(A,B,k)
 
     mf = @(x,y) mat_vec(A,B,x,y);    
-    [U,S,V] = svds(mf,[size(A,1),size(B,2)],k);
+    sz = [size(A,1), size(B,2)];
+    if isempty(B)
+        sz(2) = size(A,2);
+    end
+    [U,S,V] = svds(mf,sz,k);
 
-    W = zeros(size(A,1),k);
-    H = zeros(size(B,2),k);
+    W = zeros(sz(1),k);
+    H = zeros(sz(2),k);
 
     for j = 1:k
         x = U(:,j);
@@ -46,16 +53,29 @@
         W(:,j) = sqrt(S(j,j)*sigma)*u;
         H(:,j) = sqrt(S(j,j)*sigma)*v;
     end
-    avg = mean(A,1)*mean(B,2);
+    if isempty(B)
+        avg = mean(A(:));
+    else
+        avg = mean(A,1)*mean(B,2);
+    end
     W(W==0) = avg;
     H(H==0) = avg;
 
     function y = mat_vec(A,B,x,method)
 
         if strcmp(method,'notransp')
-            y = A*(B*x);
+            if isempty(B)
+                y = A*x;
+            else
+                y = A*(B*x);
+            end
         elseif strcmp(method,'transp')
-            y = B'*(A'*x);
+            if isempty(B)
+                y = A'*x;
+            else
+                y = B'*(A'*x);
+                
+            end
         end
     end
 end