Merge 7888a34 into 3d3d020

.dvc/.gitignore

@@ -8,3 +8,42 @@
 /cache
 /tmp/*
 
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp
+/tmp

.dvc/tmp/gdrive-user-credentials.json

@@ -1 +1 @@
-{"access_token": "ya29.a0AfH6SMAKnyeAb_neSiUApH-xr8LiYI-OxbU3zjZVdIDlm-PpSpZv58nPHw86980-0uwoYmkMEvAwCxbqFCz1QHwhu6Z9uTdkQNrSnJwfSpYGVQHO-TtFLLHXESEpZppkBVRhDrQgyZxezsidomVTYI8rYJyVu4dfqMU", "client_id": "710796635688-iivsgbgsb6uv1fap6635dhvuei09o66c.apps.googleusercontent.com", "client_secret": "a1Fz59uTpVNeG_VGuSKDLJXv", "refresh_token": "1//0fpJK6wboLxL7CgYIARAAGA8SNwF-L9IrKmhg8LLnJuWRCiN_YIE88ImD9Y0yV6iCIw220vP3gD4nISe-K5aq8l_r0m1qcUW1D_w", "token_expiry": "2020-06-17T02:04:27Z", "token_uri": "https://oauth2.googleapis.com/token", "user_agent": null, "revoke_uri": "https://oauth2.googleapis.com/revoke", "id_token": null, "id_token_jwt": null, "token_response": {"access_token": "ya29.a0AfH6SMAKnyeAb_neSiUApH-xr8LiYI-OxbU3zjZVdIDlm-PpSpZv58nPHw86980-0uwoYmkMEvAwCxbqFCz1QHwhu6Z9uTdkQNrSnJwfSpYGVQHO-TtFLLHXESEpZppkBVRhDrQgyZxezsidomVTYI8rYJyVu4dfqMU", "expires_in": 3599, "refresh_token": "1//0fpJK6wboLxL7CgYIARAAGA8SNwF-L9IrKmhg8LLnJuWRCiN_YIE88ImD9Y0yV6iCIw220vP3gD4nISe-K5aq8l_r0m1qcUW1D_w", "scope": "https://www.googleapis.com/auth/drive https://www.googleapis.com/auth/drive.appdata", "token_type": "Bearer"}, "scopes": ["https://www.googleapis.com/auth/drive", "https://www.googleapis.com/auth/drive.appdata"], "token_info_uri": "https://oauth2.googleapis.com/tokeninfo", "invalid": false, "_class": "OAuth2Credentials", "_module": "oauth2client.client"}
+{"access_token": "ya29.a0AfH6SMA4on6-UlQWL1FDRFX1fhQoNnAQ5v6fnXVb0assyzvK3KabjC3LDymkf7KCo8EFFyL5OIZhQnRYgW0ULFaBtKL_Cw6nOEzabtw0-nUJT7iTN36oANBD1LVNXcArLbQM7teTmjbBS7-KiZRkWOgcPMndbcwZwLlg", "client_id": "710796635688-iivsgbgsb6uv1fap6635dhvuei09o66c.apps.googleusercontent.com", "client_secret": "a1Fz59uTpVNeG_VGuSKDLJXv", "refresh_token": "1//0fpJK6wboLxL7CgYIARAAGA8SNwF-L9IrKmhg8LLnJuWRCiN_YIE88ImD9Y0yV6iCIw220vP3gD4nISe-K5aq8l_r0m1qcUW1D_w", "token_expiry": "2020-07-02T16:37:51Z", "token_uri": "https://oauth2.googleapis.com/token", "user_agent": null, "revoke_uri": "https://oauth2.googleapis.com/revoke", "id_token": null, "id_token_jwt": null, "token_response": {"access_token": "ya29.a0AfH6SMA4on6-UlQWL1FDRFX1fhQoNnAQ5v6fnXVb0assyzvK3KabjC3LDymkf7KCo8EFFyL5OIZhQnRYgW0ULFaBtKL_Cw6nOEzabtw0-nUJT7iTN36oANBD1LVNXcArLbQM7teTmjbBS7-KiZRkWOgcPMndbcwZwLlg", "expires_in": 3599, "scope": "https://www.googleapis.com/auth/drive.appdata https://www.googleapis.com/auth/drive", "token_type": "Bearer"}, "scopes": ["https://www.googleapis.com/auth/drive.appdata", "https://www.googleapis.com/auth/drive"], "token_info_uri": "https://oauth2.googleapis.com/tokeninfo", "invalid": false, "_class": "OAuth2Credentials", "_module": "oauth2client.client"}

setup.py

@@ -30,6 +30,7 @@
       'tqdm',
       'scanimage-tiff-reader',
       'pyqtgraph',
+      'cached-property',
     ],
     tests_require=[
       'pytest',

suite2p/detection/__init__.py

@@ -1,2 +1,2 @@
 from .detect import main_detect
-from .stats import roi_stats
+from .stats import roi_stats, ROI

suite2p/detection/chan2detect.py

@@ -1,6 +1,6 @@
 import numpy as np
 from scipy.ndimage import gaussian_filter
-from .masks import create_cell_masks, create_neuropil_masks, make_masks
+from .masks import create_cell_masks, create_cell_pix, create_neuropil_masks
 
 '''
 identify cells with channel 2 brightness (aka red cells)
@@ -53,10 +53,18 @@ def detect(ops, stats):
 
     # compute pixels in cell and in area around cell (including overlaps)
     # (exclude pixels from other cells)
-    cell_pix, cell_masks0, neuropil_masks = make_masks(ops=ops, stats=stats)
+    cell_pix = create_cell_pix(stats, Ly=ops['Ly'], Lx=ops['Lx'], allow_overlap=ops['allow_overlap'])
+    cell_masks0 = create_cell_masks(stats, Ly=ops['Ly'], Lx=ops['Lx'], allow_overlap=ops['allow_overlap'])
+    neuropil_masks = create_neuropil_masks(
+        ypixs=[stat['ypix'] for stat in stats],
+        xpixs=[stat['xpix'] for stat in stats],
+        cell_pix=cell_pix,
+        inner_neuropil_radius=ops['inner_neuropil_radius'],
+        min_neuropil_pixels=ops['min_neuropil_pixels'],
+    )
     cell_masks = np.zeros((len(stats), Ly * Lx), np.float32)
-    for n in range(len(stats)):
-        cell_masks[n, cell_masks0[n][0]] = cell_masks0[n][1]
+    for cell_mask, cell_mask0 in zip(cell_masks, cell_masks0):
+        cell_mask[cell_mask0[0]] = cell_mask0[1]
 
     inpix = cell_masks @ mimg2.flatten()
     extpix = neuropil_masks @ mimg2.flatten()

suite2p/detection/detect.py

@@ -2,8 +2,8 @@
 import numpy as np
 from pathlib import Path
 from . import sourcery, sparsedetect, chan2detect
-from .stats import roi_stats
-from .masks import get_overlaps, count_overlaps, remove_overlappers, make_masks
+from .stats import ROI
+from .masks import create_cell_masks, create_neuropil_masks, create_cell_pix
 
 
 def main_detect(ops):
@@ -12,20 +12,29 @@ def main_detect(ops):
     else:
         d0 = ops['diameter']
         dy, dx = (d0, d0) if isinstance(d0, int) else d0
-    stat = select_rois(dy=dy, dx=dx, Ly=ops['Ly'], Lx=ops['Lx'], max_overlap=ops['max_overlap'], sparse_mode=ops['sparse_mode'], ops=ops)
+    stats = select_rois(dy=dy, dx=dx, Ly=ops['Ly'], Lx=ops['Lx'], max_overlap=ops['max_overlap'], sparse_mode=ops['sparse_mode'], ops=ops)
     # extract fluorescence and neuropil
     t0 = time.time()
-    cell_pix, cell_masks, neuropil_masks = make_masks(ops, stat)
+    cell_pix = create_cell_pix(stats, Ly=ops['Ly'], Lx=ops['Lx'], allow_overlap=ops['allow_overlap'])
+    cell_masks = create_cell_masks(stats, Ly=ops['Ly'], Lx=ops['Lx'], allow_overlap=ops['allow_overlap'])
+    neuropil_masks = create_neuropil_masks(
+        ypixs=[stat['ypix'] for stat in stats],
+        xpixs=[stat['xpix'] for stat in stats],
+        cell_pix=cell_pix,
+        inner_neuropil_radius=ops['inner_neuropil_radius'],
+        min_neuropil_pixels=ops['min_neuropil_pixels'],
+    )
+
     print('Masks made in %0.2f sec.' % (time.time() - t0))
 
-    ic = np.ones(len(stat), np.bool)
+    ic = np.ones(len(stats), 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)
+        ops, redcell = chan2detect.detect(ops, stats)
         np.save(Path(ops['save_path']).joinpath('redcell.npy'), redcell[ic])
-    return cell_pix, cell_masks, neuropil_masks, stat, ops
+    return cell_pix, cell_masks, neuropil_masks, stats, ops
 
 
 def select_rois(dy: int, dx: int, Ly: int, Lx: int, max_overlap: float, sparse_mode: bool, ops):
@@ -36,21 +45,34 @@ def select_rois(dy: int, dx: int, Ly: int, Lx: int, max_overlap: float, sparse_m
         ops, stats = sourcery.sourcery(ops)
     print('Found %d ROIs, %0.2f sec' % (len(stats), time.time() - t0))
 
-    stats = roi_stats(dy=dy, dx=dx, stats=stats)
+    rois = [ROI(ypix=stat['ypix'], xpix=stat['xpix'], lam=stat['lam'], dx=dx, dy=dy) for stat in stats]
 
-    ypixs = [stat['ypix'] for stat in stats]
-    xpixs = [stat['xpix'] for stat in stats]
-    overlap_masks = get_overlaps(
-        overlaps=count_overlaps(Ly=Ly, Lx=Lx, ypixs=ypixs, xpixs=xpixs),
-        ypixs=ypixs,
-        xpixs=xpixs,
-    )
-    for stat, overlap_mask in zip(stats, overlap_masks):
-        stat['overlap'] = overlap_mask
+    mrs_normeds = ROI.get_mean_r_squared_normed_all(rois=rois)
+    npix_normeds = ROI.get_n_pixels_normed_all(rois=rois)
+    n_overlaps = ROI.get_overlap_count_image(rois=rois, Ly=Ly, Lx=Lx)
+    keep_rois = ROI.filter_overlappers(rois=rois, overlap_image=n_overlaps, max_overlap=max_overlap)
+
+    good_stats = []
+    for keep_roi, roi, mrs_normed, npix_normed, stat in zip(keep_rois, rois, mrs_normeds, npix_normeds, stats):
+        if keep_roi:
+            stat.update({
+                'mrs': mrs_normed,
+                'mrs0': roi.mean_r_squared0,
+                'compact': roi.mean_r_squared_compact,
+                'med': list(roi.median_pix),
+                'npix': roi.n_pixels,
+                'npix_norm': npix_normed,
+                'footprint': 0 if 'footprint' not in stat else stat['footprint'],
+                'overlap': n_overlaps[roi.ypix, roi.xpix] > 1,
+            })
+            if 'radius' not in stat:
+                stat.update({
+                    'radius': roi.radius,
+                    'aspect_ratio': roi.aspect_ratio,
+                })
+            good_stats.append(stat)
 
-    ix = remove_overlappers(ypixs=ypixs, xpixs=xpixs, max_overlap=max_overlap, Ly=Ly, Lx=Lx)
-    stats = [stats[i] for i in ix]
-    print('After removing overlaps, %d ROIs remain' % (len(stats)))
-    return stats
+    print('After removing overlaps, %d ROIs remain' % (len(good_stats)))
+    return good_stats
 
 

suite2p/detection/masks.py

@@ -1,40 +1,31 @@
-from typing import List
+from typing import List, Tuple
 from itertools import count
 import numpy as np
 
 from suite2p.detection.sparsedetect import extendROI
 
-def count_overlaps(Ly: int, Lx: int, ypixs, xpixs) -> np.ndarray:
-    overlap = np.zeros((Ly, Lx))
-    for xpix, ypix in zip(xpixs, ypixs):
-        overlap[ypix, xpix] += 1
-    return overlap
-
-
-def get_overlaps(overlaps, ypixs: List[np.ndarray], xpixs: List[np.ndarray]) -> List[np.ndarray]:
-    """computes overlapping pixels from ROIs"""
-    return [overlaps[ypix, xpix] > 1 for ypix, xpix in zip(ypixs, xpixs)]
 
+def create_cell_pix(stats, Ly, Lx, allow_overlap=False) -> np.ndarray:
+    """Returns Ly x Lx array of whether it contains a cell (1) or not (0)."""
+    cell_pix = np.zeros((Ly, Lx))
+    for stat in stats:
+        mask = ... if allow_overlap else ~stat['overlap']
+        ypix = stat['ypix'][mask]
+        xpix = stat['xpix'][mask]
+        lam = stat['lam'][mask]
+        if xpix.size:
+            cell_pix[ypix[lam > 0], xpix[lam > 0]] = 1
 
-def remove_overlappers(ypixs, xpixs, max_overlap: float, Ly: int, Lx: int) -> List[int]:
-    """returns ROI indices are remain after removing those that overlap more than fraction max_overlap with other ROIs"""
-    overlaps = count_overlaps(Ly=Ly, Lx=Lx, ypixs=ypixs, xpixs=xpixs)
-    ix = []
-    for i, (ypix, xpix) in reversed(list(enumerate(zip(ypixs, xpixs)))):  # todo: is there an ordering effect here that affects which rois will be removed and which will stay?
-        if np.mean(overlaps[ypix, xpix] > 1) > max_overlap:  # note: fancy indexing returns a copy
-            overlaps[ypix, xpix] -= 1
-        else:
-            ix.append(i)
-    return ix[::-1]
+    return cell_pix
 
 
-def create_cell_masks(stat, Ly, Lx, allow_overlap=False):
+def create_cell_masks(stats, Ly, Lx, allow_overlap=False) -> List[Tuple[np.ndarray, np.ndarray]]:
     """ creates cell masks for ROIs in stat and computes radii
 
     Parameters
     ----------
 
-    stat : dictionary
+    stats : dictionary
         'ypix', 'xpix', 'lam'
 
     Ly : float
@@ -48,38 +39,21 @@ def create_cell_masks(stat, Ly, Lx, allow_overlap=False):
 
     Returns
     -------
-    
-    cell_pix : 2D array
-        size [Ly x Lx] where 1 if pixel belongs to cell
-    
     cell_masks : list 
         len ncells, each has tuple of pixels belonging to each cell and weights
 
     """
-
-    ncells = len(stat)
-    cell_pix = np.zeros((Ly,Lx))
     cell_masks = []
+    for stat in stats:
+        mask = ... if allow_overlap else ~stat['overlap']
+        ypix = stat['ypix'][mask]
+        xpix = stat['xpix'][mask]
+        cell_mask = np.ravel_multi_index((ypix, xpix), (Ly, Lx))
+        lam = stat['lam'][mask]
+        lam_normed = lam / lam.sum() if lam.size > 0 else np.empty(0)
+        cell_masks.append((cell_mask, lam_normed))
 
-    for n in range(ncells):
-        if allow_overlap:
-            overlap = np.zeros((stat[n]['npix'],), bool)
-        else:
-            overlap = stat[n]['overlap']
-        ypix = stat[n]['ypix'][~overlap]
-        xpix = stat[n]['xpix'][~overlap]
-        lam  = stat[n]['lam'][~overlap]
-        if xpix.size:
-            # add pixels of cell to cell_pix (pixels to exclude in neuropil computation)
-            cell_pix[ypix[lam>0],xpix[lam>0]] += 1
-            ipix = np.ravel_multi_index((ypix, xpix), (Ly,Lx)).astype('int')
-            #utils.sub2ind((Ly,Lx), ypix, xpix)
-            cell_masks.append((ipix, lam/lam.sum()))
-        else:
-            cell_masks.append((np.zeros(0).astype('int'), np.zeros(0)))
-
-    cell_pix = np.minimum(1, cell_pix)
-    return cell_pix, cell_masks
+    return cell_masks
 
 
 def create_neuropil_masks(ypixs, xpixs, cell_pix, inner_neuropil_radius, min_neuropil_pixels):
@@ -119,18 +93,6 @@ def create_neuropil_masks(ypixs, xpixs, cell_pix, inner_neuropil_radius, min_neu
         neuropil_mask[ypix1[ix], xpix1[ix]] = 1.
         neuropil_mask[ypix, xpix] = 0
 
-    return neuropil_masks / np.sum(neuropil_masks, axis=(1, 2), keepdims=True)
-
-
-def make_masks(ops, stats):
-    Ly, Lx = ops['Ly'], ops['Lx']
-    cell_pix, cell_masks = create_cell_masks(stats, Ly=Ly, Lx=Lx, allow_overlap=ops['allow_overlap'])
-    neuropil_masks = create_neuropil_masks(
-        ypixs=[stat['ypix'] for stat in stats],
-        xpixs=[stat['xpix'] for stat in stats],
-        cell_pix=cell_pix,
-        inner_neuropil_radius=ops['inner_neuropil_radius'],
-        min_neuropil_pixels=ops['min_neuropil_pixels'],
-    )
+    neuropil_masks /= np.sum(neuropil_masks, axis=(1, 2), keepdims=True)
     neuropil_masks = np.reshape(neuropil_masks, (-1, Ly * Lx))
-    return cell_pix, cell_masks, neuropil_masks
+    return neuropil_masks

suite2p/detection/sourcery.py

@@ -6,6 +6,7 @@
 from scipy.ndimage import gaussian_filter
 from matplotlib.colors import hsv_to_rgb
 
+from .stats import fitMVGaus
 from . import utils
 
 
@@ -271,7 +272,7 @@ def get_stat(ops, stats, Ucell, codes, frac=0.5):
         stat['med'] = [np.median(stat['ypix']), np.median(stat['xpix'])]
         stat['npix'] = xpix.size
         if 'radius' not in stat:
-            ry, rx = utils.fitMVGaus(ypix / d0[0], xpix / d0[1], lam, 2).radii
+            ry, rx = fitMVGaus(ypix / d0[0], xpix / d0[1], lam, 2).radii
             stat['radius'] = ry * d0.mean()
             stat['aspect_ratio'] = 2 * ry/(.01 + ry + rx)