Merge 7d2d964 into f57b9bf

data/test_data/detection/detect_output_1p1c0.npy.dvc

@@ -0,0 +1,7 @@
+md5: 0cb8789ec5068a987cb27d8cff2b1515
+outs:
+- md5: 58d29362cf4ab7ef3450dbe70f68385b
+  path: detect_output_1p1c0.npy
+  cache: true
+  metric: false
+  persist: false

data/test_data/detection/detect_output_2p2c0.npy.dvc

@@ -0,0 +1,7 @@
+md5: 29a50cca9549f174e60e6b7e7d9f5407
+outs:
+- md5: d6f0f1bb1f9abd8191543a4e74246c62
+  path: detect_output_2p2c0.npy
+  cache: true
+  metric: false
+  persist: false

data/test_data/detection/detect_output_2p2c1.npy.dvc

@@ -0,0 +1,7 @@
+md5: 203801255c0a2f6eebf0307e48aaf9cc
+outs:
+- md5: 81416ecadd0c0c29191a6f172355cd4d
+  path: detect_output_2p2c1.npy
+  cache: true
+  metric: false
+  persist: false

suite2p/detection/__init__.py

@@ -0,0 +1 @@
+from .detect import roi_stats, main_detect

suite2p/extraction/chan2detect.py → suite2p/detection/chan2detect.py

@@ -43,7 +43,6 @@ def correct_bleedthrough(Ly, Lx, nblks, mimg, mimg2):
     return mimg2
 
 def detect(ops, stat):
-    #ops2 = ops.copy()
     mimg = ops['meanImg'].copy()
     mimg2 = ops['meanImg_chan2'].copy()
 

suite2p/detection/detect.py

@@ -0,0 +1,114 @@
+import time
+import numpy as np
+from pathlib import Path
+from . import sourcery, sparsedetect, masks, chan2detect, utils
+
+
+def main_detect(ops, stat=None):
+    stat = select_rois(ops, stat)
+    # extract fluorescence and neuropil
+    cell_pix, cell_masks, neuropil_masks = make_masks(ops, stat)
+    ic = np.ones(len(stat), np.bool)
+    # if second channel, detect bright cells in second channel
+    if 'meanImg_chan2' in ops:
+        if 'chan2_thres' not in ops:
+            ops['chan2_thres'] = 0.65
+        ops, redcell = chan2detect.detect(ops, stat)
+        np.save(Path(ops['save_path']).joinpath('redcell.npy'), redcell[ic])
+    return cell_pix, cell_masks, neuropil_masks, stat, ops
+
+
+def select_rois(ops, stat=None):
+    t0 = time.time()
+    if stat is None:
+        if ops['sparse_mode']:
+            ops, stat = sparsedetect.sparsery(ops)
+        else:
+            ops, stat = sourcery.sourcery(ops)
+        print('Found %d ROIs, %0.2f sec' % (len(stat), time.time() - t0))
+    stat = roi_stats(ops, stat)
+
+    stat = masks.get_overlaps(stat, ops)
+    stat, ix = masks.remove_overlappers(stat, ops, ops['Ly'], ops['Lx'])
+    print('After removing overlaps, %d ROIs remain' % (len(stat)))
+    return stat
+
+
+def make_masks(ops, stat):
+    t0=time.time()
+    cell_pix, cell_masks = masks.create_cell_masks(stat, ops['Ly'], ops['Lx'], ops['allow_overlap'])
+    neuropil_masks = masks.create_neuropil_masks(ops, stat, cell_pix)
+    Ly=ops['Ly']
+    Lx=ops['Lx']
+    neuropil_masks = np.reshape(neuropil_masks, (-1,Ly*Lx))
+    print('Masks made in %0.2f sec.'%(time.time()-t0))
+    return cell_pix, cell_masks, neuropil_masks
+
+
+def roi_stats(ops, stat):
+    """ computes statistics of ROIs
+
+    Parameters
+    ----------
+    ops : dictionary
+        'aspect', 'diameter'
+
+    stat : dictionary
+        'ypix', 'xpix', 'lam'
+
+    Returns
+    -------
+    stat : dictionary
+        adds 'npix', 'npix_norm', 'med', 'footprint', 'compact', 'radius', 'aspect_ratio'
+
+    """
+    if 'aspect' in ops:
+        d0 = np.array([int(ops['aspect']*10), 10])
+    else:
+        d0 = ops['diameter']
+        if isinstance(d0, int):
+            d0 = [d0,d0]
+
+    rs = masks.circle_mask(np.array([30, 30]))
+    rsort = np.sort(rs.flatten())
+
+    ncells = len(stat)
+    mrs = np.zeros((ncells,))
+    for k in range(0,ncells):
+        stat0 = stat[k]
+        ypix = stat0['ypix']
+        xpix = stat0['xpix']
+        lam = stat0['lam']
+        # compute footprint of ROI
+        y0 = np.median(ypix)
+        x0 = np.median(xpix)
+
+        # compute compactness of ROI
+        r2 = ((ypix-y0))**2 + ((xpix-x0))**2
+        r2 = r2**.5
+        stat0['mrs']  = np.mean(r2)
+        mrs[k] = stat0['mrs']
+        stat0['mrs0'] = np.mean(rsort[:r2.size])
+        stat0['compact'] = stat0['mrs'] / (1e-10+stat0['mrs0'])
+        stat0['med']  = [np.median(stat0['ypix']), np.median(stat0['xpix'])]
+        stat0['npix'] = xpix.size
+
+        if 'footprint' not in stat0:
+            stat0['footprint'] = 0
+        if 'med' not in stat:
+            stat0['med'] = [np.median(stat0['ypix']), np.median(stat0['xpix'])]
+        if 'radius' not in stat0:
+            radius = utils.fitMVGaus(ypix / d0[0], xpix / d0[1], lam, 2)[2]
+            stat0['radius'] = radius[0] * d0.mean()
+            stat0['aspect_ratio'] = 2 * radius[0]/(.01 + radius[0] + radius[1])
+
+    npix = np.array([stat[n]['npix'] for n in range(len(stat))]).astype('float32')
+    npix /= np.mean(npix[:100])
+
+    mmrs = np.nanmedian(mrs[:100])
+    for n in range(len(stat)):
+        stat[n]['mrs'] = stat[n]['mrs'] / (1e-10+mmrs)
+        stat[n]['npix_norm'] = npix[n]
+    stat = np.array(stat)
+
+    return stat

