Merge pull request #1228 from MouseLand/add_chan2_manual_labeling

Add chan2 manual labeling

Author: carsen-stringer

suite2p/detection/detect.py

@@ -166,8 +166,7 @@ def detection_wrapper(f_reg, diameter=[12., 12.], tau=1., fs=30, meanImg_chan2=N
 
     if mov is None:
         nbins = settings["nbins"]
-        bin_size = int(max(1, n_frames // nbins, np.round(tau * fs)))
-        #bin_size = int(max(1, np.round(tau * fs)))
+        bin_size = settings.get("bin_size") or int(max(1, n_frames // nbins, np.round(tau * fs)))
         logger.info("Binning movie in chunks of %2.2d frames" % bin_size)
         mov = bin_movie(f_reg, bin_size, yrange=yrange, xrange=xrange,
                         badframes=badframes, nbins=nbins)

suite2p/detection/sparsedetect.py

@@ -565,6 +565,8 @@ def sparsery(mov, sdmov, highpass_neuropil,
     scale, estimate_mode = find_best_scale(I=I, spatial_scale=spatial_scale)
 
     spatscale_pix = 3 * 2**scale
+    if isinstance(spatscale_pix, np.ndarray):
+        spatscale_pix = spatscale_pix.item()
     mask_window = int(((spatscale_pix * 1.5) // 2) * 2)
     Th2 = threshold_scaling * 5 * max(
         1, scale)  # threshold for accepted peaks (scale it by spatial scale)

suite2p/gui/drawroi.py

@@ -8,7 +8,7 @@
 import numpy as np
 import pyqtgraph as pg
 from qtpy import QtGui, QtCore
-from qtpy.QtWidgets import QPushButton, QLabel, QLineEdit, QMainWindow, QGridLayout, QButtonGroup, QMessageBox, QWidget
+from qtpy.QtWidgets import QPushButton, QLabel, QLineEdit, QMainWindow, QGridLayout, QButtonGroup, QMessageBox, QWidget, QVBoxLayout
 from matplotlib.colors import hsv_to_rgb
 from scipy import stats
 from scipy.ndimage import rotate
@@ -121,8 +121,10 @@ def __init__(self, bid, Text, parent=None):
         self.show()
 
     def press(self, parent, bid):
+        for b in range(len(parent.views)):
+            parent.viewbtns.button(b).setStyleSheet(parent.styleUnpressed)
+        parent.viewbtns.button(bid).setStyleSheet(parent.stylePressed)
         parent.img0.setImage(parent.masked_images[:, :, :, bid])
-
         parent.win.show()
         parent.show()
 
@@ -141,7 +143,7 @@ def __init__(self, parent):
         # layout = QtGui.QFormLayout()
         self.cwidget.setLayout(self.l0)
         self.stylePressed = ("QPushButton {Text-align: left; "
-                             "background-color: rgb(100,50,100); "
+                             "background-color: rgb(100,100,100); "
                              "color:white;}")
         self.styleUnpressed = ("QPushButton {Text-align: left; "
                                "background-color: rgb(50,50,50); "
@@ -218,17 +220,28 @@ def __init__(self, parent):
             "W: mean img", "E: mean img (enhanced)", "R: correlation map",
             "T: max projection"
         ]
+        self.has_chan2 = "meanImg_chan2" in parent.ops
+        if self.has_chan2:
+            self.views.append("Y: mean img chan2")
         b = 0
         self.viewbtns = QButtonGroup(self)
+        view_container = QWidget()
+        view_vbox = QVBoxLayout()
+        view_vbox.setContentsMargins(0, 0, 0, 0)
+        view_vbox.setSpacing(4)
+        view_container.setLayout(view_vbox)
         for names in self.views:
             btn = ViewButton(b, "&" + names, self)
             self.viewbtns.addButton(btn, b)
-            self.l0.addWidget(btn, b, 4, 1, 1)
+            view_vbox.addWidget(btn)
             btn.setEnabled(True)
             b += 1
-        b = 0
-        self.viewbtns.button(b).setChecked(True)
-        self.viewbtns.button(b).setStyleSheet(self.stylePressed)
+        view_vbox.addStretch()
+        self.l0.addWidget(view_container, 0, 4, 3, 1)
+        for b in range(len(self.views)):
+            self.viewbtns.button(b).setStyleSheet(self.styleUnpressed)
+        self.viewbtns.button(0).setChecked(True)
+        self.viewbtns.button(0).setStyleSheet(self.stylePressed)
 
         self.l0.addWidget(QLabel("neuropil"), 13, 13, 1, 1)
 
@@ -309,37 +322,38 @@ def close_GUI(self):
         self.close()
 
     def normalize_img_add_masks(self):
-        masked_image = np.zeros(
-            ((self.Ly, self.Lx, 3, 4)))  # 3 for RGB and 4 for buttons
-        for i in np.arange(4):  # 4 because 4 buttons
+        nviews = len(self.views)
+        masked_image = np.zeros((self.Ly, self.Lx, 3, nviews))
+        yr = slice(self.parent.ops["yrange"][0], self.parent.ops["yrange"][1])
+        xr = slice(self.parent.ops["xrange"][0], self.parent.ops["xrange"][1])
+        for i in np.arange(nviews):
             if i == 0:
-                mimg = np.zeros((self.Ly, self.Lx), np.float32)
-                mimg[self.parent.ops["yrange"][0]:self.parent.ops["yrange"][1],
-                     self.parent.ops["xrange"][0]:self.parent.ops["xrange"][1]] = self.parent.ops["meanImg"][
-                         self.parent.ops["yrange"][0]:self.parent.ops["yrange"][1],
-                         self.parent.ops["xrange"][0]:self.parent.ops["xrange"][1]]
-
+                src = self.parent.ops["meanImg"]
             elif i == 1:
-                mimg = np.zeros((self.Ly, self.Lx), np.float32)
-                mimg[self.parent.ops["yrange"][0]:self.parent.ops["yrange"][1],
-                     self.parent.ops["xrange"][0]:self.parent.ops["xrange"][1]] = self.parent.ops["meanImgE"][
-                         self.parent.ops["yrange"][0]:self.parent.ops["yrange"][1],
-                         self.parent.ops["xrange"][0]:self.parent.ops["xrange"][1]]
+                src = self.parent.ops["meanImgE"]
             elif i == 2:
-                mimg = np.zeros((self.Ly, self.Lx), np.float32)
-                mimg[self.parent.ops["yrange"][0]:self.parent.ops["yrange"][1],
-                     self.parent.ops["xrange"][0]:self.parent.ops["xrange"][1]] = self.parent.ops["Vcorr"]
-
+                src = self.parent.ops["Vcorr"]
+            elif i == 3:
+                src = self.parent.ops.get("max_proj", None)
+            elif i == 4 and self.has_chan2:
+                src = self.parent.ops["meanImg_chan2"]
             else:
-                mimg = np.zeros((self.Ly, self.Lx), np.float32)
-                if "max_proj" in self.parent.ops:
-                    mimg[self.parent.ops["yrange"][0]:self.parent.ops["yrange"][1],
-                         self.parent.ops["xrange"][0]:self.parent.ops["xrange"][1]] = self.parent.ops["max_proj"]
-
-            mimg1 = np.percentile(mimg, 1)
-            mimg99 = np.percentile(mimg, 99)
-            mimg = (mimg - mimg1) / (mimg99 - mimg1)
-            mimg = np.maximum(0, np.minimum(1, mimg))
+                src = None
+
+            mimg = np.zeros((self.Ly, self.Lx), np.float32)
+            if src is not None:
+                mimg1 = np.percentile(src, 1)
+                mimg99 = np.percentile(src, 99)
+                if mimg99 > mimg1:
+                    src = (src - mimg1) / (mimg99 - mimg1)
+                else:
+                    src = np.zeros_like(src)
+                src = np.clip(src, 0, 1).astype(np.float32)
+                if src.shape[0] == self.Ly and src.shape[1] == self.Lx:
+                    mimg = src
+                else:
+                    mimg[yr, xr] = src
+
             masked_image[:, :, :, i] = self.create_masks_of_cells(mimg)
 
         return masked_image
@@ -387,6 +401,9 @@ def keyPressEvent(self, event):
             elif event.key() == QtCore.Qt.Key_T:
                 self.viewbtns.button(3).setChecked(True)
                 self.viewbtns.button(3).press(self, 3)
+            elif event.key() == QtCore.Qt.Key_Y and self.has_chan2:
+                self.viewbtns.button(4).setChecked(True)
+                self.viewbtns.button(4).press(self, 4)
 
     def add_ROI(self, pos=None):
         self.iROI = len(self.ROIs)
@@ -557,13 +574,13 @@ def remove(self, parent):
         parent.win.show()
         parent.show()
 
-    def rotate_ROI(self, parent, ellipse, xrange, yrange, posx, posy):
+    def rotate_ROI(self, parent, ellipse, xrange, yrange, center_x, center_y):
         #Rotates ROI depending on Rotatehandle degree
         ellipse = rotate(ellipse, angle=math.floor(self.ROI.angle()), order=0)
         ellipse = np.flip(ellipse, axis=0)
-        xrange = (np.arange(-1 * int(ellipse.shape[1] - 1), 1) + int(posx)).astype(np.int32)
-        yrange = (np.arange(-1 * int(ellipse.shape[0] - 1), 1) + int(posy)).astype(np.int32)
-        yrange += int(np.floor(ellipse.shape[0] / 2)) + 1
+        w, h = ellipse.shape[1], ellipse.shape[0]
+        xrange = (np.arange(-(w // 2), w - w // 2) + int(center_x)).astype(np.int32)
+        yrange = (np.arange(-(h // 2), h - h // 2) + int(center_y)).astype(np.int32)
         return ellipse, xrange, yrange
 
     def position(self, parent):
@@ -579,13 +596,16 @@ def position(self, parent):
         yrange = (np.arange(-1 * int(sizey), 1) + int(posy)).astype(np.int32)
         yrange += int(np.floor(sizey / 2)) + 1
         # what is ellipse circling?
-        br = self.ROI.boundingRect()
         ellipse = np.zeros((yrange.size, xrange.size), "bool")
         x, y = np.meshgrid(np.arange(0, xrange.size, 1), np.arange(0, yrange.size, 1))
         ellipse = ((y - br.center().y())**2 / (br.height() / 2)**2 +
                    (x - br.center().x())**2 / (br.width() / 2)**2) <= 1
+        center_scene = self.ROI.mapToScene(br.center())
+        center_view = parent.p0.mapSceneToView(center_scene)
+        center_x = center_view.x()
+        center_y = center_view.y()
         if self.ROI.angle() not in (0, 180, -180):
-            ellipse, xrange, yrange = self.rotate_ROI(parent, ellipse, xrange, yrange, posx, posy)
+            ellipse, xrange, yrange = self.rotate_ROI(parent, ellipse, xrange, yrange, center_x, center_y)
         #ensures that ROI is not placed outside of movie coordinates
         ellipse = ellipse[:, np.logical_and(xrange >= 0, xrange < parent.Lx)]
         xrange = xrange[np.logical_and(xrange >= 0, xrange < parent.Lx)]

suite2p/gui/io.py

@@ -457,29 +457,12 @@ def save_iscell(parent):
 
 def save_mat(parent):
     print("saving to mat")
-    matpath = os.path.join(parent.basename, "Fall.mat")
-    if "date_proc" in parent.ops:
-        parent.ops["date_proc"] = []
-    scipy.io.savemat(
-        matpath, {
-            "stat":
-                parent.stat,
-            "settings":
-                parent.ops,
-            "F":
-                parent.Fcell,
-            "Fneu":
-                parent.Fneu,
-            "spks":
-                parent.Spks,
-            "iscell":
-                np.concatenate(
-                    (parent.iscell[:, np.newaxis], parent.probcell[:, np.newaxis]),
-                    axis=1),
-            "redcell":
-                np.concatenate((np.expand_dims(parent.redcell, axis=1),
-                                np.expand_dims(parent.probredcell, axis=1)), axis=1)
-        })
+    iscell = np.concatenate(
+        (parent.iscell[:, np.newaxis], parent.probcell[:, np.newaxis]), axis=1)
+    redcell = np.concatenate(
+        (np.expand_dims(parent.redcell, axis=1), np.expand_dims(parent.probredcell, axis=1)), axis=1)
+    ops = {**parent.ops, "save_path": parent.basename}
+    io.save_mat(ops, parent.stat, parent.Fcell, parent.Fneu, parent.Spks, iscell, redcell)
 
 
 def save_merge(parent):

suite2p/io/h5.py

@@ -82,7 +82,7 @@ def h5py_to_binary(dbs, settings, reg_file, reg_file_chan2):
                     else:
                         irange = np.arange(
                             ik / ncp, min(ik / ncp + nbatch / ncp, nframes_all / ncp),
-                            1)
+                            1).astype(int)
                         if irange.size == 0:
                             break
                         im = f[key][irange, ...]

suite2p/io/save.py

@@ -60,11 +60,15 @@ def save_mat(ops, stat, F, Fneu, spks, iscell, redcell,
     stat = np.array(stat, dtype=object)
 
 
-    # Check for None values in ops_matlab and replace with empty arrays
-    for k, v in ops_matlab.items():
-        if v is None:
-            logger.warning(f"ops_matlab['{k}'] is None, replacing with empty array")
-            ops_matlab[k] = np.array([])
+    # Check for None values in ops_matlab and replace with empty arrays in recursive manner
+    def replace_none(d):
+        for k, v in d.items():
+            if v is None:
+                logger.warning(f"'{k}' is None, replacing with empty array")
+                d[k] = np.array([])
+            elif isinstance(v, dict):
+                replace_none(v)
+    replace_none(ops_matlab)
 
     # Handle None variables by replacing with empty arrays
     if redcell is None:
@@ -201,6 +205,12 @@ def combined(save_folder, save=True):
     plane_folders = natsorted([
         f.path for f in os.scandir(save_folder) if f.is_dir() and f.name[:5] == "plane"
     ])
+    top_db_path = os.path.join(save_folder, "db.npy")
+    if os.path.exists(top_db_path):
+        top_db = np.load(top_db_path, allow_pickle=True).item()
+        ignore_flyback = set(top_db.get("ignore_flyback") or [])
+        if ignore_flyback:
+            plane_folders = [f for i, f in enumerate(plane_folders) if i not in ignore_flyback]
     dbs = [
         np.load(os.path.join(f, "db.npy"), allow_pickle=True).item()
         for f in plane_folders

suite2p/parameters.py

@@ -544,7 +544,7 @@
             "min": 1,
             "max": np.inf,
             "default": None,
-            "description": "Size of bins for cell detection (default is tau * fs).",
+            "description": "Size of bins for cell detection. If None, computed as max(1, nframes // nbins, round(tau * fs)), where nframes is the total number of frames in the registered movie.",
         },
         "highpass_time": {
             "gui_name": "Highpass time",

suite2p/registration/register.py

@@ -487,11 +487,11 @@ def normalize_reference_image(refImg):
     return refImg, rmin, rmax
 
 
-def register_frames(f_align_in, refImg, f_align_out=None, batch_size=100, 
-                    bidiphase=0, norm_frames=True, smooth_sigma=1.15, spatial_taper=3.45, 
-                    block_size=(128,128), nonrigid=True, maxregshift=0.1, 
+def register_frames(f_align_in, refImg, f_align_out=None, batch_size=100,
+                    bidiphase=0, norm_frames=True, smooth_sigma=1.15, spatial_taper=3.45,
+                    block_size=(128,128), nonrigid=True, maxregshift=0.1,
                     smooth_sigma_time=0, snr_thresh=1.2, maxregshiftNR=5,
-                    device=torch.device("cuda"), tif_root=None, apply_shifts=True,
+                    subpixel=10, device=torch.device("cuda"), tif_root=None, apply_shifts=True,
                     upsample_meanImg=False):
     """
     Register frames to a reference image using rigid and optionally nonrigid shifts.
@@ -578,11 +578,11 @@ def register_frames(f_align_in, refImg, f_align_out=None, batch_size=100,
     else:
         nZ = 1
 
-    refAndMasks = compute_filters_and_norm(refImg, norm_frames=norm_frames, 
+    refAndMasks = compute_filters_and_norm(refImg, norm_frames=norm_frames,
                                            spatial_smooth=smooth_sigma,
-                                           spatial_taper=spatial_taper, 
-                                           block_size=block_size if nonrigid else None, 
-                                           device=device)
+                                           spatial_taper=spatial_taper,
+                                           block_size=block_size if nonrigid else None,
+                                           subpixel=subpixel, device=device)
     blocks = refAndMasks[-3] if nZ==1 else refAndMasks[0][-3]
     rmin = refAndMasks[-2] if nZ==1 else [refAndMasks[z][-2] for z in range(nZ)]
     rmax = refAndMasks[-1] if nZ==1 else [refAndMasks[z][-1] for z in range(nZ)]
@@ -829,8 +829,8 @@ def assign_reg_io(f_reg, f_raw=None, f_reg_chan2=None,
         if f_align_in.shape[0] != f_alt_in.shape[0]:
             raise ValueError("number of frames in f_align_in and f_alt_in must match")
 
+    tif_root_align, tif_root_alt = None, None
     if save_path:
-        tif_root_align, tif_root_alt = None, None
         if reg_tif:
             tifroot = os.path.join(save_path, "reg_tif")
             os.makedirs(tifroot, exist_ok=True)
@@ -954,8 +954,9 @@ def registration_wrapper(f_reg, f_raw=None, f_reg_chan2=None, f_raw_chan2=None,
                                 spatial_taper=settings["spatial_taper"], block_size=settings["block_size"],
                                 nonrigid=settings["nonrigid"],
                                 maxregshift=settings["maxregshift"], smooth_sigma_time=settings["smooth_sigma_time"],
-                                    snr_thresh=settings["snr_thresh"], maxregshiftNR=settings["maxregshiftNR"],
-                                    device=device, upsample_meanImg=settings.get("upsample_meanImg", False))
+                                snr_thresh=settings["snr_thresh"], maxregshiftNR=settings["maxregshiftNR"],
+                                subpixel=settings["subpixel"],
+                                device=device, upsample_meanImg=settings.get("upsample_meanImg", False))
         rmin, rmax, mean_img, offsets_all, blocks, mean_img_ups, counts_ups, meanImg_ups = outputs
         yoff, xoff, corrXY, yoff1, xoff1, corrXY1, zest, cmax_all = offsets_all
 

suite2p/run_s2p.py

@@ -378,6 +378,9 @@ def run_s2p(db={}, settings=default_settings(), server={}):
         plane_folders = natsorted([
             f.path for f in os.scandir(save_folder) if f.is_dir() and f.name[:5] == "plane"
         ])
+        ignore_flyback = set(db.get("ignore_flyback") or [])
+        if ignore_flyback:
+            plane_folders = [f for i, f in enumerate(plane_folders) if i not in ignore_flyback]
     else:
         plane_folders = []