Skip to content

suite2p.detection package#

Copyright © 2023 Howard Hughes Medical Institute, Authored by Carsen Stringer and Marius Pachitariu.

mask_centers #

mask_centers(masks)

Compute the centers and diameters for all masks in a label image.

Parameters:

Name Type Description Default
masks ndarray

2D integer label image where each unique positive value is an ROI.

 required

Returns:

Name Type Description
centers ndarray

Median center coordinates of shape (n_masks, 2), as [ymed, xmed].
diams ndarray

Estimated diameters for each mask, shape (n_masks,).
Source code in suite2p/detection/anatomical.py 
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
def mask_centers(masks):
    """
    Compute the centers and diameters for all masks in a label image.

    Parameters
    ----------
    masks : numpy.ndarray
        2D integer label image where each unique positive value is an ROI.

    Returns
    -------
    centers : numpy.ndarray
        Median center coordinates of shape (n_masks, 2), as [ymed, xmed].
    diams : numpy.ndarray
        Estimated diameters for each mask, shape (n_masks,).
    """
    centers = np.zeros((masks.max(), 2), np.int32)
    diams = np.zeros(masks.max(), np.float32)
    slices = find_objects(masks)
    for i, si in enumerate(slices):
        if si is not None:
            sr, sc = si
            ymed, xmed, diam = mask_stats(masks[sr, sc] == (i + 1))
            centers[i] = np.array([ymed, xmed])
            diams[i] = diam
    return centers, diams

mask_stats #

mask_stats(mask)

Compute the median center and diameter of a single binary mask.

Parameters:

Name Type Description Default
mask ndarray

2D binary mask for a single ROI.

 required

Returns:

Name Type Description
ymed int

Y-coordinate of the pixel closest to the median center.
xmed int

X-coordinate of the pixel closest to the median center.
diam float

Estimated diameter of the mask, computed from the number of pixels.
Source code in suite2p/detection/anatomical.py 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
def mask_stats(mask):
    """
    Compute the median center and diameter of a single binary mask.

    Parameters
    ----------
    mask : numpy.ndarray
        2D binary mask for a single ROI.

    Returns
    -------
    ymed : int
        Y-coordinate of the pixel closest to the median center.
    xmed : int
        X-coordinate of the pixel closest to the median center.
    diam : float
        Estimated diameter of the mask, computed from the number of pixels.
    """
    y, x = np.nonzero(mask)
    y = y.astype(np.int32)
    x = x.astype(np.int32)
    ymed = np.median(y)
    xmed = np.median(x)
    imin = ((x - xmed)**2 + (y - ymed)**2).argmin()
    xmed = x[imin]
    ymed = y[imin]
    diam = len(y)**0.5
    diam /= (np.pi**0.5) / 2
    return ymed, xmed, diam

masks_to_stats #

masks_to_stats(masks, weights)

Convert a label image and weight map into an array of ROI statistics dictionaries.

For each mask, extracts the pixel coordinates, pixel weights from the weight image, and computes the median center.

Parameters:

Name Type Description Default
masks ndarray

2D integer label image where each unique positive value is an ROI.

 required
weights ndarray

2D weight image of the same shape as masks (e.g. max projection), used to assign pixel weights ("lam") to each ROI.

 required

Returns:

Name Type Description
stats ndarray

Array of dictionaries, one per ROI, each containing "ypix", "xpix", "lam", "med", and "footprint".
Source code in suite2p/detection/anatomical.py 
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
def masks_to_stats(masks, weights):
    """
    Convert a label image and weight map into an array of ROI statistics dictionaries.

    For each mask, extracts the pixel coordinates, pixel weights from the weight
    image, and computes the median center.

    Parameters
    ----------
    masks : numpy.ndarray
        2D integer label image where each unique positive value is an ROI.
    weights : numpy.ndarray
        2D weight image of the same shape as `masks` (e.g. max projection),
        used to assign pixel weights ("lam") to each ROI.

    Returns
    -------
    stats : numpy.ndarray
        Array of dictionaries, one per ROI, each containing "ypix", "xpix",
        "lam", "med", and "footprint".
    """
    stats = []
    slices = find_objects(masks)
    for i, si in enumerate(slices):
        sr, sc = si
        ypix0, xpix0 = np.nonzero(masks[sr, sc] == (i + 1))
        ypix0 = ypix0.astype(int) + sr.start
        xpix0 = xpix0.astype(int) + sc.start
        ymed = np.median(ypix0)
        xmed = np.median(xpix0)
        imin = np.argmin((xpix0 - xmed)**2 + (ypix0 - ymed)**2)
        xmed = xpix0[imin]
        ymed = ypix0[imin]
        stats.append({
            "ypix": ypix0,
            "xpix": xpix0,
            "lam": weights[ypix0, xpix0],
            "med": [ymed, xmed],
            "footprint": 1
        })
    stats = np.array(stats)
    return stats

roi_detect #

roi_detect(mproj, diameter=None, settings=None, pretrained_model=None, device=torch.device('cuda'), chan2=False)

Detect ROIs in an image using Cellpose.

Runs a Cellpose model on the input image to segment cells. Optionally rescales the image if the diameter aspect ratio is not 1.

Parameters:

Name Type Description Default
mproj ndarray

2D image to segment, shape (Ly, Lx).

 required
diameter float or list of float

Expected cell diameter in pixels. If scalar, used for both axes. If list, [dy, dx]. Defaults to [30., 30.].

 None
settings dict

Detection settings dictionary. Used to get "params", "chan2_params", "cellprob_threshold", and "flow_threshold" for Cellpose.

 None
pretrained_model str

Name of the Cellpose pretrained model. Defaults to "cpsam".

 None
device torch.device, optional (default torch.device("cuda"))

Torch device, used for GPU cache cleanup after detection.

 device('cuda')
chan2 bool, optional (default False)

If True, use "chan2_params" from settings instead of "params".

 False

Returns:

Name Type Description
masks ndarray

Integer label image of shape (Ly, Lx), where each ROI has a unique positive integer label.
centers ndarray

Median center coordinates of shape (n_masks, 2).
median_diam float

Median diameter across all detected masks.
mask_diams ndarray

Diameters for each detected mask, shape (n_masks,), dtype int32.
Source code in suite2p/detection/anatomical.py 
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
def roi_detect(mproj, diameter=None, settings=None,
               pretrained_model=None, device=torch.device("cuda"), chan2=False):
    """
    Detect ROIs in an image using Cellpose.

    Runs a Cellpose model on the input image to segment cells. Optionally
    rescales the image if the diameter aspect ratio is not 1.

    Parameters
    ----------
    mproj : numpy.ndarray
        2D image to segment, shape (Ly, Lx).
    diameter : float or list of float, optional
        Expected cell diameter in pixels. If scalar, used for both axes.
        If list, [dy, dx]. Defaults to [30., 30.].
    settings : dict, optional
        Detection settings dictionary. Used to get "params", "chan2_params",
        "cellprob_threshold", and "flow_threshold" for Cellpose.
    pretrained_model : str, optional
        Name of the Cellpose pretrained model. Defaults to "cpsam".
    device : torch.device, optional (default torch.device("cuda"))
        Torch device, used for GPU cache cleanup after detection.
    chan2 : bool, optional (default False)
        If True, use "chan2_params" from settings instead of "params".

    Returns
    -------
    masks : numpy.ndarray
        Integer label image of shape (Ly, Lx), where each ROI has a
        unique positive integer label.
    centers : numpy.ndarray
        Median center coordinates of shape (n_masks, 2).
    median_diam : float
        Median diameter across all detected masks.
    mask_diams : numpy.ndarray
        Diameters for each detected mask, shape (n_masks,), dtype int32.
    """
    Lyc, Lxc = mproj.shape
    diameter = [30., 30. ] if diameter is None else diameter
    diameter = [diameter, diameter] if np.isscalar(diameter) else diameter

    rescale = diameter[1] / diameter[0]
    if rescale != 1.0:
        img = cv2.resize(img, (Lxc, int(Lyc * rescale)))
    logger.info("!NOTE! diameter set to %0.2f for cell detection with cellpose" %
                diameter[1])

    pretrained_model = "cpsam" if pretrained_model is None else pretrained_model
    model = CellposeModel(pretrained_model=pretrained_model, gpu=True if core.use_gpu() else False)
    params = settings["params"] if not chan2 else settings["chan2_params"]
    params = {} if params is None else params
    masks = model.eval(mproj, diameter=diameter[1],
                       cellprob_threshold=settings.get("cellprob_threshold", 0.0),
                       flow_threshold=settings.get("flow_threshold", 0.4),
                       **params)[0]
    shape = masks.shape
    _, masks = np.unique(np.int32(masks), return_inverse=True)
    masks = masks.reshape(shape)

    if rescale != 1.0:
        masks = cv2.resize(masks, (Lxc, Lyc), interpolation=cv2.INTER_NEAREST)
        img = cv2.resize(img, (Lxc, Lyc))

    centers, mask_diams = mask_centers(masks)
    median_diam = np.median(mask_diams)
    logger.info(">>>> %d masks detected, median diameter = %0.2f " %
          (masks.max(), median_diam))

    if device.type == "cuda":
        del model
        torch.cuda.empty_cache()

    return masks, centers, median_diam, mask_diams.astype(np.int32)

select_rois #

select_rois(mean_img, max_proj, settings, diameter=[12.0, 12.0], device=torch.device('cuda'))

Find ROIs in static images using Cellpose anatomical detection.

Prepares an image for segmentation based on the "img" setting (max projection / mean image ratio, mean image, or max projection), optionally applies spatial high-pass filtering, then runs Cellpose to detect ROIs.

Parameters:

Name Type Description Default
mean_img ndarray

2D mean image of shape (Ly, Lx).

 required
max_proj ndarray

2D maximum projection image of shape (Ly, Lx).

 required
settings dict

Detection settings dictionary. Must contain "img" (str specifying which image to segment) and optionally "highpass_spatial" (float).

 required
diameter list of float, optional (default [12., 12.])

Expected cell diameter [dy, dx] in pixels.

 [12.0, 12.0]
device torch.device, optional (default torch.device("cuda"))

Torch device for Cellpose and GPU cache cleanup.

 device('cuda')

Returns:

Name Type Description
new_settings dict

Dictionary with detection metadata including "diameter", "Vcorr", and placeholder keys "Vmax", "ihop", "Vsplit", "Vmap", "spatscale_pix".
stats ndarray

Array of ROI statistics dictionaries, each containing "ypix", "xpix", "lam", "med", and "footprint".
Source code in suite2p/detection/anatomical.py 
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
def select_rois(mean_img, max_proj, settings, 
                diameter=[12., 12.],
                device=torch.device("cuda")):
    """
    Find ROIs in static images using Cellpose anatomical detection.

    Prepares an image for segmentation based on the "img" setting (max
    projection / mean image ratio, mean image, or max projection), optionally
    applies spatial high-pass filtering, then runs Cellpose to detect ROIs.

    Parameters
    ----------
    mean_img : numpy.ndarray
        2D mean image of shape (Ly, Lx).
    max_proj : numpy.ndarray
        2D maximum projection image of shape (Ly, Lx).
    settings : dict
        Detection settings dictionary. Must contain "img" (str specifying
        which image to segment) and optionally "highpass_spatial" (float).
    diameter : list of float, optional (default [12., 12.])
        Expected cell diameter [dy, dx] in pixels.
    device : torch.device, optional (default torch.device("cuda"))
        Torch device for Cellpose and GPU cache cleanup.

    Returns
    -------
    new_settings : dict
        Dictionary with detection metadata including "diameter", "Vcorr",
        and placeholder keys "Vmax", "ihop", "Vsplit", "Vmap",
        "spatscale_pix".
    stats : numpy.ndarray
        Array of ROI statistics dictionaries, each containing "ypix",
        "xpix", "lam", "med", and "footprint".
    """
    Lyc, Lxc = mean_img.shape
    if settings["img"] == 'max_proj / meanImg':
        img = np.log(np.maximum(1e-3, max_proj / np.maximum(1e-3, mean_img)))
        weights = max_proj
    elif settings["img"] == 'meanImg':
        img = mean_img
        weights = 0.1 + np.clip(
            (mean_img - np.percentile(mean_img, 1)) /
            (np.percentile(mean_img, 99) - np.percentile(mean_img, 1)), 0, 1)
    else:
        img = max_proj.copy()
        weights = max_proj

    t0 = time.time()

    if settings.get("highpass_spatial", 0):
        img = np.clip(normalize99(img), 0, 1)
        img -= gaussian_filter(img, diameter[1] * settings["highpass_spatial"])
        img -= gaussian_filter(img, diameter[1] * settings["highpass_spatial"])

    masks, centers, median_diam, mask_diams = roi_detect(
        img, diameter=diameter, settings=settings, device=device)

    stats = masks_to_stats(masks, weights)
    logger.info("Detected %d ROIs, %0.2f sec" % (len(stats), time.time() - t0))

    new_settings = {
        "diameter": median_diam,
        "Vmax": 0,
        "ihop": 0,
        "Vsplit": 0,
        "Vcorr": img,
        "Vmap": 0,
        "spatscale_pix": 0
    }

    return new_settings, stats

Copyright © 2023 Howard Hughes Medical Institute, Authored by Carsen Stringer and Marius Pachitariu.

cellpose_overlap #

cellpose_overlap(stats, mimg2, diameter, chan2_threshold=0.25, device=torch.device('cuda'), settings=None)

Classify cells as red by computing overlap with Cellpose-detected masks.

Runs Cellpose on the red channel mean image to detect anatomical masks, then computes the intersection-over-union (IOU) between each ROI and the Cellpose masks. ROIs with IOU above the threshold are classified as red.

Parameters:

Name Type Description Default
stats ndarray

Array of ROI statistics dictionaries, each containing "ypix" and "xpix".

 required
mimg2 ndarray

Red channel mean image of shape (Ly, Lx).

 required
diameter float or list of float

Expected cell diameter in pixels for Cellpose detection.

 required
chan2_threshold float, optional (default 0.25)

IOU threshold for red cell classification.

 0.25
device torch.device, optional (default torch.device("cuda"))

Torch device for Cellpose and GPU cache cleanup.

 device('cuda')
settings dict

Detection settings dictionary passed to Cellpose.

 None

Returns:

Name Type Description
redstats ndarray

Array of shape (n_cells, 2) where column 0 is the binary red cell label and column 1 is the maximum IOU with Cellpose masks.
masks ndarray

Integer label image of Cellpose-detected masks in the red channel, shape (Ly, Lx).
Source code in suite2p/detection/chan2detect.py 
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
def cellpose_overlap(stats, mimg2, diameter, chan2_threshold=0.25, device=torch.device("cuda"),
                     settings=None):
    """
    Classify cells as red by computing overlap with Cellpose-detected masks.

    Runs Cellpose on the red channel mean image to detect anatomical masks,
    then computes the intersection-over-union (IOU) between each ROI and the
    Cellpose masks. ROIs with IOU above the threshold are classified as red.

    Parameters
    ----------
    stats : numpy.ndarray
        Array of ROI statistics dictionaries, each containing "ypix" and "xpix".
    mimg2 : numpy.ndarray
        Red channel mean image of shape (Ly, Lx).
    diameter : float or list of float
        Expected cell diameter in pixels for Cellpose detection.
    chan2_threshold : float, optional (default 0.25)
        IOU threshold for red cell classification.
    device : torch.device, optional (default torch.device("cuda"))
        Torch device for Cellpose and GPU cache cleanup.
    settings : dict, optional
        Detection settings dictionary passed to Cellpose.

    Returns
    -------
    redstats : numpy.ndarray
        Array of shape (n_cells, 2) where column 0 is the binary red cell
        label and column 1 is the maximum IOU with Cellpose masks.
    masks : numpy.ndarray
        Integer label image of Cellpose-detected masks in the red channel,
        shape (Ly, Lx).
    """
    from . import anatomical
    masks = anatomical.roi_detect(mimg2, diameter=diameter, device=device, settings=settings, chan2=True)[0]
    Ly, Lx = masks.shape
    redstats = np.zeros((len(stats), 2),
                        np.float32)  #changed the size of preallocated space
    for i in range(len(stats)):
        smask = np.zeros((Ly, Lx), np.uint16)
        ypix0, xpix0 = stats[i]["ypix"], stats[i]["xpix"]
        smask[ypix0, xpix0] = 1
        ious = utils.mask_ious(masks, smask)[0]
        if ious.size > 0:
            iou = ious.max()
        else:
            iou = 0.0
        redstats[
            i,
        ] = np.array([iou > chan2_threshold, iou])  #this had the wrong dimension
    return redstats, masks

correct_bleedthrough #

correct_bleedthrough(Ly, Lx, nblks, mimg, mimg2)

Subtract bleedthrough of the green channel into the red channel.

Uses non-rigid regression with nblks x nblks spatial blocks to estimate and remove the green-to-red bleedthrough, producing a corrected red channel mean image.

Parameters:

Name Type Description Default
Ly int

Height of the image in pixels.

 required
Lx int

Width of the image in pixels.

 required
nblks int

Number of spatial blocks along each axis for piecewise regression.

 required
mimg ndarray

Green channel mean image of shape (Ly, Lx).

 required
mimg2 ndarray

Red channel mean image of shape (Ly, Lx).

 required

Returns:

Name Type Description
mimg2 ndarray

