fixed bugs that were causing inaccurate results when running in patches

- spatial overlap of patches is accounted for
- fixed a bug when reading Y from memory to compute YrA

Author: epnev

run_CNMF_patches.m

@@ -137,6 +137,9 @@
 end
 
 cnt = 0;
+B = sparse(prod(sizY(1:end-1)),length(patches));
+MASK = zeros(sizY(1:end-1));
+F = zeros(length(patches),sizY(end));
 for i = 1:length(patches)
     for k = 1:K
         if k <= size(RESULTS(i).A,2)
@@ -152,6 +155,9 @@
         end
     end
     if length(sizY) == 3
+        b_temp = sparse(sizY(1),sizY(2));
+        b_temp(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4)) = reshape(RESULTS(i).b,patches{i}(2)-patches{i}(1)+1,patches{i}(4)-patches{i}(3)+1);
+        MASK(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4)) = MASK(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4)) + 1;
         P.sn(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4)) = reshape(RESULTS(i).P.sn,patches{i}(2)-patches{i}(1)+1,patches{i}(4)-patches{i}(3)+1);
         P.active_pixels(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4)) = P.active_pixels(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4)) + ...
             reshape(RESULTS(i).P.active_pixels,patches{i}(2)-patches{i}(1)+1,patches{i}(4)-patches{i}(3)+1);
@@ -167,8 +173,12 @@
     P.c1 = [P.c1;RESULTS(i).P.c1];
     P.gn = [P.gn;RESULTS(i).P.gn];
     P.neuron_sn = [P.neuron_sn;RESULTS(i).P.neuron_sn];
+    B(:,i) = b_temp(:);
+    F(i,:) = RESULTS(i).f;
 end
 A(:,cnt+1:end) = [];
+A = spdiags(1./MASK(:),0,prod(sizY(1:end-1)),prod(sizY(1:end-1)))*A;
+B = spdiags(1./MASK(:),0,prod(sizY(1:end-1)),prod(sizY(1:end-1)))*B;
 C = cell2mat({RESULTS(:).C}');
 S = cell2mat({RESULTS(:).S}');
 ff = find(sum(A,1)==0);
@@ -207,9 +217,11 @@
 end
 fprintf(' done. \n');
 %% classify components
-ff = classify_components(Am,Pm,options);
-A = Am(:,ff);
-C = Cm(ff,:);
+% ff = classify_components(Am,Pm,options);
+% A = Am(:,ff);
+% C = Cm(ff,:);
+A = Am;
+C = Cm;
 
 %% update spatial components
 fprintf('Updating spatial components...');
@@ -226,15 +238,21 @@
 %     cnt = cnt + fm(i)-fp(i)+1;
 % end
 
-bsum = zeros(length(patches),1);
-for i = 1:length(patches)
-    bsum(i) = sum(RESULTS(i).b);
-end
-bsum = bsum/sum(bsum);
-f_p = cell2mat({RESULTS(:).f}');
-fin = mean(spdiags(bsum,0,length(patches),length(patches))*f_p);
-fin = medfilt1(fin,11);
+% bsum = zeros(length(patches),1);
+% for i = 1:length(patches)
+%     bsum(i) = sum(RESULTS(i).b);
+% end
+% bsum = bsum/sum(bsum);
+% f_p = cell2mat({RESULTS(:).f}');
+% fin = mean(spdiags(bsum,0,length(patches),length(patches))*f_p);
+% fin = medfilt1(fin,11);
 
+fin = mean(F);
+for iter = 1:10
+    bin = max(B*(F*fin')/norm(fin)^2,0);
+    fin = max((bin'*B)*F/norm(bin)^2,0);
+end
+%%
 % PROCESS PATCHES 
 options.d1 = sizY(1);
 options.d2 = sizY(2);

update_temporal_components.m

@@ -142,18 +142,19 @@
     mc = min(d,15);  % number of constraints to be considered
     LD = 10*ones(mc,K);
 else
+    step = 5e3;
     nA = sum(A.^2);
-    AA = A'*A/spdiags(nA(:),0,length(nA),length(nA));
+    AA = (A'*A)/spdiags(nA(:),0,length(nA),length(nA));
     if memmaped
         YA = zeros(T,length(nA));
         tic;
-        for i = 1:5e3:d
-            YA = YA + double(Y.Yr(i:min(i+1e4-1,d),:))'*A(i:min(i+1e4-1,d),:);
+        for i = 1:step:d
+            YA = YA + double(Y.Yr(i:min(i+step-1,d),:))'*A(i:min(i+step-1,d),:);
         end
         toc
         YA = YA/spdiags(nA(:),0,length(nA),length(nA));
     else
-        YA = Y'*A/spdiags(nA(:),0,length(nA),length(nA));
+        YA = (Y'*A)/spdiags(nA(:),0,length(nA),length(nA));
     end
     YrA = (YA - Cin'*AA);
     if strcmpi(method,'constrained_foopsi') || strcmpi(method,'MCEM_foopsi')