suite2p/extraction/masks.py → suite2p/detection/masks.py

@@ -1,6 +1,6 @@
 import numpy as np
 
-from .sparsedetect import extendROI
+from suite2p.detection.sparsedetect import extendROI
 
 
 def get_overlaps(stat, ops):

suite2p/extraction/sourcery.py → suite2p/detection/sourcery.py

@@ -22,19 +22,15 @@ def getSVDdata(ops):
     sdmov = utils.get_sdmov(mov, ops)
     mov /= sdmov
     mov = np.reshape(mov, (-1,Lyc*Lxc))
-    if 1:
-        # compute covariance of binned frames
-        cov = mov @ mov.transpose() / mov.shape[1]
-        cov = cov.astype('float32')
+    # compute covariance of binned frames
+    cov = mov @ mov.transpose() / mov.shape[1]
+    cov = cov.astype('float32')
 
-        nsvd_for_roi = min(ops['nbinned'], int(cov.shape[0]/2))
-        u, s, v = np.linalg.svd(cov)
+    nsvd_for_roi = min(ops['nbinned'], int(cov.shape[0]/2))
+    u, s, v = np.linalg.svd(cov)
 
-        u = u[:, :nsvd_for_roi]
-        U = u.transpose() @ mov
-    else:
-        U = mov
-        u = []
+    u = u[:, :nsvd_for_roi]
+    U = u.transpose() @ mov
     U = np.reshape(U, (-1,Lyc,Lxc))
     U = np.transpose(U, (1, 2, 0)).copy()
     return ops, U, sdmov, u

suite2p/extraction/utils.py → suite2p/detection/utils.py

@@ -1,7 +1,7 @@
+import numpy as np
 import math
 import time
 
-import numpy as np
 from scipy.ndimage import gaussian_filter
 
 
@@ -50,7 +50,6 @@ def bin_movie(ops):
     nimgbatch = bin_size * (nimgbatch // bin_size)
     nbytesread = np.int64(Ly*Lx*nimgbatch*2)
     mov = np.zeros((ops['nbinned'], Lyc, Lxc), np.float32)
-    max_proj = np.zeros((Lyc, Lxc), np.float32)
     ix = 0
     idata = 0
     # load and bin data
@@ -83,8 +82,6 @@ def bin_movie(ops):
     max_proj = mov.max(axis=0)
     print('Binned movie [%d,%d,%d], %0.2f sec.'%(mov.shape[0], mov.shape[1], mov.shape[2], time.time()-t0))
 
-    #nimgbatch = min(mov.shape[0] , max(int(500/nt0), int(240./nt0 * ops['fs'])))
-
     # data is high-pass filtered
     if ops['high_pass']<10:
         # slow high-pass

suite2p/extraction/__init__.py

@@ -1,3 +1,2 @@
 from .dcnv import preprocess, oasis
-from .extract import extract_traces, detect_and_extract, enhanced_mean_image, roi_stats
-from .masks import create_cell_masks, create_neuropil_masks, get_overlaps
+from .extract import extract_traces, extract, enhanced_mean_image