Corrected red channel mean image with bleedthrough subtracted, clipped to non-negative values, shape (Ly, Lx).
Source code in suite2p/detection/chan2detect.py 
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
def correct_bleedthrough(Ly, Lx, nblks, mimg, mimg2):
    """
    Subtract bleedthrough of the green channel into the red channel.

    Uses non-rigid regression with nblks x nblks spatial blocks to estimate
    and remove the green-to-red bleedthrough, producing a corrected red
    channel mean image.

    Parameters
    ----------
    Ly : int
        Height of the image in pixels.
    Lx : int
        Width of the image in pixels.
    nblks : int
        Number of spatial blocks along each axis for piecewise regression.
    mimg : numpy.ndarray
        Green channel mean image of shape (Ly, Lx).
    mimg2 : numpy.ndarray
        Red channel mean image of shape (Ly, Lx).

    Returns
    -------
    mimg2 : numpy.ndarray
        Corrected red channel mean image with bleedthrough subtracted,
        clipped to non-negative values, shape (Ly, Lx).
    """
    sT = np.round((Ly + Lx) / (nblks * 2) * 0.25)
    mask = np.zeros((Ly, Lx, nblks, nblks), np.float32)
    weights = np.zeros((nblks, nblks), np.float32)
    yb = np.linspace(0, Ly, nblks + 1).astype(int)
    xb = np.linspace(0, Lx, nblks + 1).astype(int)
    for iy in range(nblks):
        for ix in range(nblks):
            ny = np.arange(yb[iy], yb[iy + 1]).astype(int)
            nx = np.arange(xb[ix], xb[ix + 1]).astype(int)
            mask[:, :, iy, ix] = quadrant_mask(Ly, Lx, ny, nx, sT)
            x = mimg[np.ix_(ny, nx)].flatten()
            x2 = mimg2[np.ix_(ny, nx)].flatten()
            # predict chan2 from chan1
            a = (x * x2).sum() / (x * x).sum()
            weights[iy, ix] = a
    mask /= mask.sum(axis=-1).sum(axis=-1)[:, :, np.newaxis, np.newaxis]
    mask *= weights
    mask *= mimg[:, :, np.newaxis, np.newaxis]
    mimg2 -= mask.sum(axis=-1).sum(axis=-1)
    mimg2 = np.maximum(0, mimg2)
    return mimg2

detect #

detect(meanImg, meanImg_chan2, stats, diameter, cellpose_chan2=False, chan2_threshold=0.65, device=torch.device('cuda'), settings=None, inner_neuropil_radius=2, min_neuropil_pixels=350)

Identify red cells using the second channel (e.g. tdTomato).

Attempts Cellpose-based overlap detection first if cellpose_chan2 is True. Falls back to intensity-ratio detection if Cellpose fails or is disabled.

Parameters:

Name Type Description Default
meanImg ndarray

Green channel mean image of shape (Ly, Lx).

 required
meanImg_chan2 ndarray

Red channel mean image of shape (Ly, Lx).

 required
stats ndarray

Array of ROI statistics dictionaries, each containing "ypix", "xpix", and "lam".

 required
diameter float or list of float

Expected cell diameter in pixels for Cellpose detection.

 required
cellpose_chan2 bool, optional (default False)

If True, attempt Cellpose-based red cell detection first.

 False
chan2_threshold float, optional (default 0.65)

Threshold for red cell classification. Meaning depends on method: IOU threshold for Cellpose, intensity ratio for fallback.

 0.65
device torch.device, optional (default torch.device("cuda"))

Torch device for Cellpose and GPU cache cleanup.

 device('cuda')
settings dict

Detection settings dictionary passed to Cellpose.

 None
inner_neuropil_radius int, optional (default 2)

Radius in pixels of the inner exclusion zone for neuropil masks, used in intensity-ratio fallback.

 2
min_neuropil_pixels int, optional (default 350)

Minimum number of neuropil pixels, used in intensity-ratio fallback.

 350

Returns:

Name Type Description
masks ndarray or None

Cellpose-detected masks in the red channel if Cellpose was used, otherwise None.
redstats ndarray

Array of shape (n_cells, 2) where column 0 is the binary red cell label and column 1 is the red probability.
Source code in suite2p/detection/chan2detect.py 
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
def detect(meanImg, meanImg_chan2, stats, diameter, cellpose_chan2=False, chan2_threshold=0.65,
           device=torch.device("cuda"), settings=None, inner_neuropil_radius=2,
           min_neuropil_pixels=350):
    """
    Identify red cells using the second channel (e.g. tdTomato).

    Attempts Cellpose-based overlap detection first if ``cellpose_chan2`` is
    True. Falls back to intensity-ratio detection if Cellpose fails or is
    disabled.

    Parameters
    ----------
    meanImg : numpy.ndarray
        Green channel mean image of shape (Ly, Lx).
    meanImg_chan2 : numpy.ndarray
        Red channel mean image of shape (Ly, Lx).
    stats : numpy.ndarray
        Array of ROI statistics dictionaries, each containing "ypix", "xpix",
        and "lam".
    diameter : float or list of float
        Expected cell diameter in pixels for Cellpose detection.
    cellpose_chan2 : bool, optional (default False)
        If True, attempt Cellpose-based red cell detection first.
    chan2_threshold : float, optional (default 0.65)
        Threshold for red cell classification. Meaning depends on method:
        IOU threshold for Cellpose, intensity ratio for fallback.
    device : torch.device, optional (default torch.device("cuda"))
        Torch device for Cellpose and GPU cache cleanup.
    settings : dict, optional
        Detection settings dictionary passed to Cellpose.
    inner_neuropil_radius : int, optional (default 2)
        Radius in pixels of the inner exclusion zone for neuropil masks,
        used in intensity-ratio fallback.
    min_neuropil_pixels : int, optional (default 350)
        Minimum number of neuropil pixels, used in intensity-ratio fallback.

    Returns
    -------
    masks : numpy.ndarray or None
        Cellpose-detected masks in the red channel if Cellpose was used,
        otherwise None.
    redstats : numpy.ndarray
        Array of shape (n_cells, 2) where column 0 is the binary red cell
        label and column 1 is the red probability.
    """
    mimg = meanImg.copy()
    mimg2 = meanImg_chan2.copy()

    redstats = None
    masks = None
    if cellpose_chan2:
        try:
            logger.info(">>>> CELLPOSE estimating masks in anatomical channel")
            redstats, masks = cellpose_overlap(stats, mimg2, diameter=diameter,
                                               chan2_threshold=chan2_threshold,
                                               device=device, settings=settings)
        except:
            logger.info(
                "ERROR importing or running cellpose, continuing with intensity-based anatomical estimates"
            )

    if redstats is None:
        # subtract bleedthrough of green into red channel
        # non-rigid regression with nblks x nblks pieces
        nblks = 3
        #Ly, Lx = settings["Ly"], settings["Lx"]
        #mimg2_corr = correct_bleedthrough(Ly, Lx, nblks, mimg, mimg2)
        redstats = intensity_ratio(mimg2, stats, chan2_threshold=chan2_threshold,
                                  inner_neuropil_radius=inner_neuropil_radius,
                                  min_neuropil_pixels=min_neuropil_pixels)

    return masks, redstats

intensity_ratio #

intensity_ratio(mimg2, stats, chan2_threshold=0.65, inner_neuropil_radius=2, min_neuropil_pixels=350)

Classify cells as red using the intensity ratio of cell to surround.

Computes the ratio of channel 2 intensity inside each cell mask to the intensity in the surrounding neuropil, excluding other cells. Cells with a ratio above the threshold are classified as red.

Parameters:

Name Type Description Default
mimg2 ndarray

Red channel mean image of shape (Ly, Lx).

 required
stats ndarray

Array of ROI statistics dictionaries, each containing "ypix", "xpix", and "lam".

 required
chan2_threshold float, optional (default 0.65)

Threshold on the intensity ratio for red cell classification.

 0.65
inner_neuropil_radius int, optional (default 2)

Radius in pixels of the inner exclusion zone around each cell for neuropil mask creation.

 2
min_neuropil_pixels int, optional (default 350)

Minimum number of pixels in each neuropil mask.

 350

Returns:

Name Type Description
redcell ndarray

Array of shape (n_cells, 2) where column 0 is the binary red cell label and column 1 is the red probability (intensity ratio).
Source code in suite2p/detection/chan2detect.py 
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
def intensity_ratio(mimg2, stats, chan2_threshold=0.65, inner_neuropil_radius=2,
                    min_neuropil_pixels=350):
    """
    Classify cells as red using the intensity ratio of cell to surround.

    Computes the ratio of channel 2 intensity inside each cell mask to the
    intensity in the surrounding neuropil, excluding other cells. Cells with
    a ratio above the threshold are classified as red.

    Parameters
    ----------
    mimg2 : numpy.ndarray
        Red channel mean image of shape (Ly, Lx).
    stats : numpy.ndarray
        Array of ROI statistics dictionaries, each containing "ypix", "xpix",
        and "lam".
    chan2_threshold : float, optional (default 0.65)
        Threshold on the intensity ratio for red cell classification.
    inner_neuropil_radius : int, optional (default 2)
        Radius in pixels of the inner exclusion zone around each cell for
        neuropil mask creation.
    min_neuropil_pixels : int, optional (default 350)
        Minimum number of pixels in each neuropil mask.

    Returns
    -------
    redcell : numpy.ndarray
        Array of shape (n_cells, 2) where column 0 is the binary red cell
        label and column 1 is the red probability (intensity ratio).
    """
    Ly, Lx = mimg2.shape
    cell_pix = masks.create_cell_pix(stats, Ly=Ly, Lx=Lx)
    cell_masks0 = [
        masks.create_cell_mask(stat, Ly=Ly, Lx=Lx, allow_overlap=True) for stat in stats
    ]
    neuropil_ipix = 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=inner_neuropil_radius,
        min_neuropil_pixels=min_neuropil_pixels
    )
    cell_masks = np.zeros((len(stats), Ly * Lx), np.float32)
    neuropil_masks = np.zeros((len(stats), Ly * Lx), np.float32)
    for cell_mask, cell_mask0, neuropil_mask, neuropil_mask0 in zip(
            cell_masks, cell_masks0, neuropil_masks, neuropil_ipix):
        cell_mask[cell_mask0[0]] = cell_mask0[1]
        neuropil_mask[neuropil_mask0.astype(np.int64)] = 1. / len(neuropil_mask0)

    inpix = cell_masks @ mimg2.flatten()
    extpix = neuropil_masks @ mimg2.flatten()
    inpix = np.maximum(1e-3, inpix)
    redprob = inpix / (inpix + extpix)
    redcell = redprob > chan2_threshold

    return np.stack((redcell, redprob), axis=-1)

quadrant_mask #

quadrant_mask(Ly, Lx, ny, nx, sT)

Create a smoothed binary mask for a rectangular block of the image.

Sets a rectangular region to 1 and applies Gaussian smoothing to create soft edges for blending in the bleedthrough correction.

Parameters:

Name Type Description Default
Ly int

Height of the image in pixels.

 required
Lx int

Width of the image in pixels.

 required
ny ndarray

Y-indices defining the block rows.

 required
nx ndarray

X-indices defining the block columns.

 required
sT float

Standard deviation for Gaussian smoothing of the mask.

 required

Returns:

Name Type Description
mask ndarray

Smoothed mask of shape (Ly, Lx), dtype float32.
Source code in suite2p/detection/chan2detect.py 
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
def quadrant_mask(Ly, Lx, ny, nx, sT):
    """
    Create a smoothed binary mask for a rectangular block of the image.

    Sets a rectangular region to 1 and applies Gaussian smoothing to create
    soft edges for blending in the bleedthrough correction.

    Parameters
    ----------
    Ly : int
        Height of the image in pixels.
    Lx : int
        Width of the image in pixels.
    ny : numpy.ndarray
        Y-indices defining the block rows.
    nx : numpy.ndarray
        X-indices defining the block columns.
    sT : float
        Standard deviation for Gaussian smoothing of the mask.

    Returns
    -------
    mask : numpy.ndarray
        Smoothed mask of shape (Ly, Lx), dtype float32.
    """
    mask = np.zeros((Ly, Lx), np.float32)
    mask[np.ix_(ny, nx)] = 1
    mask = gaussian_filter(mask, sT)
    return mask

Copyright © 2023 Howard Hughes Medical Institute, Authored by Carsen Stringer and Marius Pachitariu.

pca_denoise #

pca_denoise(mov, block_size, n_comps_frac)

Denoise a movie using block-wise PCA reconstruction.

Splits the movie into spatial blocks, projects each block onto its top PCA components, reconstructs, and blends the blocks with spatial tapering.

Parameters:

Name Type Description Default
mov ndarray

Movie of shape (nframes, Ly, Lx). Modified in-place (mean subtracted then restored).

 required
block_size list of int

Block size [Lyb, Lxb] for spatial tiling.

 required
n_comps_frac float

Fraction of the smaller block dimension used to set the number of PCA components (number of PCs n_comps = min(Lyb, Lxb) * n_comps_frac).

 required

Returns:

Name Type Description
reconstruction ndarray

