diff --git a/caiman/source_extraction/cnmf/online_cnmf.py b/caiman/source_extraction/cnmf/online_cnmf.py index d6a6275c1..98fef5459 100644 --- a/caiman/source_extraction/cnmf/online_cnmf.py +++ b/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)