fix blurring of residual when adding component in OnACID

caiman/source_extraction/cnmf/online_cnmf.py

@@ -797,7 +797,7 @@ def update_num_components(t, sv, Ab, Cf, Yres_buf, Y_buf, rho_buf,
                 groups = update_order(Ab, Ain, groups)[0]
             else:
                 groups = update_order(Ab_dense[indeces], ain, groups)[0]
-                Ab_dense = np.hstack((Ab_dense,Ain))
+                Ab_dense = np.hstack((Ab_dense, Ain))
             # faster version of scipy.sparse.hstack
             csc_append(Ab, Ain_csc)
             ind_A.append(Ab.indices[Ab.indptr[M]:Ab.indptr[M + 1]])
@@ -821,34 +821,25 @@ def update_num_components(t, sv, Ab, Cf, Yres_buf, Y_buf, rho_buf,
             M = M + 1
 
             Yres_buf[:, indeces] -= np.outer(cin, ain)
-            # vb = imblur(np.reshape(Ain, dims, order='F'), sig=gSig,
-            #             siz=gSiz, nDimBlur=2).ravel()
-            # restrict blurring to region where component is located
-#            vb = np.reshape(Ain, dims, order='F')
-            slices = tuple(slice(max(0, ijs[0] - 2*sg), min(d, ijs[1] + 2*sg))
-                           for ijs, sg, d in zip(ijSig, gSiz//2, dims))  # is 2 enough?
-
-            slice_within = tuple(slice(ijs[0] - sl.start, ijs[1] - sl.start)
-                           for ijs, sl in zip(ijSig, slices))
 
+            # restrict blurring to region where component is located
+            # update bigger region than neural patch to avoid boundary effects
+            slices_update = tuple(slice(max(0, ijs[0] - sg // 2), min(d, ijs[1] + sg // 2))
+                                  for ijs, sg, d in zip(ijSig, gSiz, dims))
+            # filter even bigger region to avoid boundary effects
+            slices_filter = tuple(slice(max(0, ijs[0] - sg), min(d, ijs[1] + sg))
+                                  for ijs, sg, d in zip(ijSig, gSiz, dims))
 
             ind_vb = np.ravel_multi_index(
-               np.ix_(*[np.arange(ij[0], ij[1])
-                      for ij in ijSig]), dims, order='C').ravel()
-
-
-
-
-            vb_buf = [imblur(np.maximum(0,vb.reshape(dims,order='F')[slices][slice_within]), sig=gSig, siz=gSiz, nDimBlur=len(dims)) for vb in Yres_buf]
-
-            vb_buf2 = np.stack([vb.ravel() for vb in vb_buf])
-
-#            ind_vb = np.ravel_multi_index(
-#                    np.ix_(*[np.arange(s.start, s.stop)
-#                           for s in slices_small]), dims).ravel()
-
-            rho_buf[:, ind_vb] = vb_buf2**2
-
+                np.ix_(*[np.arange(sl.start, sl.stop)
+                         for sl in slices_update]), dims, order='C').ravel()
+
+            rho_buf[:, ind_vb] = np.stack([imblur(
+                vb.reshape(dims, order='F')[slices_filter], sig=gSig, siz=gSiz,
+                nDimBlur=len(dims))[tuple([slice(
+                    slices_update[i].start - slices_filter[i].start,
+                    slices_update[i].stop - slices_filter[i].start)
+                    for i in range(len(dims))])].ravel() for vb in Yres_buf])**2
 
             sv[ind_vb] = np.sum(rho_buf[:, ind_vb], 0)
 #            sv = np.sum([imblur(vb.reshape(dims,order='F'), sig=gSig, siz=gSiz, nDimBlur=len(dims))**2 for vb in Yres_buf], 0).reshape(-1)