Denoised movie of shape (nframes, Ly, Lx).
Source code in suite2p/detection/denoise.py 
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
def pca_denoise(mov, block_size, n_comps_frac):
    """
    Denoise a movie using block-wise PCA reconstruction.

    Splits the movie into spatial blocks, projects each block onto its
    top PCA components, reconstructs, and blends the blocks with spatial
    tapering.

    Parameters
    ----------
    mov : numpy.ndarray
        Movie of shape (nframes, Ly, Lx). Modified in-place (mean
        subtracted then restored).
    block_size : list of int
        Block size [Lyb, Lxb] for spatial tiling.
    n_comps_frac : float
        Fraction of the smaller block dimension used to set the number
        of PCA components (number of PCs n_comps = min(Lyb, Lxb) * n_comps_frac).

    Returns
    -------
    reconstruction : numpy.ndarray
        Denoised movie of shape (nframes, Ly, Lx).
    """
    t0 = time.time()
    nframes, Ly, Lx = mov.shape
    yblock, xblock, nb, block_size = make_blocks(Ly, Lx, block_size=block_size)[:4]

    mov_mean = mov.mean(axis=0)
    mov -= mov_mean

    nblocks = len(yblock)
    Lyb, Lxb = block_size
    n_comps = int(min(min(Lyb * Lxb, nframes), min(Lyb, Lxb) * n_comps_frac))
    maskMul = spatial_taper(Lyb // 4, Lyb, Lxb).numpy()
    norm = np.zeros((Ly, Lx), np.float32)
    reconstruction = np.zeros_like(mov)
    block_re = np.zeros((nblocks, nframes, Lyb * Lxb))
    for i in range(nblocks):
        block = mov[:, yblock[i][0]:yblock[i][-1],
                    xblock[i][0]:xblock[i][-1]].reshape(-1, Lyb * Lxb)
        model = PCA(n_components=n_comps, random_state=0).fit(block)
        block_re[i] = (block @ model.components_.T) @ model.components_
        norm[yblock[i][0]:yblock[i][-1], xblock[i][0]:xblock[i][-1]] += maskMul

    block_re = block_re.reshape(nblocks, nframes, Lyb, Lxb)
    block_re *= maskMul
    for i in range(nblocks):
        reconstruction[:, yblock[i][0]:yblock[i][-1],
                       xblock[i][0]:xblock[i][-1]] += block_re[i]
    reconstruction /= norm
    logger.info("Binned movie denoised (for cell detection only) in %0.2f sec." %
          (time.time() - t0))
    reconstruction += mov_mean
    return reconstruction

Copyright © 2023 Howard Hughes Medical Institute, Authored by Carsen Stringer and Marius Pachitariu.

bin_movie #

bin_movie(f_reg, bin_size, yrange=None, xrange=None, badframes=None, nbins=5000)

Temporally bin the registered movie.

Reads batches of frames, excludes bad frames, crops to the valid region, and averages groups of bin_size frames to produce a downsampled movie.

Parameters:

Name Type Description Default
f_reg ndarray or BinaryFile

Registered movie of shape (n_frames, Ly, Lx).

 required
bin_size int

Number of frames to average per bin.

 required
yrange list of int

Two-element list [y_start, y_end] defining the Y crop range.

 None
xrange list of int

Two-element list [x_start, x_end] defining the X crop range.

 None
badframes ndarray

Boolean array of shape (n_frames,) where True marks frames to exclude.

 None
nbins int, optional (default 5000)

Maximum number of output binned frames.

 5000

Returns:

Name Type Description
mov ndarray

Binned movie of shape (num_binned_frames, Lyc, Lxc), dtype float32.
Source code in suite2p/detection/detect.py 
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
def bin_movie(f_reg, bin_size, yrange=None, xrange=None, badframes=None, nbins=5000):
    """
    Temporally bin the registered movie.

    Reads batches of frames, excludes bad frames, crops to the valid
    region, and averages groups of ``bin_size`` frames to produce a
    downsampled movie.

    Parameters
    ----------
    f_reg : numpy.ndarray or BinaryFile
        Registered movie of shape (n_frames, Ly, Lx).
    bin_size : int
        Number of frames to average per bin.
    yrange : list of int, optional
        Two-element list [y_start, y_end] defining the Y crop range.
    xrange : list of int, optional
        Two-element list [x_start, x_end] defining the X crop range.
    badframes : numpy.ndarray, optional
        Boolean array of shape (n_frames,) where True marks frames to exclude.
    nbins : int, optional (default 5000)
        Maximum number of output binned frames.

    Returns
    -------
    mov : numpy.ndarray
        Binned movie of shape (num_binned_frames, Lyc, Lxc), dtype float32.
    """
    n_frames = f_reg.shape[0]
    good_frames = ~badframes if badframes is not None else np.ones(n_frames, dtype=bool)
    batch_size = min(good_frames.sum(), 500)
    Lyc = yrange[1] - yrange[0]
    Lxc = xrange[1] - xrange[0]

    # Number of binned frames is rounded down when binning frames
    num_binned_frames = min(nbins, n_frames // bin_size)
    mov = np.zeros((num_binned_frames, Lyc, Lxc), np.float32)
    curr_bin_number = 0
    t0 = time.time()

    # Iterate over n_frames to maintain binning over TIME
    tstarts = np.arange(0, n_frames, batch_size)
    n_batches = min(nbins // (batch_size // bin_size), len(tstarts))
    tstarts = tstarts[np.linspace(0, len(tstarts) - 1, n_batches, dtype="int")]

    tqdm_out = TqdmToLogger(logger, level=logging.INFO)
    for n in trange(n_batches, mininterval=10, file=tqdm_out):
        tstart = tstarts[n]
        tend = min(tstart + batch_size, n_frames)
        data = f_reg[tstart : tend]

        # exclude badframes
        good_indices = good_frames[tstart : tend]
        if good_indices.mean() > 0.5:
            data = data[good_indices]

        # crop to valid region
        if yrange is not None and xrange is not None:
            data = data[:, slice(*yrange), slice(*xrange)]

        # bin in time
        if data.shape[0] > bin_size:
            # Downsample by binning via reshaping and taking mean of each bin
            # only if current batch size exceeds or matches bin_size
            n_d = data.shape[0]
            data = data[:(n_d // bin_size) * bin_size]
            data = data.reshape(-1, bin_size, Lyc, Lxc).astype(np.float32).mean(axis=1)
        else:
            # Current batch size is below bin_size (could have many bad frames in this batch)
            # Downsample taking the mean of batch to get a single bin
            data = data.mean(axis=0)[np.newaxis, :, :]
        # Only fill in binned data if not exceeding the number of bins mov has
        if mov.shape[0] > curr_bin_number:
            # Fill in binned data
            nb = data.shape[0]
            mov[curr_bin_number:curr_bin_number + nb] = data
            curr_bin_number += nb
    mov = mov[:curr_bin_number]
    logger.info("Binned movie of size [%d,%d,%d] created in %0.2f sec." %
          (mov.shape[0], mov.shape[1], mov.shape[2], time.time() - t0))
    return mov

detection_wrapper #

detection_wrapper(f_reg, diameter=[12.0, 12.0], tau=1.0, fs=30, meanImg_chan2=None, yrange=None, xrange=None, badframes=None, mov=None, preclassify=0.0, classifier_path=None, settings=default_settings()['detection'], device=torch.device('cuda'))

Run the full ROI detection pipeline on a registered movie.

Bins the movie in time, optionally denoises and high-pass filters, detects ROIs using the selected algorithm (sparsery, sourcery, or cellpose), computes ROI statistics, and optionally preclassifies and detects red cells in a second channel.

Parameters:

Name Type Description Default
f_reg ndarray or BinaryFile

Registered movie of shape (n_frames, Ly, Lx).

 required
diameter list of float, optional (default [12., 12.])

Expected cell diameter [dy, dx] in pixels.

 [12.0, 12.0]
tau float, optional (default 1.)

Timescale of the indicator in seconds, used to set bin size.

 1.0
fs float, optional (default 30)

Sampling rate in Hz, used with tau to set bin size.

 30
meanImg_chan2 ndarray

Mean image of the second channel, shape (Ly, Lx). If provided, red cell detection is performed.

 None
yrange list of int

Two-element list [y_start, y_end] defining the Y crop range.

 None
xrange list of int

Two-element list [x_start, x_end] defining the X crop range.

 None
badframes ndarray

Boolean array of shape (n_frames,) marking frames to exclude.

 None
mov ndarray

Pre-binned movie of shape (nbinned, Lyc, Lxc). If provided, skips the binning step.

 None
preclassify float, optional (default 0.)

If positive, apply a classifier and remove ROIs with probability below this threshold before final statistics.

 0.0
classifier_path str

Path to a saved classifier file. If None and preclassify > 0, uses the default user classifier.

 None
settings dict

Detection settings dictionary.

 default_settings()['detection']
device torch.device, optional (default torch.device("cuda"))

Torch device for cellpose-based detection.

 device('cuda')

Returns:

Name Type Description
new_settings dict

Dictionary with detection metadata including "meanImg_crop", "max_proj", "diameter", and algorithm-specific keys.
stat ndarray

Array of ROI statistics dictionaries, each containing "ypix", "xpix", "lam", "med", and other computed statistics.
redcell ndarray or None

Array of shape (n_cells, 2) with red cell labels and probabilities, or None if no second channel was provided.
Source code in suite2p/detection/detect.py 
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
def detection_wrapper(f_reg, diameter=[12., 12.], tau=1., fs=30, meanImg_chan2=None,
                      yrange=None, xrange=None, badframes=None, mov=None, 
                      preclassify=0., classifier_path=None, 
                      settings=default_settings()["detection"],
                      device=torch.device("cuda")):
    """
    Run the full ROI detection pipeline on a registered movie.

    Bins the movie in time, optionally denoises and high-pass filters,
    detects ROIs using the selected algorithm (sparsery, sourcery, or
    cellpose), computes ROI statistics, and optionally preclassifies and
    detects red cells in a second channel.

    Parameters
    ----------
    f_reg : numpy.ndarray or BinaryFile
        Registered movie of shape (n_frames, Ly, Lx).
    diameter : list of float, optional (default [12., 12.])
        Expected cell diameter [dy, dx] in pixels.
    tau : float, optional (default 1.)
        Timescale of the indicator in seconds, used to set bin size.
    fs : float, optional (default 30)
        Sampling rate in Hz, used with tau to set bin size.
    meanImg_chan2 : numpy.ndarray, optional
        Mean image of the second channel, shape (Ly, Lx). If provided,
        red cell detection is performed.
    yrange : list of int, optional
        Two-element list [y_start, y_end] defining the Y crop range.
    xrange : list of int, optional
        Two-element list [x_start, x_end] defining the X crop range.
    badframes : numpy.ndarray, optional
        Boolean array of shape (n_frames,) marking frames to exclude.
    mov : numpy.ndarray, optional
        Pre-binned movie of shape (nbinned, Lyc, Lxc). If provided,
        skips the binning step.
    preclassify : float, optional (default 0.)
        If positive, apply a classifier and remove ROIs with probability
        below this threshold before final statistics.
    classifier_path : str, optional
        Path to a saved classifier file. If None and preclassify > 0,
        uses the default user classifier.
    settings : dict, optional
        Detection settings dictionary.
    device : torch.device, optional (default torch.device("cuda"))
        Torch device for cellpose-based detection.

    Returns
    -------
    new_settings : dict
        Dictionary with detection metadata including "meanImg_crop",
        "max_proj", "diameter", and algorithm-specific keys.
    stat : numpy.ndarray
        Array of ROI statistics dictionaries, each containing "ypix",
        "xpix", "lam", "med", and other computed statistics.
    redcell : numpy.ndarray or None
        Array of shape (n_cells, 2) with red cell labels and probabilities,
        or None if no second channel was provided.
    """

    n_frames, Ly, Lx = f_reg.shape
    yrange = [0, Ly] if yrange is None else yrange
    xrange = [0, Lx] if xrange is None else xrange

    if mov is None:
        nbins = settings["nbins"]
        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)
    else:
        if mov.shape[1] != yrange[-1] - yrange[0]:
            raise ValueError("mov.shape[1] is not same size as yrange")
        elif mov.shape[2] != xrange[-1] - xrange[0]:
            raise ValueError("mov.shape[2] is not same size as xrange")

    if settings.get("inverted_activity", False):
        mov -= mov.min()
        mov *= -1
        mov -= mov.min()

    if settings.get("denoise", False):
        mov = pca_denoise(
            mov, block_size=settings["block_size"],
            n_comps_frac=0.5)

    meanImg = mov.mean(axis=0) 

    mov = utils.temporal_high_pass_filter(mov=mov, width=settings["highpass_time"])
    max_proj = mov.max(axis=0) 

    t0 = time.time()
    if settings["algorithm"] == "cellpose":
        if anatomical.CELLPOSE_INSTALLED:
            logger.info(f">>>> CELLPOSE finding masks in {settings['cellpose_settings']['img']}")
            new_settings, stat = anatomical.select_rois(meanImg, max_proj, settings=settings["cellpose_settings"],
                                          diameter=diameter,
                                          device=device)
        else:
            logger.info("Warning: cellpose did not import ", anatomical.cellpose_error)
            logger.info("cannot use anatomical mode, will use functional detection instead")

    if settings["algorithm"] != "cellpose" or not anatomical.CELLPOSE_INSTALLED:
        settings["algorithm"] = "sparsery" if settings["algorithm"] == "cellpose" else settings["algorithm"]
        sdmov = utils.standard_deviation_over_time(mov, batch_size=1000)
        if settings["algorithm"] == "sparsery":
            new_settings, stat = sparsedetect.sparsery(
                mov=mov, sdmov=sdmov,
                threshold_scaling=settings["threshold_scaling"],
                **settings["sparsery_settings"]
            )
        else:
            new_settings, stat = sourcery.sourcery(mov=mov, sdmov=sdmov, diameter=diameter,
                                              threshold_scaling=settings["threshold_scaling"],
                                              **settings["sourcery_settings"])
    logger.info("Detected %d ROIs, %0.2f sec" % (len(stat), time.time() - t0))
    stat = np.array(stat)

    if len(stat) == 0:
        raise ValueError(
            "no ROIs were found -- check registered binary and maybe try changing spatial scale / diameter / threshold_scaling"
        )

    # move ROIs to original coordinates
    ymin, xmin = int(yrange[0]), int(xrange[0])
    for s in stat:
        s["ypix"] += ymin; s["xpix"] += xmin; 
        if "med" in s:
            s["med"][0] += ymin; s["med"][1] += xmin

    if preclassify > 0.:
        if classifier_path is None:
            logger.info(f"NOTE: Applying user classifier at {str(user_classfile)}")
            classifier_path = user_classfile

        n0 = len(stat)
        stat = roi_stats(stat, Ly, Lx, diameter=diameter, 
                        max_overlap=None,
                        do_soma_crop=settings.get("soma_crop", 1), 
                        npix_norm_min=settings.get("npix_norm_min", 0.),
                        npix_norm_max=settings.get("npix_norm_max", 100.),
                        median=settings["algorithm"]=="cellpose")
        #import pdb; pdb.set_trace()
        iscell = classify(stat=stat, classfile=classifier_path)
        ic = (iscell[:, 1] > preclassify).flatten().astype("bool")
        stat = stat[ic]

        if len(stat) == 0:
            raise ValueError("no ROIs passed preclassify -- turn off preclassify or lower threshold_scaling")

        logger.info(f"preclassify threshold {preclassify}, {(~ic).sum()} ROIs removed")

    stat = roi_stats(stat, Ly, Lx, diameter=diameter, 
                        max_overlap=settings.get("max_overlap", 0.75),
                        do_soma_crop=settings.get("soma_crop", 1), 
                        npix_norm_min=settings.get("npix_norm_min", 0.),
                        npix_norm_max=settings.get("npix_norm_max", 100.),
                        median=settings["algorithm"]=="cellpose")
    logger.info("After removing by overlaps and npix_norm, %d ROIs remain" % (len(stat)))

    # if second channel, detect bright cells in second channel
    if meanImg_chan2 is not None:
        extraction_defaults = default_settings()["extraction"]
        redmasks, redcell = chan2detect.detect(meanImg, meanImg_chan2, stat, diameter=diameter,
                                    cellpose_chan2=settings.get("cellpose_chan2", False),
                                    chan2_threshold=settings.get("chan2_threshold", 0.65),
                                    settings=settings['cellpose_settings'],
                                    inner_neuropil_radius=extraction_defaults["inner_neuropil_radius"],
                                    min_neuropil_pixels=extraction_defaults["min_neuropil_pixels"],
                                    device=device)
        new_settings["chan2_masks"] = redmasks
    else:
        redcell = None

    new_settings["meanImg_crop"] = meanImg
    new_settings["max_proj"] = max_proj
    new_settings["diameter"] = diameter

    return new_settings, stat, redcell

Copyright © 2023 Howard Hughes Medical Institute, Authored by Carsen Stringer and Marius Pachitariu.

compute_gt_matches #

compute_gt_matches(img, masks, stat_func, ops=None, reg_file=None, threshold=0.5)

anatomical img and masks matched to functional ROIs in stat_func

Source code in suite2p/detection/metrics.py 
13
14
15
16
17
18
19
20
21
22
23
def compute_gt_matches(img, masks, stat_func, ops=None, reg_file=None, threshold=0.5):
    """ anatomical img and masks matched to functional ROIs in stat_func """
    Ly, Lx = masks.shape
    stat_anat, iorig = extend_anatomical(img, masks, ops=ops, reg_file=reg_file)
    iou, iout, preds, ap = match_func_anat(stat_func, stat_anat, Ly, Lx, threshold)

    chosen_cells = iout > threshold
    func_ids = preds[chosen_cells] - 1
    overlaps = iout[chosen_cells]

    return stat_anat, iorig, iou, func_ids, overlaps

match_func_anat #

match_func_anat(stat_func, stat_anat, Ly, Lx, threshold=0.5)

match functional ROIs to anatomical ROIs by correlation

Source code in suite2p/detection/metrics.py 
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
def match_func_anat(stat_func, stat_anat, Ly, Lx, threshold=0.5):
    """ match functional ROIs to anatomical ROIs by correlation"""
    iou = np.zeros((len(stat_anat), len(stat_func)))
    ly = 15
    for i, sf in enumerate(stat_func):
        if sf["ypix"].size < 20:
            continue
        ypix, xpix, lam = sf["ypix"].copy(), sf["xpix"].copy(), sf["lam"].copy()
        lam /= (lam**2).sum()**0.5
        # box around ROI
        ymed, xmed = sf["med"][0], sf["med"][1]
        inds = (slice(max(0, ymed - ly),
                      min(ymed + ly, Ly)), slice(max(0, xmed - ly), min(xmed + ly, Lx)))
        mf = np.zeros((Ly, Lx), np.float32)
        mf[ypix, xpix] = lam
        mfc = mf[inds].flatten()
        mfc /= (mfc**2).sum()**0.5

        # matched anatomical masks (will not compute IOU for all masks)
        for j, sa in enumerate(stat_anat):
            ypix_a, xpix_a = sa["ypix"], sa["xpix"]
            if (np.logical_and(ypix_a > inds[0].start, ypix_a < inds[0].stop).sum() > 0
                    and np.logical_and(xpix_a > inds[1].start, xpix_a
                                       < inds[1].stop).sum() > 0):
                lam_a = sa["lam"].copy()
                lam_a /= (lam_a**2).sum()**0.5
                ma = np.zeros((Ly, Lx), np.float32)
                ma[ypix_a, xpix_a] = lam_a
                mac = ma[inds].flatten()
                mac /= ((mac**2).sum()**0.5 + 1e-10)
                iou[j, i] = (mac * mfc).sum()
        if i % 1000 == 0:
            print("%d ROIs processed" % i)
    print("%d ROIs processed" % i)

    n_true = len(stat_anat)
    n_pred = len(stat_func)
    iout, preds = match_masks(iou)
    tp = (iout > threshold).sum()
    print((iout > threshold).sum())
    fn = n_true - tp
    fp = n_pred - tp
    ap = tp / (fn + tp + fp)
    print("TP: %d, FN: %d, FP: %d, AP: %0.3f" % (tp, fn, fp, ap))

    return iou, iout, preds, ap

Copyright © 2023 Howard Hughes Medical Institute, Authored by Carsen Stringer and Marius Pachitariu.

connected_region #

connected_region(stat, connected, Ly, Lx)

Optionally restrict each ROI to its largest connected component.

Parameters:

Name Type Description Default
stat list of dict

List of ROI statistics dictionaries, each containing "ypix", "xpix", and "lam".

 required
connected bool

If True, apply connected component extraction to each ROI.

 required
Ly int

Height of the image in pixels.

 required
Lx int

Width of the image in pixels.

 required

Returns:

Name Type Description
stat list of dict

Updated list of ROI dictionaries.
Source code in suite2p/detection/sourcery.py 
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
def connected_region(stat, connected, Ly, Lx):
    """
    Optionally restrict each ROI to its largest connected component.

    Parameters
    ----------
    stat : list of dict
        List of ROI statistics dictionaries, each containing "ypix", "xpix",
        and "lam".
    connected : bool
        If True, apply connected component extraction to each ROI.
    Ly : int
        Height of the image in pixels.
    Lx : int
        Width of the image in pixels.

    Returns
    -------
    stat : list of dict
        Updated list of ROI dictionaries.
    """
    if connected:
        for j in range(len(stat)):
            stat[j] = get_connected(Ly, Lx, stat[j])
    return stat

create_neuropil_basis #

create_neuropil_basis(diameter, Ly, Lx)

Compute Fourier-based neuropil basis functions for the image.

Creates a set of 2D Fourier basis functions tiled across the image, used to model spatially varying neuropil contamination.

Parameters:

Name Type Description Default
diameter float or list of float

Cell diameter used to determine the tiling density.

 required
Ly int

Height of the image in pixels.

 required
Lx int

Width of the image in pixels.

 required

Returns:

Name Type Description
S ndarray

Normalized neuropil basis functions of shape (Ly, Lx, nbasis).
Source code in suite2p/detection/sourcery.py 
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
def create_neuropil_basis(diameter, Ly, Lx):
    """
    Compute Fourier-based neuropil basis functions for the image.

    Creates a set of 2D Fourier basis functions tiled across the image,
    used to model spatially varying neuropil contamination.

    Parameters
    ----------
    diameter : float or list of float
        Cell diameter used to determine the tiling density.
    Ly : int
        Height of the image in pixels.
    Lx : int
        Width of the image in pixels.

    Returns
    -------
    S : numpy.ndarray
        Normalized neuropil basis functions of shape (Ly, Lx, nbasis).
    """

    ratio_neuropil = 6.
    tile_factor = 1.

    ntilesY = 1 + 2 * int(
        np.ceil(tile_factor * Ly / (ratio_neuropil * diameter[0] / 2)) / 2)
    ntilesX = 1 + 2 * int(
        np.ceil(tile_factor * Lx / (ratio_neuropil * diameter[1] / 2)) / 2)
    ntilesY = np.maximum(2, ntilesY)
    ntilesX = np.maximum(2, ntilesX)
    yc = np.linspace(1, Ly, ntilesY)
    xc = np.linspace(1, Lx, ntilesX)
    ys = np.arange(0, Ly)
    xs = np.arange(0, Lx)

    Kx = np.ones((Lx, ntilesX), "float32")
    Ky = np.ones((Ly, ntilesY), "float32")
    if 1:
        # basis functions are fourier modes
        for k in range(int((ntilesX - 1) / 2)):
            Kx[:, 2 * k + 1] = np.sin(2 * math.pi * (xs + 0.5) * (1 + k) / Lx)
            Kx[:, 2 * k + 2] = np.cos(2 * math.pi * (xs + 0.5) * (1 + k) / Lx)
        for k in range(int((ntilesY - 1) / 2)):
            Ky[:, 2 * k + 1] = np.sin(2 * math.pi * (ys + 0.5) * (1 + k) / Ly)
            Ky[:, 2 * k + 2] = np.cos(2 * math.pi * (ys + 0.5) * (1 + k) / Ly)
    else:
        for k in range(ntilesX):
            Kx[:, k] = np.cos(math.pi * (xs + 0.5) * k / Lx)
        for k in range(ntilesY):
            Ky[:, k] = np.cos(math.pi * (ys + 0.5) * k / Ly)

    S = np.zeros((ntilesY, ntilesX, Ly, Lx), np.float32)
    for kx in range(ntilesX):
        for ky in range(ntilesY):
            S[ky, kx, :, :] = np.outer(Ky[:, ky], Kx[:, kx])

    S = np.reshape(S, (ntilesY * ntilesX, Ly * Lx))
    S = S / np.reshape(np.sum(S**2, axis=-1)**0.5, (-1, 1))
    S = np.transpose(S, (1, 0)).copy()
    S = np.reshape(S, (Ly, Lx, -1))
    return S

extendROI #

extendROI(ypix, xpix, Ly, Lx, niter=1)

Expand an ROI by one pixel in each cardinal direction.

Adds the 4-connected neighbors of the current pixel set and removes any that fall outside the image boundaries. Repeated for niter iterations.

Parameters:

Name Type Description Default
ypix ndarray

Y-coordinates of the current ROI pixels.

 required
xpix ndarray

X-coordinates of the current ROI pixels.

 required
Ly int

Height of the image in pixels.

 required
Lx int

Width of the image in pixels.

 required
niter int, optional (default 1)

Number of expansion iterations.

 1

Returns:

Name Type Description
ypix ndarray

Expanded Y-coordinates.
xpix ndarray

Expanded X-coordinates.
Source code in suite2p/detection/sourcery.py 
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
def extendROI(ypix, xpix, Ly, Lx, niter=1):
    """
    Expand an ROI by one pixel in each cardinal direction.

    Adds the 4-connected neighbors of the current pixel set and removes
    any that fall outside the image boundaries. Repeated for `niter`
    iterations.

    Parameters
    ----------
    ypix : numpy.ndarray
        Y-coordinates of the current ROI pixels.
    xpix : numpy.ndarray
        X-coordinates of the current ROI pixels.
    Ly : int
        Height of the image in pixels.
    Lx : int
        Width of the image in pixels.
    niter : int, optional (default 1)
        Number of expansion iterations.

    Returns
    -------
    ypix : numpy.ndarray
        Expanded Y-coordinates.
    xpix : numpy.ndarray
        Expanded X-coordinates.
    """
    for k in range(niter):
        yx = ((ypix, ypix, ypix, ypix - 1, ypix + 1), (xpix, xpix + 1, xpix - 1, xpix,
                                                       xpix))
        yx = np.array(yx)
        yx = yx.reshape((2, -1))
        yu = np.unique(yx, axis=1)
        ix = np.all((yu[0] >= 0, yu[0] < Ly, yu[1] >= 0, yu[1] < Lx), axis=0)
        ypix, xpix = yu[:, ix]
    return ypix, xpix

getSVDdata #

getSVDdata(mov, diameter)

Compute SVD spatial components from temporally binned movie data.

Smooths each frame with a 2D Gaussian filter, computes the temporal covariance matrix, and returns the top SVD spatial components reshaped as images.

Parameters:

Name Type Description Default
mov ndarray

Temporally binned movie of shape (nbins, Ly, Lx).

 required
diameter float or list of float

Cell diameter used to set the Gaussian smoothing sigma (sigma = diameter / 10).

 required

Returns:

Name Type Description
U ndarray

Spatial SVD components of shape (Ly, Lx, nsvd).
u ndarray

Temporal SVD components of shape (nbins, nsvd).
Source code in suite2p/detection/sourcery.py 
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
def getSVDdata(mov: np.ndarray, diameter):
    """
    Compute SVD spatial components from temporally binned movie data.

    Smooths each frame with a 2D Gaussian filter, computes the temporal covariance
    matrix, and returns the top SVD spatial components reshaped as images.

    Parameters
    ----------
    mov : numpy.ndarray
        Temporally binned movie of shape (nbins, Ly, Lx).
    diameter : float or list of float
        Cell diameter used to set the Gaussian smoothing sigma (sigma = diameter / 10).

    Returns
    -------
    U : numpy.ndarray
        Spatial SVD components of shape (Ly, Lx, nsvd).
    u : numpy.ndarray
        Temporal SVD components of shape (nbins, nsvd).
    """
    nbins, Lyc, Lxc = np.shape(mov)

    sig = diameter / 10.  # PICK UP
    for j in range(nbins):
        mov[j, :, :] = gaussian_filter(mov[j, :, :], sig)

    # compute noise variance across frames
    mov = np.reshape(mov, (-1, Lyc * Lxc))
    # compute covariance of binned frames
    cov = mov @ mov.transpose() / mov.shape[1]
    cov = cov.astype("float32")

    nsvd_for_roi = min(nbins, int(cov.shape[0] / 2))
    u, s, v = np.linalg.svd(cov)

    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 U, u

getSVDproj #

getSVDproj(mov, u, diameter, smooth_masks=False)

Project temporally binned movie data onto existing SVD temporal components.

Optionally smooths frames before projection. Used during refinement to recompute spatial components from updated movie data.

Parameters:

Name Type Description Default
mov ndarray

Temporally binned movie of shape (nbins, Ly, Lx).

 required
u ndarray

Temporal SVD components of shape (nbins, nsvd).

 required
diameter float or list of float

Cell diameter used to set the Gaussian smoothing sigma.

 required
smooth_masks bool, optional (default False)

If True, smooth each frame before projection.

 False

Returns:

Name Type Description
U ndarray

Spatial SVD projections of shape (Ly, Lx, nsvd).
Source code in suite2p/detection/sourcery.py 
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
def getSVDproj(mov: np.ndarray, u, diameter, smooth_masks=False):
    """
    Project temporally binned movie data onto existing SVD temporal components.

    Optionally smooths frames before projection. Used during refinement to
    recompute spatial components from updated movie data.

    Parameters
    ----------
    mov : numpy.ndarray
        Temporally binned movie of shape (nbins, Ly, Lx).
    u : numpy.ndarray
        Temporal SVD components of shape (nbins, nsvd).
    diameter : float or list of float
        Cell diameter used to set the Gaussian smoothing sigma.
    smooth_masks : bool, optional (default False)
        If True, smooth each frame before projection.

    Returns
    -------
    U : numpy.ndarray
        Spatial SVD projections of shape (Ly, Lx, nsvd).
    """
    nbins, Lyc, Lxc = np.shape(mov)
    if smooth_masks:
        sig = np.maximum([.5, .5], diameter / 20.)
        for j in range(nbins):
            mov[j, :, :] = gaussian_filter(mov[j, :, :], sig)

    mov = np.reshape(mov, (-1, Lyc * Lxc))

    U = u.transpose() @ mov
    U = U.transpose().copy().reshape((Lyc, Lxc, -1))

    return U

getStU #

getStU(diameter, U)

Compute neuropil basis functions and their covariances with the SVD spatial components.

Parameters:

Name Type Description Default
diameter float or list of float

Cell diameter used for neuropil basis construction.

 required
U ndarray

Spatial SVD components of shape (Ly, Lx, nsvd).

 required

Returns:

Name Type Description
S ndarray

Neuropil basis functions of shape (Ly, Lx, nbasis).
StU ndarray

Cross-covariance of neuropil basis with SVD components, shape (nbasis, nsvd).
StS ndarray

Auto-covariance of neuropil basis, shape (nbasis, nbasis).
Source code in suite2p/detection/sourcery.py 
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
def getStU(diameter, U):
    """
    Compute neuropil basis functions and their covariances with the SVD spatial components.

    Parameters
    ----------
    diameter : float or list of float
        Cell diameter used for neuropil basis construction.
    U : numpy.ndarray
        Spatial SVD components of shape (Ly, Lx, nsvd).

    Returns
    -------
    S : numpy.ndarray
        Neuropil basis functions of shape (Ly, Lx, nbasis).
    StU : numpy.ndarray
        Cross-covariance of neuropil basis with SVD components, shape (nbasis, nsvd).
    StS : numpy.ndarray
        Auto-covariance of neuropil basis, shape (nbasis, nbasis).
    """
    Lyc, Lxc, nbins = np.shape(U)
    S = create_neuropil_basis(diameter, Lyc, Lxc)
    # compute covariance of neuropil masks with spatial masks
    StU = S.reshape((Lyc * Lxc, -1)).transpose() @ U.reshape((Lyc * Lxc, -1))
    StS = S.reshape((Lyc * Lxc, -1)).transpose() @ S.reshape((Lyc * Lxc, -1))
    #U = np.reshape(U, (-1,Lyc,Lxc))
    return S, StU, StS

getVmap #

getVmap(Ucell, sig)

Compute the variance ratio map from SVD spatial components.

Smooths the spatial components and computes the ratio of smoothed variance to total variance, producing a map that highlights regions with locally correlated activity.

Parameters:

Name Type Description Default
Ucell ndarray

Residual SVD spatial components of shape (Ly, Lx, nsvd).

 required
sig numpy.ndarray or list of float

Gaussian smoothing sigma [sy, sx] in pixels.

 required

Returns:

Name Type Description
log_variances ndarray

Variance ratio map of shape (Ly, Lx), dtype float64.
us ndarray

Smoothed spatial components of shape (Ly, Lx, nsvd).
Source code in suite2p/detection/sourcery.py 
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
def getVmap(Ucell, sig):
    """
    Compute the variance ratio map from SVD spatial components.

    Smooths the spatial components and computes the ratio of smoothed
    variance to total variance, producing a map that highlights regions
    with locally correlated activity.

    Parameters
    ----------
    Ucell : numpy.ndarray
        Residual SVD spatial components of shape (Ly, Lx, nsvd).
    sig : numpy.ndarray or list of float
        Gaussian smoothing sigma [sy, sx] in pixels.

    Returns
    -------
    log_variances : numpy.ndarray
        Variance ratio map of shape (Ly, Lx), dtype float64.
    us : numpy.ndarray
        Smoothed spatial components of shape (Ly, Lx, nsvd).
    """
    us = gaussian_filter(Ucell, [sig[0], sig[1], 0.], mode="wrap")
    # compute log variance at each location
    log_variances = (us**2).mean(axis=-1) / gaussian_filter(
        (Ucell**2).mean(axis=-1), sig, mode="wrap")
    return log_variances.astype("float64"), us

get_connected #

get_connected(Ly, Lx, stat)

Extract the connected component of an ROI starting from its brightest pixel.

Grows outward from the pixel with maximum weight, keeping only pixels that are contiguous and have non-zero weight.

Parameters:

Name Type Description Default
Ly int

Height of the image in pixels.

 required
Lx int

Width of the image in pixels.

 required
stat dict

ROI statistics dictionary containing "ypix", "xpix", and "lam". Modified in-place.

 required

Returns:

Name Type Description
stat dict

Updated ROI dictionary with connected pixels only.
Source code in suite2p/detection/sourcery.py 
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
def get_connected(Ly, Lx, stat):
    """
    Extract the connected component of an ROI starting from its brightest pixel.

    Grows outward from the pixel with maximum weight, keeping only pixels
    that are contiguous and have non-zero weight.

    Parameters
    ----------
    Ly : int
        Height of the image in pixels.
    Lx : int
        Width of the image in pixels.
    stat : dict
        ROI statistics dictionary containing "ypix", "xpix", and "lam".
        Modified in-place.

    Returns
    -------
    stat : dict
        Updated ROI dictionary with connected pixels only.
    """
    ypix, xpix, lam = stat["ypix"], stat["xpix"], stat["lam"]
    i0 = lam.argmax()
    mask = np.zeros((Ly, Lx))
    mask[ypix, xpix] = lam
    ypix, xpix = ypix[i0], xpix[i0]
    nsel = 1
    while 1:
        ypix, xpix = extendROI(ypix, xpix, Ly, Lx)
        ix = mask[ypix, xpix] > 1e-10
        ypix, xpix = ypix[ix], xpix[ix]
        if len(ypix) <= nsel:
            break
        nsel = len(ypix)
    lam = mask[ypix, xpix]
    stat["ypix"], stat["xpix"], stat["lam"] = ypix, xpix, lam
    return stat

iter_extend #

iter_extend(ypix, xpix, Ucell, code, refine=-1, change_codes=False)

Iteratively extend an ROI by projecting onto the SVD code vector.

Starting from seed pixels, repeatedly expands the region and keeps pixels whose projection onto the code vector exceeds a threshold. Stops when no new pixels are added or the region exceeds 10000 pixels.

Parameters:

Name Type Description Default
ypix ndarray or int

Initial Y-coordinates of the ROI seed.

 required
xpix ndarray or int

Initial X-coordinates of the ROI seed.

 required
Ucell ndarray

Residual SVD spatial components of shape (Ly, Lx, nsvd).

 required
code ndarray

SVD code vector for this ROI, shape (nsvd,).

 required
refine (int, optional(default - 1))

Refinement stage indicator. Negative means initial detection.

 -1
change_codes bool, optional (default False)

If True and refine < 0, update the code vector during extension.

 False

Returns:

Name Type Description
ypix ndarray

Y-coordinates of the extended ROI.
xpix ndarray

X-coordinates of the extended ROI.
lam ndarray

Normalized pixel weights (projections onto code vector).
ix ndarray

Boolean mask of pixels kept in the final iteration.
code ndarray

Updated code vector, shape (nsvd,).
Source code in suite2p/detection/sourcery.py 
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
def iter_extend(ypix, xpix, Ucell, code, refine=-1, change_codes=False):
    """
    Iteratively extend an ROI by projecting onto the SVD code vector.

    Starting from seed pixels, repeatedly expands the region and keeps
    pixels whose projection onto the code vector exceeds a threshold.
    Stops when no new pixels are added or the region exceeds 10000 pixels.

    Parameters
    ----------
    ypix : numpy.ndarray or int
        Initial Y-coordinates of the ROI seed.
    xpix : numpy.ndarray or int
        Initial X-coordinates of the ROI seed.
    Ucell : numpy.ndarray
        Residual SVD spatial components of shape (Ly, Lx, nsvd).
    code : numpy.ndarray
        SVD code vector for this ROI, shape (nsvd,).
    refine : int, optional (default -1)
        Refinement stage indicator. Negative means initial detection.
    change_codes : bool, optional (default False)
        If True and refine < 0, update the code vector during extension.

    Returns
    -------
    ypix : numpy.ndarray
        Y-coordinates of the extended ROI.
    xpix : numpy.ndarray
        X-coordinates of the extended ROI.
    lam : numpy.ndarray
        Normalized pixel weights (projections onto code vector).
    ix : numpy.ndarray
        Boolean mask of pixels kept in the final iteration.
    code : numpy.ndarray
        Updated code vector, shape (nsvd,).
    """
    Lyc, Lxc, nsvd = Ucell.shape
    npix = 0
    iter = 0
    while npix < 10000:
        npix = ypix.size
        ypix, xpix = extendROI(ypix, xpix, Lyc, Lxc, 1)
        usub = Ucell[ypix, xpix, :]
        lam = usub @ np.expand_dims(code, axis=1)
        lam = np.squeeze(lam, axis=1)
        # ix = lam>max(0, np.mean(lam)/3)
        ix = lam > max(0, lam.max() / 5.0)
        if ix.sum() == 0:
            break
        ypix, xpix, lam = ypix[ix], xpix[ix], lam[ix]
        lam = lam / np.sum(lam**2 + 1e-10)**.5
        if refine < 0 and change_codes:
            code = lam @ usub[ix, :]
        if iter == 0:
            sgn = 1.
            #sgn = np.sign(ix.sum()-npix)
        if np.sign(sgn * (ix.sum() - npix)) <= 0:
            break
        else:
            npix = ypix.size
        iter += 1
    return ypix, xpix, lam, ix, code

localMax #

localMax(V, footprint, thres)

Find local maxima of a correlation map above a threshold.

Uses a maximum filter with the given footprint to identify pixels that are local maxima and exceed the threshold.

Parameters:

Name Type Description Default
V ndarray

2D correlation map of shape (Ly, Lx).

 required
footprint ndarray

2D boolean footprint for the maximum filter (usually circular).

 required
thres float

Minimum value for a local maximum to be included.

 required

Returns:

Name Type Description
i ndarray

Y-indices of local maxima, dtype int32.
j ndarray

X-indices of local maxima, dtype int32.
Source code in suite2p/detection/sourcery.py 
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
def localMax(V, footprint, thres):
    """
    Find local maxima of a correlation map above a threshold.

    Uses a maximum filter with the given footprint to identify pixels that
    are local maxima and exceed the threshold.

    Parameters
    ----------
    V : numpy.ndarray
        2D correlation map of shape (Ly, Lx).
    footprint : numpy.ndarray
        2D boolean footprint for the maximum filter (usually circular).
    thres : float
        Minimum value for a local maximum to be included.

    Returns
    -------
    i : numpy.ndarray
        Y-indices of local maxima, dtype int32.
    j : numpy.ndarray
        X-indices of local maxima, dtype int32.
    """
    maxV = filters.maximum_filter(V, footprint=footprint, mode="reflect")
    imax = V > np.maximum(thres, maxV - 1e-10)
    i, j = imax.nonzero()
    i = i.astype(np.int32)
    j = j.astype(np.int32)
    return i, j

morphOpen #

morphOpen(V, footprint)

Compute the morphological opening of a correlation map.

Applies a minimum filter followed by a maximum filter (negated minimum of negation) using the given footprint to remove small bright features.

Parameters:

Name Type Description Default
V ndarray

2D correlation map of shape (Ly, Lx).

 required
footprint ndarray

2D boolean footprint for the morphological operation.

 required

Returns:

Name Type Description
vrem ndarray

Morphologically opened image of shape (Ly, Lx).
Source code in suite2p/detection/sourcery.py 
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
def morphOpen(V, footprint):
    """
    Compute the morphological opening of a correlation map.

    Applies a minimum filter followed by a maximum filter (negated minimum
    of negation) using the given footprint to remove small bright features.

    Parameters
    ----------
    V : numpy.ndarray
        2D correlation map of shape (Ly, Lx).
    footprint : numpy.ndarray
        2D boolean footprint for the morphological operation.

    Returns
    -------
    vrem : numpy.ndarray
        Morphologically opened image of shape (Ly, Lx).
    """
    vrem = filters.minimum_filter(V, footprint=footprint)
    vrem = -filters.minimum_filter(-vrem, footprint=footprint)
    return vrem

sourcery #

sourcery(mov, sdmov, diameter, threshold_scaling=1.0, connected=True, max_iterations=20, smooth_masks=False)

Detect ROIs using the Sourcery algorithm (SVD-based iterative detection).

Computes SVD components of the movie, builds neuropil basis functions and subtracts them, and iteratively detects ROIs by finding local maxima in the variance ratio map and extending them via similarity in SVD projection.

Parameters:

Name Type Description Default
mov ndarray

Temporally binned movie of shape (nbins, Ly, Lx). Divided by sdmov in-place before processing.

 required
sdmov ndarray

Standard deviation of the movie, shape (Ly * Lx,) or (Ly, Lx), used for normalization.

 required
diameter float or list of float

Expected cell diameter in pixels.

 required
threshold_scaling float, optional (default 1.0)

Scaling factor for the peak detection threshold.

 1.0
connected bool, optional (default True)

If True, restrict each detected ROI to its largest connected component.

 True
max_iterations int, optional (default 20)

Maximum number of detection and refinement iterations.

 20
smooth_masks bool, optional (default False)

If True, spatially smooth frames before SVD projection during refinement.

 False

Returns:

Name Type Description
new_settings dict

Dictionary with detection metadata including "diameter", "Vcorr", and placeholder keys "Vmax", "ihop", "Vsplit", "Vmap", "spatscale_pix".
stat ndarray

Array of ROI statistics dictionaries, each containing "ypix", "xpix", and "lam".
Source code in suite2p/detection/sourcery.py 
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
def sourcery(mov: np.ndarray, sdmov, diameter, threshold_scaling=1.0,
             connected=True, max_iterations=20, smooth_masks=False):
    """
    Detect ROIs using the Sourcery algorithm (SVD-based iterative detection).

    Computes SVD components of the movie, builds neuropil basis functions 
    and subtracts them, and iteratively detects ROIs by finding local 
    maxima in the variance ratio map and extending them via similarity 
    in SVD projection.

    Parameters
    ----------
    mov : numpy.ndarray
        Temporally binned movie of shape (nbins, Ly, Lx). Divided by
        `sdmov` in-place before processing.
    sdmov : numpy.ndarray
        Standard deviation of the movie, shape (Ly * Lx,) or (Ly, Lx),
        used for normalization.
    diameter : float or list of float
        Expected cell diameter in pixels.
    threshold_scaling : float, optional (default 1.0)
        Scaling factor for the peak detection threshold.
    connected : bool, optional (default True)
        If True, restrict each detected ROI to its largest connected component.
    max_iterations : int, optional (default 20)
        Maximum number of detection and refinement iterations.
    smooth_masks : bool, optional (default False)
        If True, spatially smooth frames before SVD projection during refinement.

    Returns
    -------
    new_settings : dict
        Dictionary with detection metadata including "diameter", "Vcorr",
        and placeholder keys "Vmax", "ihop", "Vsplit", "Vmap",
        "spatscale_pix".
    stat : numpy.ndarray
        Array of ROI statistics dictionaries, each containing "ypix",
        "xpix", and "lam".
    """
    mov /= sdmov

    change_codes = True
    t0 = time.time()

    U, u = getSVDdata(mov=mov, diameter=diameter)  # get SVD components
    S, StU, StS = getStU(diameter, U)
    Ly, Lx, nsvd = U.shape
    d0 = diameter
    sig = np.ceil(d0 / 4)  # smoothing constant
    # make array of radii values of size (2*d0+1,2*d0+1)
    rs, dy, dx = circleMask(d0)
    nsvd = U.shape[-1]
    nbasis = S.shape[-1]
    codes = np.zeros((0, nsvd), np.float32)
    LtU = np.zeros((0, nsvd), np.float32)
    LtS = np.zeros((0, nbasis), np.float32)
    L = np.zeros((Ly, Lx, 0), np.float32)
    # regress maps onto basis functions and subtract neuropil contribution
    neu = np.linalg.solve(StS, StU).astype("float32")
    Ucell = U - (S.reshape((-1, nbasis)) @ neu).reshape(U.shape)

    it = 0
    ncells = 0
    refine = -1

    # initialize
    ypix, xpix, lam = [], [], []
    logger.info(f"max_iterations = {max_iterations}; will stop when no more peaks above threshold, or max_iterations reached")
    for it in range(max_iterations):
        if refine < 0:
            V, us = getVmap(Ucell, sig)
            if it == 0:
                vrem = morphOpen(V, rs <= 1.)
            V = V - vrem  # make V more uniform
            if it == 0:
                V = V.astype("float64")
                # find indices of all maxima in +/- 1 range
                maxV = filters.maximum_filter(V, footprint=np.ones((3, 3)),
                                              mode="reflect")
                imax = V > (maxV - 1e-10)
                peaks = V[imax]
                # use the median of these peaks to decide if ROI is accepted
                thres = threshold_scaling * np.median(peaks[peaks > 1e-4])
                Vcorr = V.copy()
            V = np.minimum(V, Vcorr)

            # add extra ROIs here
            n = ncells
            while n < ncells + 200:
                ind = np.argmax(V)
                i, j = np.unravel_index(ind, V.shape)
                if V[i, j] < thres:
                    break
                yp, xp, la, ix, code = iter_extend(i, j, Ucell, us[i, j, :],
                                                   change_codes=change_codes)
                codes = np.append(codes, np.expand_dims(code, axis=0), axis=0)
                ypix.append(yp)
                xpix.append(xp)
                lam.append(la)
                Ucell[ypix[n], xpix[n], :] -= np.outer(lam[n], codes[n, :])

                yp, xp = extendROI(yp, xp, Ly, Lx, int(np.mean(d0)))
                V[yp, xp] = 0
                n += 1
            newcells = len(ypix) - ncells
            if it == 0:
                Nfirst = newcells
            L = np.append(L, np.zeros((Ly, Lx, newcells), "float32"), axis=-1)
            LtU = np.append(LtU, np.zeros((newcells, nsvd), "float32"), axis=0)
            LtS = np.append(LtS, np.zeros((newcells, nbasis), "float32"), axis=0)
            for n in range(ncells, len(ypix)):
                L[ypix[n], xpix[n], n] = lam[n]
                LtU[n, :] = lam[n] @ U[ypix[n], xpix[n], :]
                LtS[n, :] = lam[n] @ S[ypix[n], xpix[n], :]
            ncells += newcells

            # regression with neuropil
            LtL = L.reshape((-1, ncells)).transpose() @ L.reshape((-1, ncells))
            cellcode = np.concatenate((LtL, LtS), axis=1)
            neucode = np.concatenate((LtS.transpose(), StS), axis=1)
            codes = np.concatenate((cellcode, neucode), axis=0)
            Ucode = np.concatenate((LtU, StU), axis=0)
            codes = np.linalg.solve(codes + 1e-3 * np.eye((codes.shape[0])),
                                    Ucode).astype("float32")
            neu = codes[ncells:, :]
            codes = codes[:ncells, :]

        Ucell = U - (S.reshape((-1, nbasis)) @ neu + L.reshape(
            (-1, ncells)) @ codes).reshape(U.shape)
        # reestimate masks
        n, k = 0, 0
        while n < len(ypix):
            Ucell[ypix[n], xpix[n], :] += np.outer(lam[n], codes[k, :])
            ypix[n], xpix[n], lam[n], ix, codes[n, :] = iter_extend(
                ypix[n], xpix[n], Ucell, codes[k, :], refine, change_codes=change_codes)
            k += 1
            if ix.sum() == 0:
                logger.info("dropped ROI with no pixels")
                del ypix[n], xpix[n], lam[n]
                continue
            Ucell[ypix[n], xpix[n], :] -= np.outer(lam[n], codes[n, :])
            n += 1
        codes = codes[:n, :]
        ncells = len(ypix)
        L = np.zeros((Ly, Lx, ncells), "float32")
        LtU = np.zeros((ncells, nsvd), "float32")
        LtS = np.zeros((ncells, nbasis), "float32")
        for n in range(ncells):
            L[ypix[n], xpix[n], n] = lam[n]
            if refine < 0:
                LtU[n, :] = lam[n] @ U[ypix[n], xpix[n], :]
                LtS[n, :] = lam[n] @ S[ypix[n], xpix[n], :]
        err = (Ucell**2).mean()
        t1 = time.time() - t0
        logger.info(f"iter {it},\tROIs: {ncells},\terr: {err:0.4f}, \ttime: {t1:0.2f} sec")

        if refine == 0:
            break
        if refine == 2:
            # good place to get connected regions
            stat = [{
                "ypix": ypix[n],
                "lam": lam[n],
                "xpix": xpix[n]
            } for n in range(ncells)]
            stat = connected_region(stat, connected, Ly, Lx)
            ypix = [stat[n]["ypix"] for n in range(len(stat))]
            xpix = [stat[n]["xpix"] for n in range(len(stat))]
            lam = [stat[n]["lam"] for n in range(len(stat))]
            ncells = len(ypix)
        if refine > 0:
            Ucell = Ucell + (S.reshape((-1, nbasis)) @ neu).reshape(U.shape)
        if refine < 0 and (newcells < Nfirst / 10 or it == max_iterations - 1):
            refine = 3
            U = getSVDproj(mov, u, diameter, smooth_masks)
            Ucell = U
        if refine >= 0:
            StU = S.reshape((Ly * Lx, -1)).transpose() @ Ucell.reshape(
                (Ly * Lx, -1))
            #StU = np.reshape(S, (Lyc*Lxc,-1)).transpose() @ np.reshape(Ucell, (Lyc*Lxc, -1))
            neu = np.linalg.solve(StS, StU).astype("float32")
        refine -= 1
    Ucell = U - (S.reshape((-1, nbasis)) @ neu).reshape(U.shape)

    sdmov = np.reshape(sdmov, (Ly, Lx))
    stat = [{
        "ypix": ypix[n],
        "lam": lam[n] * sdmov[ypix[n], xpix[n]],
        "xpix": xpix[n]
    } for n in range(ncells)]

    stat = connected_region(stat, connected, Ly, Lx)
    # Remove empty cells
    stat = [s for s in stat if len(s["ypix"]) != 0]
    stat = np.array(stat)

    new_settings = {
        "diameter": diameter,
        "Vmax": 0,
        "ihop": 0,
        "Vsplit": 0,
        "Vcorr": Vcorr,
        "Vmap": 0,
        "spatscale_pix": 0
    }
    return new_settings, stat

Copyright © 2023 Howard Hughes Medical Institute, Authored by Carsen Stringer and Marius Pachitariu.

add_square #

add_square(yi, xi, lx, Ly, Lx)

Create a square mask of pixels around a peak, normalized to unit norm.

Parameters:

Name Type Description Default
yi int

Y-coordinate of the center pixel.

 required
xi int

X-coordinate of the center pixel.

 required
lx int

Side length of the square in pixels.

 required
Ly int

Height of the full image.

 required
Lx int

Width of the full image.

 required

Returns:

Name Type Description
y0 ndarray

Y-coordinates of the square pixels, flattened.
x0 ndarray

X-coordinates of the square pixels, flattened.
mask ndarray

Pixel weights normalized to unit L2 norm, flattened.
Source code in suite2p/detection/sparsedetect.py 
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
def add_square(yi, xi, lx, Ly, Lx):
    """
    Create a square mask of pixels around a peak, normalized to unit norm.

    Parameters
    ----------
    yi : int
        Y-coordinate of the center pixel.
    xi : int
        X-coordinate of the center pixel.
    lx : int
        Side length of the square in pixels.
    Ly : int
        Height of the full image.
    Lx : int
        Width of the full image.

    Returns
    -------
    y0 : numpy.ndarray
        Y-coordinates of the square pixels, flattened.
    x0 : numpy.ndarray
        X-coordinates of the square pixels, flattened.
    mask : numpy.ndarray
        Pixel weights normalized to unit L2 norm, flattened.
    """
    lhf = int((lx - 1) / 2)
    ipix = np.tile(np.arange(-lhf, -lhf + lx, dtype=np.int32), reps=(lx, 1))
    x0 = xi + ipix
    y0 = yi + ipix.T
    mask = np.ones_like(ipix, dtype=np.float32)
    ix = np.all((y0 >= 0, y0 < Ly, x0 >= 0, x0 < Lx), axis=0)
    x0 = x0[ix]
    y0 = y0[ix]
    mask = mask[ix]
    mask = mask / norm(mask)
    return y0.flatten(), x0.flatten(), mask.flatten()

downsample #

downsample(mov, taper_edge=True)

Downsample a movie by 2x in both spatial dimensions.

Averages adjacent pixels along Y then X. If a dimension has odd size, the last pixel is halved when taper_edge is True.

Parameters:

Name Type Description Default
mov ndarray

Movie of shape (n_frames, Ly, Lx).

 required
taper_edge bool, optional (default True)

If True, halve the edge pixel when the dimension is odd.

 True

Returns:

Name Type Description
mov2 ndarray

Downsampled movie of shape (n_frames, ceil(Ly/2), ceil(Lx/2)).
Source code in suite2p/detection/sparsedetect.py 
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
def downsample(mov, taper_edge=True):
    """
    Downsample a movie by 2x in both spatial dimensions.

    Averages adjacent pixels along Y then X. If a dimension has odd size,
    the last pixel is halved when ``taper_edge`` is True.

    Parameters
    ----------
    mov : numpy.ndarray
        Movie of shape (n_frames, Ly, Lx).
    taper_edge : bool, optional (default True)
        If True, halve the edge pixel when the dimension is odd.

    Returns
    -------
    mov2 : numpy.ndarray
        Downsampled movie of shape (n_frames, ceil(Ly/2), ceil(Lx/2)).
    """
    n_frames, Ly, Lx = mov.shape

    # bin along Y
    movd = np.zeros((n_frames, int(np.ceil(Ly / 2)), Lx), "float32")
    movd[:, :Ly // 2, :] = np.mean([mov[:, 0:-1:2, :], mov[:, 1::2, :]], axis=0)
    if Ly % 2 == 1:
        movd[:, -1, :] = mov[:, -1, :] / 2 if taper_edge else mov[:, -1, :]

    # bin along X
    mov2 = np.zeros((n_frames, int(np.ceil(Ly / 2)), int(np.ceil(Lx / 2))), "float32")
    mov2[:, :, :Lx // 2] = np.mean([movd[:, :, 0:-1:2], movd[:, :, 1::2]], axis=0)
    if Lx % 2 == 1:
        mov2[:, :, -1] = movd[:, :, -1] / 2 if taper_edge else movd[:, :, -1]

    return mov2

estimate_spatial_scale #

estimate_spatial_scale(I)

Estimate the dominant spatial scale from multi-scale correlation maps.

Finds the scale index that appears most frequently among the top 50 brightest local maxima in the max-projected correlation image.

Parameters:

Name Type Description Default
I ndarray

Multi-scale correlation maps of shape (n_scales, Ly, Lx).

 required

Returns:

Name Type Description
im int

Index of the estimated best spatial scale.
Source code in suite2p/detection/sparsedetect.py 
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
def estimate_spatial_scale(I):
    """
    Estimate the dominant spatial scale from multi-scale correlation maps.

    Finds the scale index that appears most frequently among the top 50
    brightest local maxima in the max-projected correlation image.

    Parameters
    ----------
    I : numpy.ndarray
        Multi-scale correlation maps of shape (n_scales, Ly, Lx).

    Returns
    -------
    im : int
        Index of the estimated best spatial scale.
    """
    I0 = I.max(axis=0)
    imap = np.argmax(I, axis=0).flatten()
    ipk = np.abs(I0 - maximum_filter(I0, size=(11, 11))).flatten() < 1e-4
    isort = np.argsort(I0.flatten()[ipk])[::-1]
    im, _ = mode(imap[ipk][isort[:50]], keepdims=True)
    return im.item()

extendROI #

extendROI(ypix, xpix, Ly, Lx, niter=1)

Extend ROI pixel coordinates by growing the mask in 4-connected directions.

Parameters:

Name Type Description Default
ypix ndarray

Y-coordinates of the mask pixels.

 required
xpix ndarray

X-coordinates of the mask pixels.

 required
Ly int

Height of the image.

 required
Lx int

Width of the image.

 required
niter int, optional (default 1)

Number of dilation iterations.

 1

Returns:

Name Type Description
ypix ndarray

Extended Y-coordinates.
xpix ndarray

Extended X-coordinates.
Source code in suite2p/detection/sparsedetect.py 
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
def extendROI(ypix, xpix, Ly, Lx, niter=1):
    """
    Extend ROI pixel coordinates by growing the mask in 4-connected directions.

    Parameters
    ----------
    ypix : numpy.ndarray
        Y-coordinates of the mask pixels.
    xpix : numpy.ndarray
        X-coordinates of the mask pixels.
    Ly : int
        Height of the image.
    Lx : int
        Width of the image.
    niter : int, optional (default 1)
        Number of dilation iterations.

    Returns
    -------
    ypix : numpy.ndarray
        Extended Y-coordinates.
    xpix : numpy.ndarray
        Extended X-coordinates.
    """
    for k in range(niter):
        yx = ((ypix, ypix, ypix, ypix - 1, ypix + 1), (xpix, xpix + 1, xpix - 1, xpix,
                                                       xpix))
        yx = np.array(yx)
        yx = yx.reshape((2, -1))
        yu = np.unique(yx, axis=1)
        ix = np.all((yu[0] >= 0, yu[0] < Ly, yu[1] >= 0, yu[1] < Lx), axis=0)
        ypix, xpix = yu[:, ix]
    return ypix, xpix

extend_mask #

extend_mask(ypix, xpix, lam, Ly, Lx)

Extend a mask into the 8 surrounding pixels of each pixel.

Each pixel spreads its weight equally to itself and its 8 neighbors. Overlapping contributions are summed.

Parameters:

Name Type Description Default
ypix ndarray

Y-coordinates of the mask pixels.

 required
xpix ndarray

X-coordinates of the mask pixels.

 required
lam ndarray

Pixel weights.

 required
Ly int

Height of the image.

 required
Lx int

Width of the image.

 required

Returns:

Name Type Description
ypix1 ndarray

Y-coordinates of the extended mask.
xpix1 ndarray

X-coordinates of the extended mask.
lam1 ndarray

Accumulated pixel weights for the extended mask.
Source code in suite2p/detection/sparsedetect.py 
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
def extend_mask(ypix, xpix, lam, Ly, Lx):
    """
    Extend a mask into the 8 surrounding pixels of each pixel.

    Each pixel spreads its weight equally to itself and its 8 neighbors.
    Overlapping contributions are summed.

    Parameters
    ----------
    ypix : numpy.ndarray
        Y-coordinates of the mask pixels.
    xpix : numpy.ndarray
        X-coordinates of the mask pixels.
    lam : numpy.ndarray
        Pixel weights.
    Ly : int
        Height of the image.
    Lx : int
        Width of the image.

    Returns
    -------
    ypix1 : numpy.ndarray
        Y-coordinates of the extended mask.
    xpix1 : numpy.ndarray
        X-coordinates of the extended mask.
    lam1 : numpy.ndarray
        Accumulated pixel weights for the extended mask.
    """
    nel = len(xpix)
    yx = ((ypix, ypix, ypix, ypix - 1, ypix - 1, ypix - 1, ypix + 1, ypix + 1,
           ypix + 1), (xpix, xpix + 1, xpix - 1, xpix, xpix + 1, xpix - 1, xpix,
                       xpix + 1, xpix - 1))
    yx = np.array(yx)
    yx = yx.reshape((2, -1))
    yu, ind = np.unique(yx, axis=1, return_inverse=True)
    LAM = np.zeros(yu.shape[1])
    for j in range(len(ind)):
        LAM[ind[j]] += lam[j % nel] / 3
    ix = np.all((yu[0] >= 0, yu[0] < Ly, yu[1] >= 0, yu[1] < Lx), axis=0)
    ypix1, xpix1 = yu[:, ix]
    lam1 = LAM[ix]
    return ypix1, xpix1, lam1

find_best_scale #

find_best_scale(I, spatial_scale)

Determine the best spatial scale, either forced or estimated from data.

If spatial_scale is positive, clamps it to [1, 4] and returns it as forced. Otherwise estimates the scale from the multi-scale correlation maps.

Parameters:

Name Type Description Default
I ndarray

Multi-scale correlation maps of shape (n_scales, Ly, Lx).

 required
spatial_scale int

User-specified spatial scale. If positive, used directly (forced). If zero or negative, the scale is estimated from the data.

 required

Returns:

Name Type Description
scale int

Best spatial scale index.
estimate_mode EstimateMode

Whether the scale was EstimateMode.Forced or EstimateMode.Estimated.
Source code in suite2p/detection/sparsedetect.py 
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
def find_best_scale(I, spatial_scale):
    """
    Determine the best spatial scale, either forced or estimated from data.

    If ``spatial_scale`` is positive, clamps it to [1, 4] and returns it as
    forced. Otherwise estimates the scale from the multi-scale correlation
    maps.

    Parameters
    ----------
    I : numpy.ndarray
        Multi-scale correlation maps of shape (n_scales, Ly, Lx).
    spatial_scale : int
        User-specified spatial scale. If positive, used directly (forced).
        If zero or negative, the scale is estimated from the data.

    Returns
    -------
    scale : int
        Best spatial scale index.
    estimate_mode : EstimateMode
        Whether the scale was ``EstimateMode.Forced`` or
        ``EstimateMode.Estimated``.
    """
    if spatial_scale > 0:
        return max(1, min(4, spatial_scale)), EstimateMode.Forced
    else:
        scale = estimate_spatial_scale(I=I)
        if scale > 0:
            return scale, EstimateMode.Estimated
        else:
            warn(
                "Spatial scale estimation failed.  Setting spatial scale to 1 in order to continue."
            )
            return 1, EstimateMode.Forced

iter_extend #

iter_extend(ypix, xpix, mov, Lyc, Lxc, active_frames)

Iteratively extend a mask based on pixel activity on active frames.

Repeatedly grows the ROI by one pixel on each side, keeping only pixels whose mean activity on active frames exceeds 1/5 of the maximum pixel. Stops when the mask stops growing or reaches 10000 pixels.

Parameters:

Name Type Description Default
ypix ndarray

Y-coordinates of the initial mask pixels.

 required
xpix ndarray

X-coordinates of the initial mask pixels.

 required
mov ndarray

Binned residual movie of shape (nbinned, Lyc * Lxc).

 required
Lyc int

Height of the movie frame.

 required
Lxc int

Width of the movie frame.

 required
active_frames ndarray

Indices of frames used to compute activity.

 required

Returns:

Name Type Description
ypix ndarray

Y-coordinates of the extended mask.
xpix ndarray

X-coordinates of the extended mask.
lam ndarray

Pixel weights normalized to unit L2 norm.
Source code in suite2p/detection/sparsedetect.py 
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
def iter_extend(ypix, xpix, mov, Lyc, Lxc, active_frames):
    """
    Iteratively extend a mask based on pixel activity on active frames.

    Repeatedly grows the ROI by one pixel on each side, keeping only pixels
    whose mean activity on active frames exceeds 1/5 of the maximum pixel. 
    Stops when the mask stops growing or reaches 10000 pixels.

    Parameters
    ----------
    ypix : numpy.ndarray
        Y-coordinates of the initial mask pixels.
    xpix : numpy.ndarray
        X-coordinates of the initial mask pixels.
    mov : numpy.ndarray
        Binned residual movie of shape (nbinned, Lyc * Lxc).
    Lyc : int
        Height of the movie frame.
    Lxc : int
        Width of the movie frame.
    active_frames : numpy.ndarray
        Indices of frames used to compute activity.

    Returns
    -------
    ypix : numpy.ndarray
        Y-coordinates of the extended mask.
    xpix : numpy.ndarray
        X-coordinates of the extended mask.
    lam : numpy.ndarray
        Pixel weights normalized to unit L2 norm.
    """
    npix = 0
    iter = 0
    while npix < 10000:
        npix = ypix.size
        # extend ROI by 1 pixel on each side
        ypix, xpix = extendROI(ypix, xpix, Lyc, Lxc, 1)
        # activity in proposed ROI on ACTIVE frames
        usub = mov[np.ix_(active_frames, ypix * Lxc + xpix)]
        lam = usub.mean(axis=0)
        ix = lam > max(0, lam.max() / 5.0)
        if ix.sum() == 0:
            break
        ypix, xpix, lam = ypix[ix], xpix[ix], lam[ix]
        if iter == 0:
            sgn = 1.
        if np.sign(sgn * (ix.sum() - npix)) <= 0:
            break
        else:
            npix = ypix.size
        iter += 1
    lam = lam / np.sum(lam**2)**.5
    return ypix, xpix, lam

multiscale_mask #

multiscale_mask(ypix0, xpix0, lam0, Lyp, Lxp)

Downsample a mask to all spatial scales used in sparse detection.

Given pixel coordinates and weights at the original resolution, creates downsampled versions by successively halving coordinates and accumulating weights, then extends each mask into surrounding pixels.

Parameters:

Name Type Description Default
ypix0 ndarray

Y-coordinates of the mask pixels at full resolution.

 required
xpix0 ndarray

X-coordinates of the mask pixels at full resolution.

 required
lam0 ndarray

Pixel weights at full resolution.

 required
Lyp ndarray

Heights of the downsampled images at each scale, shape (n_scales,).

 required
Lxp ndarray

Widths of the downsampled images at each scale, shape (n_scales,).

 required

Returns:

Name Type Description
ys list of numpy.ndarray

Y-coordinates of the mask at each spatial scale.
xs list of numpy.ndarray

X-coordinates of the mask at each spatial scale.
lms list of numpy.ndarray

Pixel weights at each spatial scale.
Source code in suite2p/detection/sparsedetect.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
def multiscale_mask(ypix0, xpix0, lam0, Lyp, Lxp):
    """
    Downsample a mask to all spatial scales used in sparse detection.

    Given pixel coordinates and weights at the original resolution, creates
    downsampled versions by successively halving coordinates and accumulating
    weights, then extends each mask into surrounding pixels.

    Parameters
    ----------
    ypix0 : numpy.ndarray
        Y-coordinates of the mask pixels at full resolution.
    xpix0 : numpy.ndarray
        X-coordinates of the mask pixels at full resolution.
    lam0 : numpy.ndarray
        Pixel weights at full resolution.
    Lyp : numpy.ndarray
        Heights of the downsampled images at each scale, shape (n_scales,).
    Lxp : numpy.ndarray
        Widths of the downsampled images at each scale, shape (n_scales,).

    Returns
    -------
    ys : list of numpy.ndarray
        Y-coordinates of the mask at each spatial scale.
    xs : list of numpy.ndarray
        X-coordinates of the mask at each spatial scale.
    lms : list of numpy.ndarray
        Pixel weights at each spatial scale.
    """
    xs = [xpix0]
    ys = [ypix0]
    lms = [lam0]
    for j in range(1, len(Lyp)):
        ipix, ind = np.unique(
            np.int32(xs[j - 1] / 2) + np.int32(ys[j - 1] / 2) * Lxp[j],
            return_inverse=True)
        LAM = np.zeros(len(ipix))
        for i in range(len(xs[j - 1])):
            LAM[ind[i]] += lms[j - 1][i] / 2
        lms.append(LAM)
        ys.append(np.int32(ipix / Lxp[j]))
        xs.append(np.int32(ipix % Lxp[j]))
    for j in range(len(Lyp)):
        ys[j], xs[j], lms[j] = extend_mask(ys[j], xs[j], lms[j], Lyp[j], Lxp[j])
    return ys, xs, lms

neuropil_subtraction #

neuropil_subtraction(mov, filter_size)

Subtract a spatially low-pass filtered version of the movie to remove neuropil.

Parameters:

Name Type Description Default
mov ndarray

Movie of shape (nbinned, Ly, Lx).

 required
filter_size int

Width of the uniform spatial filter in pixels.

 required

Returns:

Name Type Description
movt ndarray

High-pass filtered movie of shape (nbinned, Ly, Lx).
Source code in suite2p/detection/sparsedetect.py 
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
def neuropil_subtraction(mov, filter_size):
    """
    Subtract a spatially low-pass filtered version of the movie to remove neuropil.

    Parameters
    ----------
    mov : numpy.ndarray
        Movie of shape (nbinned, Ly, Lx).
    filter_size : int
        Width of the uniform spatial filter in pixels.

    Returns
    -------
    movt : numpy.ndarray
        High-pass filtered movie of shape (nbinned, Ly, Lx).
    """
    nbinned, Ly, Lx = mov.shape
    c1 = uniform_filter(np.ones((Ly, Lx)), size=filter_size, mode="constant")
    movt = np.zeros_like(mov)
    for frame, framet in zip(mov, movt):
        framet[:] = frame - (uniform_filter(frame, size=filter_size, mode="constant") /
                             c1)
    return movt

sparsery #

sparsery(mov, sdmov, highpass_neuropil, spatial_scale, threshold_scaling, max_ROIs, active_percentile=0)

Detect ROIs in a movie using doubly-sparse matrix decomposition.

Subtracts neuropil, builds multi-scale correlation maps, then greedily chooses top activity peak, extends mask, optionally splits the ROI, and subtracts the detected signal from the residual movie, then extracts the highest peak from the residual and continues.

Parameters:

Name Type Description Default
mov ndarray

Binned movie of shape (nbinned, Ly, Lx).

 required
sdmov ndarray

Per-pixel standard deviation of shape (Ly, Lx), used for normalization.

 required
highpass_neuropil int

Filter size for spatial high-pass neuropil subtraction.

 required
spatial_scale int

Spatial scale setting. If positive, forced; if zero or negative, estimated from the data.

 required
threshold_scaling float

Multiplier for the activity threshold used to accept peaks.

 required
max_ROIs int

Maximum number of ROIs to detect.

 required
active_percentile float, optional (default 0)

If positive, use this percentile of the temporal projection as an alternative activity threshold.

 0

Returns:

Name Type Description
new_settings dict

Dictionary with detection metadata including "Vmax", "ihop", "Vsplit", "Vcorr", "Vmap", and "spatscale_pix".
stats list of dict

List of ROI statistics dictionaries, each containing "ypix", "xpix", "lam", "med", and "footprint".
Source code in suite2p/detection/sparsedetect.py 
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
def sparsery(mov, sdmov, highpass_neuropil,
             spatial_scale, threshold_scaling, max_ROIs,
             active_percentile=0):
    """
    Detect ROIs in a movie using doubly-sparse matrix decomposition.

    Subtracts neuropil, builds multi-scale correlation maps, then
    greedily chooses top activity peak, extends mask, optionally splits
    the ROI, and subtracts the detected signal from the residual movie, then 
    extracts the highest peak from the residual and continues.

    Parameters
    ----------
    mov : numpy.ndarray
        Binned movie of shape (nbinned, Ly, Lx).
    sdmov : numpy.ndarray
        Per-pixel standard deviation of shape (Ly, Lx), used for
        normalization.
    highpass_neuropil : int
        Filter size for spatial high-pass neuropil subtraction.
    spatial_scale : int
        Spatial scale setting. If positive, forced; if zero or negative,
        estimated from the data.
    threshold_scaling : float
        Multiplier for the activity threshold used to accept peaks.
    max_ROIs : int
        Maximum number of ROIs to detect.
    active_percentile : float, optional (default 0)
        If positive, use this percentile of the temporal projection as an
        alternative activity threshold.

    Returns
    -------
    new_settings : dict
        Dictionary with detection metadata including "Vmax", "ihop",
        "Vsplit", "Vcorr", "Vmap", and "spatscale_pix".
    stats : list of dict
        List of ROI statistics dictionaries, each containing "ypix",
        "xpix", "lam", "med", and "footprint".
    """

    mov = neuropil_subtraction(
        mov=mov / sdmov,
        filter_size=highpass_neuropil)  # subtract low-pass filtered movie

    _, Lyc, Lxc = mov.shape
    LL = np.meshgrid(np.arange(Lxc), np.arange(Lyc))
    gxy = [np.array(LL).astype("float32")]
    dmov = mov
    movu = []

    # downsample movie at various spatial scales
    Lyp, Lxp = np.zeros(5, "int32"), np.zeros(5, "int32")  # downsampled sizes
    for j in range(5):
        movu0 = square_convolution_2d(dmov, 3)
        dmov = 2 * downsample(dmov)
        gxy0 = downsample(gxy[j], False)
        gxy.append(gxy0)
        _, Lyp[j], Lxp[j] = movu0.shape
        movu.append(movu0)

    # spline over scales
    I = np.zeros((len(gxy), gxy[0].shape[1], gxy[0].shape[2]))
    for movu0, gxy0, I0 in zip(movu, gxy, I):
        gmodel = RectBivariateSpline(gxy0[1, :, 0], gxy0[0, 0, :], movu0.max(axis=0),
                                     kx=min(3, gxy0.shape[1] - 1),
                                     ky=min(3, gxy0.shape[2] - 1))
        I0[:] = gmodel(gxy[0][1, :, 0], gxy[0][0, 0, :])
    v_corr = I.max(axis=0)

    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)
    vmultiplier = max(1, mov.shape[0] / 1200)
    logger.info("NOTE: %s spatial scale ~%d pixels, time epochs %2.2f, threshold %2.2f " %
          (estimate_mode.value, spatscale_pix, vmultiplier, vmultiplier * Th2))

    # get standard deviation for pixels for all values > Th2
    v_map = [threshold_reduce(movu0, Th2) for movu0 in movu]
    movu = [movu0.reshape(movu0.shape[0], -1) for movu0 in movu]

    mov = np.reshape(mov, (-1, Lyc * Lxc))
    lxs = 3 * 2**np.arange(5)
    nscales = len(lxs)

    v_max = np.zeros(max_ROIs)
    ihop = np.zeros(max_ROIs)
    v_split = np.zeros(max_ROIs)
    V1 = deepcopy(v_map)
    stats = []
    patches = []
    seeds = []
    extract_patches = False

    logger.info(f"max_ROIs set to {max_ROIs} - will run for {max_ROIs} ROIs or until no more ROIs above threshold are found.")
    t0 = time.time()
    for tj in range(max_ROIs):
        # find peaks in stddev"s
        v0max = np.array([V1[j].max() for j in range(5)])
        imap = np.argmax(v0max)
        imax = np.argmax(V1[imap])
        yi, xi = np.unravel_index(imax, (Lyp[imap], Lxp[imap]))
        # position of peak
        yi, xi = gxy[imap][1, yi, xi], gxy[imap][0, yi, xi]
        med = [int(yi), int(xi)]

        # check if peak is larger than threshold * max(1,nbinned/1200)
        v_max[tj] = v0max.max()
        if v_max[tj] < vmultiplier * Th2:
            break
        ls = lxs[imap]

        ihop[tj] = imap

        # make square of initial pixels based on spatial scale of peak
        yi, xi = int(yi), int(xi)
        ypix0, xpix0, lam0 = add_square(yi, xi, ls, Lyc, Lxc)

        # project movie into square to get time series
        tproj = (mov[:, ypix0 * Lxc + xpix0] * lam0[0]).sum(axis=-1)
        if active_percentile > 0:
            threshold = min(Th2, np.percentile(tproj, active_percentile))
        else:
            threshold = Th2
        active_frames = np.nonzero(tproj > threshold)[0]  # frames with activity > Th2

        # get square around seed
        if extract_patches:
            mask = mov[active_frames].mean(axis=0).reshape(Lyc, Lxc)
            patches.append(utils.square_mask(mask, mask_window, yi, xi))
            seeds.append([yi, xi])

        # extend mask based on activity similarity
        for j in range(3):
            ypix0, xpix0, lam0 = iter_extend(ypix0, xpix0, mov, Lyc, Lxc, active_frames)
            tproj = mov[:, ypix0 * Lxc + xpix0] @ lam0
            active_frames = np.nonzero(tproj > threshold)[0]
            if len(active_frames) < 1:
                #if tj < max_ROIs/2: # TODO: nmasks is undefined
                #    continue
                #else:
                break

        if len(active_frames) < 1:
            #if tj < max_ROIs/2:
            #    continue
            #else:
            break

        # check if ROI should be split
        v_split[tj], ipack = two_comps(mov[:, ypix0 * Lxc + xpix0], lam0, threshold)
        if v_split[tj] > 1.25:
            lam0, xp, active_frames = ipack
            tproj[active_frames] = xp
            ix = lam0 > lam0.max() / 5
            xpix0 = xpix0[ix]
            ypix0 = ypix0[ix]
            lam0 = lam0[ix]
            ymed = np.median(ypix0)
            xmed = np.median(xpix0)
            imin = np.argmin((xpix0 - xmed)**2 + (ypix0 - ymed)**2)
            med = [ypix0[imin], xpix0[imin]]

        # update residual on raw movie
        mov[np.ix_(active_frames,
                   ypix0 * Lxc + xpix0)] -= tproj[active_frames][:, np.newaxis] * lam0
        # update filtered movie
        ys, xs, lms = multiscale_mask(ypix0, xpix0, lam0, Lyp, Lxp)
        for j in range(nscales):
            movu[j][np.ix_(active_frames, xs[j] + Lxp[j] * ys[j])] -= np.outer(
                tproj[active_frames], lms[j])
            Mx = movu[j][:, xs[j] + Lxp[j] * ys[j]]
            V1[j][ys[j], xs[j]] = (Mx**2 * np.float32(Mx > threshold)).sum(axis=0)**.5

        stats.append({
            "ypix": ypix0.astype(int),
            "xpix": xpix0.astype(int),
            "lam": lam0 * sdmov[ypix0, xpix0],
            "med": med,
            "footprint": ihop[tj]
        })

        if tj % 500 == 0:
            t1 = time.time() - t0
            logger.info(f"ROIs: {tj},\t last score: {v_max[tj]:0.4f}, \t time: {t1:0.2f}sec")


    new_settings = {
        "Vmax": v_max,
        "ihop": ihop,
        "Vsplit": v_split,
        "Vcorr": v_corr,
        "Vmap": np.asanyarray(
            v_map, dtype="object"
        ),  # needed so that scipy.io.savemat doesn"t fail in runpipeline with latest numpy (v1.24.3). dtype="object" is needed to have numpy array with elements having diff sizes
        "spatscale_pix": spatscale_pix,
    }

    return new_settings, stats

square_convolution_2d #

square_convolution_2d(mov, filter_size)

Convolve each frame with a square uniform kernel.

Parameters:

Name Type Description Default
mov ndarray

Movie of shape (nbinned, Ly, Lx).

 required
filter_size int

Width of the square uniform kernel in pixels.

 required

Returns:

Name Type Description
movt ndarray

Convolved movie of shape (nbinned, Ly, Lx), dtype float32.
Source code in suite2p/detection/sparsedetect.py 
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
def square_convolution_2d(mov, filter_size):
    """
    Convolve each frame with a square uniform kernel.

    Parameters
    ----------
    mov : numpy.ndarray
        Movie of shape (nbinned, Ly, Lx).
    filter_size : int
        Width of the square uniform kernel in pixels.

    Returns
    -------
    movt : numpy.ndarray
        Convolved movie of shape (nbinned, Ly, Lx), dtype float32.
    """
    movt = np.zeros_like(mov, dtype=np.float32)
    for frame, framet in zip(mov, movt):
        framet[:] = filter_size * uniform_filter(frame, size=filter_size,
                                                 mode="constant")
    return movt

threshold_reduce #

threshold_reduce(mov, intensity_threshold)

Compute thresholded standard deviation map across frames.

For each pixel, sums the squared values of frames exceeding the threshold, then takes the square root. Iterates frame-by-frame to reduce memory usage.

Parameters:

Name Type Description Default
mov ndarray

Movie of shape (nbinned, Ly, Lx).

 required
intensity_threshold float

Only frames where the pixel value exceeds this are included.

 required

Returns:

Name Type Description
Vt ndarray

Thresholded standard deviation map of shape (Ly, Lx), dtype float32.
Source code in suite2p/detection/sparsedetect.py 
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
def threshold_reduce(mov, intensity_threshold):
    """
    Compute thresholded standard deviation map across frames.

    For each pixel, sums the squared values of frames exceeding
    the threshold, then takes the square root. Iterates frame-by-frame
    to reduce memory usage.

    Parameters
    ----------
    mov : numpy.ndarray
        Movie of shape (nbinned, Ly, Lx).
    intensity_threshold : float
        Only frames where the pixel value exceeds this are included.

    Returns
    -------
    Vt : numpy.ndarray
        Thresholded standard deviation map of shape (Ly, Lx), dtype float32.
    """
    nbinned, Lyp, Lxp = mov.shape
    Vt = np.zeros((Lyp, Lxp), "float32")
    for t in range(nbinned):
        Vt += mov[t]**2 * (mov[t] > intensity_threshold)
    Vt = Vt**.5
    return Vt

two_comps #

two_comps(mpix0, lam, Th2)

Check if splitting an ROI into two components increases variance explained.

Projects the ROI movie onto the current mask, then tests whether a two-component split captures more variance. Returns the variance ratio and the better component.

Parameters:

Name Type Description Default
mpix0 ndarray

Binned movie for pixels in the ROI, shape (nbinned, npix).

 required
lam ndarray

Pixel weights for the current ROI.

 required
Th2 float

Intensity threshold for determining active frames.

 required

Returns:

Name Type Description
vrat float

Ratio of variance explained by two components to one component. Values above 1.25 suggest the ROI should be split.
ipick tuple

Tuple of (mu, xproj, goodframe) for the better component, where mu is the pixel weights, xproj is the temporal projection on active frames, and goodframe is the boolean active-frame mask.
Source code in suite2p/detection/sparsedetect.py 
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
def two_comps(mpix0, lam, Th2):
    """
    Check if splitting an ROI into two components increases variance explained.

    Projects the ROI movie onto the current mask, then tests whether a
    two-component split captures more variance. Returns the variance ratio
    and the better component.

    Parameters
    ----------
    mpix0 : numpy.ndarray
        Binned movie for pixels in the ROI, shape (nbinned, npix).
    lam : numpy.ndarray
        Pixel weights for the current ROI.
    Th2 : float
        Intensity threshold for determining active frames.

    Returns
    -------
    vrat : float
        Ratio of variance explained by two components to one component.
        Values above 1.25 suggest the ROI should be split.
    ipick : tuple
        Tuple of (mu, xproj, goodframe) for the better component, where
        mu is the pixel weights, xproj is the temporal projection on active
        frames, and goodframe is the boolean active-frame mask.
    """
    mpix = mpix0.copy()
    xproj = mpix @ lam
    gf0 = xproj > Th2

    mpix[gf0, :] -= np.outer(xproj[gf0], lam)
    vexp0 = np.sum(mpix0**2) - np.sum(mpix**2)

    k = np.argmax(np.sum(mpix * np.float32(mpix > 0), axis=1))
    mu = [lam * np.float32(mpix[k] < 0), lam * np.float32(mpix[k] > 0)]

    mpix = mpix0.copy()
    goodframe = []
    xproj = []
    for mu0 in mu:
        mu0[:] /= norm(mu0) + 1e-6
        xp = mpix @ mu0
        mpix[gf0, :] -= np.outer(xp[gf0], mu0)
        goodframe.append(gf0)
        xproj.append(xp[gf0])

    flag = [False, False]
    V = np.zeros(2)
    for t in range(3):
        for k in range(2):
            if flag[k]:
                continue
            mpix[goodframe[k], :] += np.outer(xproj[k], mu[k])
            xp = mpix @ mu[k]
            goodframe[k] = xp > Th2
            V[k] = np.sum(xp**2)
            if np.sum(goodframe[k]) == 0:
                flag[k] = True
                V[k] = -1
                continue
            xproj[k] = xp[goodframe[k]]
            mu[k] = np.mean(mpix[goodframe[k], :] * xproj[k][:, np.newaxis], axis=0)
            mu[k][mu[k] < 0] = 0
            mu[k] /= (1e-6 + np.sum(mu[k]**2)**.5)
            mpix[goodframe[k], :] -= np.outer(xproj[k], mu[k])
    k = np.argmax(V)
    vexp = np.sum(mpix0**2) - np.sum(mpix**2)
    vrat = vexp / vexp0
    return vrat, (mu[k], xproj[k], goodframe[k])

Copyright © 2023 Howard Hughes Medical Institute, Authored by Carsen Stringer and Marius Pachitariu.

assign_overlaps #

assign_overlaps(stats, Ly, Lx)

Assign overlap labels to each ROI based on shared pixels.

For each ROI, sets an "overlap" boolean mask indicating which of its pixels are shared with at least one other ROI.

Parameters:

Name Type Description Default
stats ndarray

Array of ROI statistics dictionaries, each containing "ypix" and "xpix".

 required
Ly int

Height of the image in pixels.

 required
Lx int

Width of the image in pixels.

 required

Returns:

Name Type Description
stats ndarray

Updated array with "overlap" key added to each ROI dictionary.
Source code in suite2p/detection/stats.py 
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
def assign_overlaps(stats, Ly, Lx):
    """
    Assign overlap labels to each ROI based on shared pixels.

    For each ROI, sets an "overlap" boolean mask indicating which of its pixels
    are shared with at least one other ROI.

    Parameters
    ----------
    stats : numpy.ndarray
        Array of ROI statistics dictionaries, each containing "ypix" and "xpix".
    Ly : int
        Height of the image in pixels.
    Lx : int
        Width of the image in pixels.

    Returns
    -------
    stats : numpy.ndarray
        Updated array with "overlap" key added to each ROI dictionary.
    """
    overlap = np.zeros((Ly, Lx), "int")
    for stat in stats:
        overlap[stat["ypix"], stat["xpix"]] += 1
    for stat in stats:
        stat["overlap"] = (overlap[stat["ypix"], stat["xpix"]] > 1).astype("bool")
    return stats

fitMVGaus #

fitMVGaus(y, x, lam0, dy, dx, thres=2.5, npts=100)

Fit a 2D Gaussian to weighted pixel coordinates and return an ellipse.

Computes the mean and covariance of the pixel distribution weighted by lam0, then generates an ellipse at thres standard deviations.

Parameters:

Name Type Description Default
y ndarray

Y-coordinates of the pixels.

 required
x ndarray

X-coordinates of the pixels.

 required
lam0 ndarray

Weights for each pixel (e.g. fluorescence intensity).

 required
dy float

Normalization factor for the y-axis (e.g. cell diameter in y).

 required
dx float

Normalization factor for the x-axis (e.g. cell diameter in x).

 required
thres float, optional (default 2.5)

Number of standard deviations for the ellipse radius.

 2.5
npts int, optional (default 100)

Number of points used to draw the ellipse.

 100

Returns:

Name Type Description
mu ndarray

Mean of the Gaussian fit, shape (2,).
cov ndarray

Covariance matrix of the Gaussian fit, shape (2, 2).
radii ndarray

Radii of the major and minor axes (sorted descending), shape (2,).
ellipse ndarray

Points on the fitted ellipse, shape (npts, 2).
Source code in suite2p/detection/stats.py 
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
def fitMVGaus(y, x, lam0, dy, dx, thres=2.5, npts: int = 100):
    """
    Fit a 2D Gaussian to weighted pixel coordinates and return an ellipse.

    Computes the mean and covariance of the pixel distribution weighted by `lam0`,
    then generates an ellipse at `thres` standard deviations.

    Parameters
    ----------
    y : numpy.ndarray
        Y-coordinates of the pixels.
    x : numpy.ndarray
        X-coordinates of the pixels.
    lam0 : numpy.ndarray
        Weights for each pixel (e.g. fluorescence intensity).
    dy : float
        Normalization factor for the y-axis (e.g. cell diameter in y).
    dx : float
        Normalization factor for the x-axis (e.g. cell diameter in x).
    thres : float, optional (default 2.5)
        Number of standard deviations for the ellipse radius.
    npts : int, optional (default 100)
        Number of points used to draw the ellipse.

    Returns
    -------
    mu : numpy.ndarray
        Mean of the Gaussian fit, shape (2,).
    cov : numpy.ndarray
        Covariance matrix of the Gaussian fit, shape (2, 2).
    radii : numpy.ndarray
        Radii of the major and minor axes (sorted descending), shape (2,).
    ellipse : numpy.ndarray
        Points on the fitted ellipse, shape (npts, 2).
    """
    y = y / dy
    x = x / dx

    # normalize pixel weights
    lam = lam0.copy()
    ix = lam > 0  #lam.max()/5
    y, x, lam = y[ix], x[ix], lam[ix]
    lam /= lam.sum()

    # mean of gaussian
    yx = np.stack((y, x))
    mu = (lam * yx).sum(axis=1)
    yx = (yx - mu[:, np.newaxis]) * lam**.5
    cov = yx @ yx.T

    # radii of major and minor axes
    radii, evec = np.linalg.eig(cov)
    radii = thres * np.maximum(0, np.real(radii))**.5

    # compute pts of ellipse
    theta = np.linspace(0, 2 * np.pi, npts)
    p = np.stack((np.cos(theta), np.sin(theta)))
    ellipse = (p.T * radii) @ evec.T + mu
    radii = np.sort(radii)[::-1]
    return mu, cov, radii, ellipse

median_pix #

median_pix(ypix, xpix)

Find the pixel closest to the median of a set of pixel coordinates.

Parameters:

Name Type Description Default
ypix ndarray

Y-coordinates of the pixels.

 required
xpix ndarray

X-coordinates of the pixels.

 required

Returns:

Name Type Description
med list of float

Two-element list [ymed, xmed] of the pixel closest to the median.
Source code in suite2p/detection/stats.py 
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
def median_pix(ypix, xpix):
    """
    Find the pixel closest to the median of a set of pixel coordinates.

    Parameters
    ----------
    ypix : numpy.ndarray
        Y-coordinates of the pixels.
    xpix : numpy.ndarray
        X-coordinates of the pixels.

    Returns
    -------
    med : list of float
        Two-element list [ymed, xmed] of the pixel closest to the median.
    """
    ymed, xmed = np.median(ypix), np.median(xpix)
    imin = ((xpix - xmed)**2 + (ypix - ymed)**2).argmin()
    xmed = xpix[imin]
    ymed = ypix[imin]
    return [ymed, xmed]

roi_stats #

roi_stats(stats, Ly, Lx, diameter=[12.0, 12.0], max_overlap=0.75, do_soma_crop=True, npix_norm_min=-1, npix_norm_max=np.inf, median=False)

Compute statistics for detected ROIs, including compactness, aspect ratio, and overlap.

For each ROI, computes the median center, soma crop, compactness, and aspect ratio from a 2D Gaussian fit. Normalizes pixel counts across ROIs, removes ROIs outside the normalized pixel range, and optionally removes ROIs with excessive overlap.

Parameters:

Name Type Description Default
stats ndarray

Array of dictionaries, each containing "ypix", "xpix", and "lam" for one detected ROI.

 required
Ly int

Height of the image in pixels.

 required
Lx int

Width of the image in pixels.

 required
diameter list of float, optional (default [12., 12.])

Expected cell diameter [dy, dx] in pixels, used for normalization.

 [12.0, 12.0]
max_overlap float, optional (default 0.75)

Maximum allowed fraction of overlapping pixels. ROIs exceeding this are removed. Set to None or 1.0 to disable.

 0.75
do_soma_crop bool, optional (default True)

If True, crop dendritic pixels before computing compactness and aspect ratio.

 True
npix_norm_min (float, optional(default - 1))

Minimum normalized pixel count. ROIs below this are removed.

 -1
npix_norm_max float, optional (default np.inf)

Maximum normalized pixel count. ROIs above this are removed.

 inf
median bool, optional (default False)

If True, use median of all ROIs for normalization. If False, use median of the 100 ROIs extracted first.

 False

Returns:

Name Type Description
stats ndarray

Updated array of ROI statistics dictionaries with added keys "med", "npix", "soma_crop", "npix_soma", "mrs", "mrs0", "compact", "radius", "aspect_ratio", "footprint", "npix_norm", "npix_norm_no_crop", and "overlap".
Source code in suite2p/detection/stats.py 
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
def roi_stats(stats, Ly: int, Lx: int, diameter=[12., 12.], max_overlap=0.75,
              do_soma_crop=True, npix_norm_min=-1, npix_norm_max=np.inf,
              median=False):
    """
    Compute statistics for detected ROIs, including compactness, aspect ratio, and overlap.

    For each ROI, computes the median center, soma crop, compactness, and aspect
    ratio from a 2D Gaussian fit. Normalizes pixel counts across ROIs, removes
    ROIs outside the normalized pixel range, and optionally removes ROIs with
    excessive overlap.

    Parameters
    ----------
    stats : numpy.ndarray
        Array of dictionaries, each containing "ypix", "xpix", and "lam" for
        one detected ROI.
    Ly : int
        Height of the image in pixels.
    Lx : int
        Width of the image in pixels.
    diameter : list of float, optional (default [12., 12.])
        Expected cell diameter [dy, dx] in pixels, used for normalization.
    max_overlap : float, optional (default 0.75)
        Maximum allowed fraction of overlapping pixels. ROIs exceeding this
        are removed. Set to None or 1.0 to disable.
    do_soma_crop : bool, optional (default True)
        If True, crop dendritic pixels before computing compactness and
        aspect ratio.
    npix_norm_min : float, optional (default -1)
        Minimum normalized pixel count. ROIs below this are removed.
    npix_norm_max : float, optional (default np.inf)
        Maximum normalized pixel count. ROIs above this are removed.
    median : bool, optional (default False)
        If True, use median of all ROIs for normalization. If False, use
        median of the 100 ROIs extracted first.

    Returns
    -------
    stats : numpy.ndarray
        Updated array of ROI statistics dictionaries with added keys "med",
        "npix", "soma_crop", "npix_soma", "mrs", "mrs0", "compact", "radius",
        "aspect_ratio", "footprint", "npix_norm", "npix_norm_no_crop", and
        "overlap".
    """

    # approx size of masks for ROI aspect ratio estimation
    dy, dx = diameter[0], diameter[1]
    d0 = np.array([float(dy), float(dx)])

    dy, dx = np.meshgrid(np.arange(-d0[0]*3, d0[0]*3 + 1) / d0[0], 
                         np.arange(-d0[1]*3, d0[1]*3 + 1) / d0[1], indexing="ij")
    rs = (dy**2 + dx**2)**0.5
    dists_disk = np.sort(rs.flatten())

    for k, stat in enumerate(stats):
        ypix, xpix, lam = stat["ypix"].copy(), stat["xpix"].copy(), stat["lam"].copy()
        med = stat.get("med", median_pix(ypix, xpix)) # median of pixels
        stat["med"], stat["npix"] = med, ypix.size

        # crop dendritic pixels out for computing compactness and aspect ratio
        if do_soma_crop:
            crop = soma_crop(ypix, xpix, lam, med)
            stat["soma_crop"] = crop
            ypix, xpix, lam = ypix[crop], xpix[crop], lam[crop]
        else:
            stat["soma_crop"] = np.ones(ypix.size, "bool")
        stat["npix_soma"] = stat["soma_crop"].sum()  

        # compute compactness of ROI
        med = np.median(ypix), np.median(xpix)
        dists = (((ypix - med[0]) / d0[0])**2 + ((xpix - med[1]) / d0[1])**2)**0.5
        stat["mrs"], stat["mrs0"] = dists.mean(), dists_disk[:ypix.size].mean()
        stat["compact"] = max(1.0, stat["mrs"] / (1e-10 + stat["mrs0"]))

        # compute aspect ratio
        if "radius" not in stat:
            radii = fitMVGaus(ypix, xpix, lam, dy=d0[0], dx=d0[1], thres=2)[2]
            stat["radius"] = radii[0] * d0.mean()
            stat["aspect_ratio"] = 2 * radii[0] / (.01 + radii[0] + radii[1])
        stat["footprint"] = stat.get("footprint", 0)

    ### compute npix_norm (normalized npix) for each ROI
    npix_soma = np.array([stat["npix_soma"] for stat in stats], dtype="float32")
    npix = np.array([stat["npix"] for stat in stats], dtype="float32")
    # use median if cellpose, otherwise use best neurons to determine normalizer
    norm_npix = np.median(npix_soma) if median else np.median(npix_soma[:100])
    npix_soma /= norm_npix + 1e-10
    norm_npix = np.median(npix) if median else np.median(npix[:100])
    npix /= norm_npix + 1e-10

    keep_rois = (npix_norm_min <= npix_soma) * (npix_soma <= npix_norm_max)
    stats = stats[keep_rois]
    npix_soma, npix = npix_soma[keep_rois], npix[keep_rois]
    nremove = (~keep_rois).sum()
    if nremove > 0:
        logger.info(f"Removed {nremove} ROIs with npix_norm < {npix_norm_min:.2f} or npix_norm > {npix_norm_max:.2f}")

    for stat, npix_soma0, npix0 in zip(stats, npix_soma, npix):
        stat["npix_norm"] = npix_soma0
        stat["npix_norm_no_crop"] = npix0       

    if max_overlap is not None and max_overlap < 1.0:
        overlap = np.zeros((Ly, Lx), "int")
        for stat in stats:
            overlap[stat["ypix"], stat["xpix"]] += 1

        keep_rois = np.zeros(len(stats), "bool")
        # remove overlapping ROIs in reversed order, because highest variance ROIs are first
        for k, stat in enumerate(stats[::-1]): 
            keep_roi = (overlap[stat["ypix"], stat["xpix"]] > 1).mean() <= max_overlap
            keep_rois[k] = keep_roi
            if not keep_roi:
                overlap[stat["ypix"], stat["xpix"]] -= 1
        keep_rois = keep_rois[::-1]
        stats = stats[keep_rois]

        for stat in stats:
            stat["overlap"] = (overlap[stat["ypix"], stat["xpix"]] > 1).astype("bool")

        nremove = (~keep_rois).sum()
        logger.info(f"Removed {nremove} ROIs with overlap > {max_overlap}")

    return stats

soma_crop #

soma_crop(ypix, xpix, lam, med)

Crop dendritic pixels from an ROI by finding the soma boundary.

Computes cumulative weighted area as a function of distance from the median center, then finds the radius where the area growth drops below a threshold.

Parameters:

Name Type Description Default
ypix ndarray

Y-coordinates of the ROI pixels.

 required
xpix ndarray

X-coordinates of the ROI pixels.

 required
lam ndarray

Weights (e.g. fluorescence) for each pixel.

 required
med list of float

Two-element list [ymed, xmed] of the ROI center.

 required

Returns:

Name Type Description
crop ndarray

Boolean array of length len(ypix), True for pixels within the soma.
Source code in suite2p/detection/stats.py 
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
def soma_crop(ypix, xpix, lam, med):
    """
    Crop dendritic pixels from an ROI by finding the soma boundary.

    Computes cumulative weighted area as a function of distance from the median
    center, then finds the radius where the area growth drops below a threshold.

    Parameters
    ----------
    ypix : numpy.ndarray
        Y-coordinates of the ROI pixels.
    xpix : numpy.ndarray
        X-coordinates of the ROI pixels.
    lam : numpy.ndarray
        Weights (e.g. fluorescence) for each pixel.
    med : list of float
        Two-element list [ymed, xmed] of the ROI center.

    Returns
    -------
    crop : numpy.ndarray
        Boolean array of length len(ypix), True for pixels within the soma.
    """
    crop = np.ones(ypix.size, "bool")
    if len(ypix) > 10:
        dists = ((ypix - med[0])**2 + (xpix - med[1])**2)**0.5
        radii = np.arange(0, dists.max(), 1)
        area = np.zeros_like(radii)
        for k, radius in enumerate(radii):
            area[k] = lam[dists < radius].sum()
        darea = np.diff(area)
        radius = radii[-1]
        threshold = darea.max() / 3
        if len(np.nonzero(darea > threshold)[0]) > 0:
            ida = np.nonzero(darea > threshold)[0][0]
            if len(np.nonzero(darea[ida:] < threshold)[0]):
                radius = radii[np.nonzero(darea[ida:] < threshold)[0][0] + ida]
        crop = dists < radius
    if crop.sum() == 0:
        crop = np.ones(ypix.size, "bool")
    return crop

Copyright © 2023 Howard Hughes Medical Institute, Authored by Carsen Stringer and Marius Pachitariu.

circleMask #

circleMask(d0)

Create a normalized distance array and return indices within a unit circle.

Parameters:

Name Type Description Default
d0 list of float

Two-element list [d0_y, d0_x] giving the radius in each dimension.

 required

Returns:

Name Type Description
rs ndarray

Normalized distance array of shape (2d0_y+1, 2d0_x+1).
dx ndarray

Normalized X-coordinates of pixels within the unit circle.
dy ndarray

Normalized Y-coordinates of pixels within the unit circle.
Source code in suite2p/detection/utils.py 
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
def circleMask(d0):
    """
    Create a normalized distance array and return indices within a unit circle.

    Parameters
    ----------
    d0 : list of float
        Two-element list [d0_y, d0_x] giving the radius in each dimension.

    Returns
    -------
    rs : numpy.ndarray
        Normalized distance array of shape (2*d0_y+1, 2*d0_x+1).
    dx : numpy.ndarray
        Normalized X-coordinates of pixels within the unit circle.
    dy : numpy.ndarray
        Normalized Y-coordinates of pixels within the unit circle.
    """
    d00 = int(np.round(d0[0]))
    d01 = int(np.round(d0[1]))
    dy = np.tile(np.arange(-d00, d00 + 1) / d00, (2 * d01 + 1, 1))
    dy = dy.transpose()
    dx = np.tile(np.arange(-d01, d01 + 1) / d01, (2 * d00 + 1, 1))

    rs = (dy**2 + dx**2)**0.5
    dx = dx[rs <= 1.]
    dy = dy[rs <= 1.]
    return rs, dx, dy

hp_gaussian_filter #

hp_gaussian_filter(mov, width)

Returns a high-pass-filtered mov by subtracting off the movie smoothed by a gaussian kernel.

Parameters:

Name Type Description Default
mov ndarray

Movie of shape (nframes, Ly, Lx).

 required
width int

The standard deviation of the Gaussian filter in time

 required

Returns:

Name Type Description
filtered_mov nImg x Ly x Lx

The filtered video
Source code in suite2p/detection/utils.py 
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
def hp_gaussian_filter(mov: np.ndarray, width: int) -> np.ndarray:
    """
    Returns a high-pass-filtered `mov` by subtracting off the movie smoothed by a gaussian kernel.

    Parameters
    ----------
    mov: numpy.ndarray
        Movie of shape (nframes, Ly, Lx).
    width: int
        The standard deviation of the Gaussian filter in time

    Returns
    -------
    filtered_mov: nImg x Ly x Lx
        The filtered video
    """
    mov = mov.copy()
    for j in range(mov.shape[1]):
        mov[:, j, :] -= gaussian_filter1d(mov[:, j, :], width, axis=0)
    return mov

hp_rolling_mean_filter #

hp_rolling_mean_filter(mov, width)

Returns high-pass-filtered mov by subtracting off the rolling mean in window of width.

Parameters:

Name Type Description Default
mov ndarray

Movie of shape (nframes, Ly, Lx).

 required
width int

The filter width in time.

 required

Returns:

Name Type Description
filtered_mov ndarray

Movie of shape (nframes, Ly, Lx), high-pass filtered in time.
Source code in suite2p/detection/utils.py 
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
def hp_rolling_mean_filter(mov: np.ndarray, width: int) -> np.ndarray:
    """
    Returns high-pass-filtered `mov` by subtracting off the rolling mean in window of `width`.

    Parameters
    ----------
    mov: numpy.ndarray
        Movie of shape (nframes, Ly, Lx).
    width: int
        The filter width in time.

    Returns
    -------
    filtered_mov: numpy.ndarray
        Movie of shape (nframes, Ly, Lx), high-pass filtered in time.

    """
    mov = mov.copy()
    for i in range(0, mov.shape[0], width):
        mov[i:i + width, :, :] -= mov[i:i + width, :, :].mean(axis=0)
    return mov

square_mask #

square_mask(mask, ly, yi, xi)

Crop a square patch from a 2D mask centered at a given position.

Parameters:

Name Type Description Default
mask ndarray

2D array to crop from, shape (Lyc, Lxc).

 required
ly int

Half-width of the square patch (output is 2ly x 2ly).

 required
yi int

Y-coordinate of the center.

 required
xi int

X-coordinate of the center.

 required

Returns:

Name Type Description
mask0 ndarray

Cropped square patch of shape (2ly, 2ly).
Source code in suite2p/detection/utils.py 
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
def square_mask(mask, ly, yi, xi):
    """
    Crop a square patch from a 2D mask centered at a given position.

    Parameters
    ----------
    mask : numpy.ndarray
        2D array to crop from, shape (Lyc, Lxc).
    ly : int
        Half-width of the square patch (output is 2*ly x 2*ly).
    yi : int
        Y-coordinate of the center.
    xi : int
        X-coordinate of the center.

    Returns
    -------
    mask0 : numpy.ndarray
        Cropped square patch of shape (2*ly, 2*ly).
    """
    Lyc, Lxc = mask.shape
    mask0 = np.zeros((2 * ly, 2 * ly), mask.dtype)
    yinds = [max(0, yi - ly), min(yi + ly, Lyc)]
    xinds = [max(0, xi - ly), min(xi + ly, Lxc)]
    mask0[max(0, ly - yi):min(2 * ly, Lyc + ly - yi),
          max(0, ly - xi):min(2 * ly, Lxc + ly - xi)] = mask[yinds[0]:yinds[1],
                                                             xinds[0]:xinds[1]]
    return mask0

standard_deviation_over_time #

standard_deviation_over_time(mov, batch_size)

Returns standard deviation of difference between pixels across time, computed in batches of batch_size.

Parameters:

Name Type Description Default
mov ndarray

Movie of shape (nframes, Ly, Lx).

 required
batch_size int

The batch size for computing the standard deviation of the difference.

 required

Returns:

Name Type Description
sdmov Ly x Lx

The standard deviation of the difference across time for each pixel.
Source code in suite2p/detection/utils.py 
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
def standard_deviation_over_time(mov: np.ndarray, batch_size: int) -> np.ndarray:
    """
    Returns standard deviation of difference between pixels across time, computed in batches of batch_size.

    Parameters
    ----------
    mov: numpy.ndarray
        Movie of shape (nframes, Ly, Lx).
    batch_size: int
        The batch size for computing the standard deviation of the difference.

    Returns
    -------
    sdmov: Ly x Lx
        The standard deviation of the difference across time for each pixel.
    """
    nbins, Ly, Lx = mov.shape
    batch_size = min(batch_size, nbins)
    sdmov = np.zeros((Ly, Lx), "float32")
    for ix in range(0, nbins, batch_size):
        sdmov += ((np.diff(mov[ix:ix + batch_size, :, :], axis=0)**2).sum(axis=0))
    sdmov = np.maximum(1e-10, np.sqrt(sdmov / nbins))
    return sdmov

temporal_high_pass_filter #

temporal_high_pass_filter(mov, width)

Returns hp-filtered mov over time, selecting an algorithm for computational performance based on the kernel width.

Parameters:

Name Type Description Default
mov ndarray

Movie of shape (nframes, Ly, Lx).

 required
width int

The filter width in time.

 required

Returns:

Name Type Description
filtered_mov ndarray

Movie of shape (nframes, Ly, Lx), high-pass filtered in time.
Source code in suite2p/detection/utils.py 
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
def temporal_high_pass_filter(mov: np.ndarray, width: int) -> np.ndarray:
    """
    Returns hp-filtered mov over time, selecting an algorithm for computational performance based on the kernel width.

    Parameters
    ----------
    mov: numpy.ndarray
        Movie of shape (nframes, Ly, Lx).
    width: int
        The filter width in time.

    Returns
    -------
    filtered_mov: numpy.ndarray
        Movie of shape (nframes, Ly, Lx), high-pass filtered in time.

    """

    return hp_gaussian_filter(mov, width) if width < 10 else hp_rolling_mean_filter(
        mov, width)  # gaussian is slower