suite2p.extraction package

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

baseline_maximin

baseline_maximin(F, win_baseline, sig_baseline, fs, batch_size=100, device=torch.device('cuda'))

Compute maximin baseline using GPU-accelerated max/min pooling.

Smooths traces with a Gaussian filter, then applies a rolling minimum followed by a rolling maximum to estimate the baseline.

Parameters:

Name	Type	Description	Default
F	ndarray	Fluorescence traces, shape (n_neurons, n_frames).	required
win_baseline	float	Window in seconds for the max/min filters.	required
sig_baseline	float	Standard deviation of the Gaussian filter in frames.	required
fs	float	Sampling rate per plane in Hz.	required
batch_size	int, optional (default 100)	Number of neurons processed per batch.	100
device	torch.device, optional (default torch.device("cuda"))	Torch device for performing operations.	device('cuda')

Returns:

Name	Type	Description
F	ndarray	Baseline-corrected fluorescence, shape (n_neurons, n_frames).

Source code in suite2p/extraction/dcnv.py

def baseline_maximin(F, win_baseline, sig_baseline,
                     fs, batch_size=100, device=torch.device('cuda')):
    """
    Compute maximin baseline using GPU-accelerated max/min pooling.

    Smooths traces with a Gaussian filter, then applies a rolling minimum
    followed by a rolling maximum to estimate the baseline.

    Parameters
    ----------
    F : numpy.ndarray
        Fluorescence traces, shape (n_neurons, n_frames).
    win_baseline : float
        Window in seconds for the max/min filters.
    sig_baseline : float
        Standard deviation of the Gaussian filter in frames.
    fs : float
        Sampling rate per plane in Hz.
    batch_size : int, optional (default 100)
        Number of neurons processed per batch.
    device : torch.device, optional (default torch.device("cuda"))
        Torch device for performing operations.

    Returns
    -------
    F : numpy.ndarray
        Baseline-corrected fluorescence, shape (n_neurons, n_frames).
    """
    win = int(win_baseline * fs)
    win += 1 if win%2==0 else 0
    ncells, n_frames = F.shape
    Flow = np.zeros((ncells, n_frames), 'float32')
    batch_size = int(batch_size)
    n_batches = int(np.ceil(ncells / batch_size))
    gwid = int(np.round(sig_baseline * 3))
    gaussian = torch.exp(- torch.arange(-gwid, gwid + 1, 1, device=device)**2 /
                          (2 * sig_baseline**2))
    gaussian /= gaussian.sum()

    for n in range(n_batches):
        nstart, nend = n * batch_size, min((n+1) * batch_size, ncells)
        data = torch.from_numpy(F[nstart : nend]).to(device, dtype=torch.float)      

        data = pad(data, (gwid, gwid), 'replicate')
        data = conv1d(data.unsqueeze(1), gaussian.unsqueeze(0).unsqueeze(0), padding=0)

        data = pad(data, (win//2, win//2), 'replicate')
        data = -max_pool1d(-data, kernel_size=win, stride=1, padding= 0)

        data = pad(data, (win//2, win//2), 'replicate')
        data = max_pool1d(data, kernel_size=win, stride=1, padding= 0)

        Flow[nstart : nend] = data.squeeze().cpu().numpy()

    F = F - Flow

    return F 

oasis

oasis(F, batch_size, tau, fs)

Compute non-negative deconvolution with no sparsity constraints.

Uses OASIS algorithm (Friedrich, Zhou & Paninski, 2017).

Parameters:

Name	Type	Description	Default
F	ndarray	Fluorescence traces, shape (n_neurons, n_frames). In the pipeline, uses neuropil-subtracted and baselined fluorescence.	required
batch_size	int	Number of neurons processed per batch.	required
tau	float	Timescale of the indicator decay in seconds.	required
fs	float	Sampling rate per plane in Hz.	required

Returns:

Name	Type	Description
S	ndarray	Deconvolved fluorescence, shape (n_neurons, n_frames).

Source code in suite2p/extraction/dcnv.py

def oasis(F, batch_size, tau, fs):
    """
    Compute non-negative deconvolution with no sparsity constraints.

    Uses OASIS algorithm (Friedrich, Zhou & Paninski, 2017).

    Parameters
    ----------
    F : numpy.ndarray
        Fluorescence traces, shape (n_neurons, n_frames). In the pipeline,
        uses neuropil-subtracted and baselined fluorescence.
    batch_size : int
        Number of neurons processed per batch.
    tau : float
        Timescale of the indicator decay in seconds.
    fs : float
        Sampling rate per plane in Hz.

    Returns
    -------
    S : numpy.ndarray
        Deconvolved fluorescence, shape (n_neurons, n_frames).
    """

    NN, NT = F.shape
    F = F.astype(np.float32)
    S = np.zeros((NN, NT), dtype=np.float32)
    n_batches = int(np.ceil(NN / batch_size))
    logger.info(f"Deconvolving {NN} neurons in {n_batches} batches")
    tqdm_out = TqdmToLogger(logger, level=logging.INFO)
    for n in trange(n_batches, file=tqdm_out):
        i = n * batch_size
        f = F[i:i + batch_size]
        v = np.zeros((f.shape[0], NT), dtype=np.float32)
        w = np.zeros((f.shape[0], NT), dtype=np.float32)
        t = np.zeros((f.shape[0], NT), dtype=np.int64)
        l = np.zeros((f.shape[0], NT), dtype=np.float32)
        s = np.zeros((f.shape[0], NT), dtype=np.float32)
        oasis_matrix(f, v, w, t, l, s, tau, fs)
        S[i:i + batch_size] = s
    return S

oasis_matrix

oasis_matrix(F, v, w, t, l, s, tau, fs)

Spike deconvolution on many neurons parallelized with prange.

Parameters:

Name	Type	Description	Default
F	ndarray	Fluorescence traces, shape (n_neurons, n_frames).	required
v	ndarray	Pool values buffer, shape (n_neurons, n_frames).	required
w	ndarray	Pool weights buffer, shape (n_neurons, n_frames).	required
t	ndarray	Pool start times buffer, shape (n_neurons, n_frames).	required
l	ndarray	Pool lengths buffer, shape (n_neurons, n_frames).	required
s	ndarray	Output spike traces, shape (n_neurons, n_frames). Modified in place.	required
tau	float	Timescale of the indicator decay in seconds.	required
fs	float	Sampling rate per plane in Hz.	required

Source code in suite2p/extraction/dcnv.py

@njit([
    "float32[:,:], float32[:,:], float32[:,:], int64[:,:], float32[:,:], float32[:,:], float32, float32"
], parallel=True, cache=True)
def oasis_matrix(F, v, w, t, l, s, tau, fs):
    """
    Spike deconvolution on many neurons parallelized with prange.

    Parameters
    ----------
    F : numpy.ndarray
        Fluorescence traces, shape (n_neurons, n_frames).
    v : numpy.ndarray
        Pool values buffer, shape (n_neurons, n_frames).
    w : numpy.ndarray
        Pool weights buffer, shape (n_neurons, n_frames).
    t : numpy.ndarray
        Pool start times buffer, shape (n_neurons, n_frames).
    l : numpy.ndarray
        Pool lengths buffer, shape (n_neurons, n_frames).
    s : numpy.ndarray
        Output spike traces, shape (n_neurons, n_frames). Modified in place.
    tau : float
        Timescale of the indicator decay in seconds.
    fs : float
        Sampling rate per plane in Hz.
    """
    for n in prange(F.shape[0]):
        oasis_trace(F[n], v[n], w[n], t[n], l[n], s[n], tau, fs)

oasis_trace

oasis_trace(F, v, w, t, l, s, tau, fs)

Spike deconvolution on a single neuron using the OASIS algorithm.

Parameters:

Name	Type	Description	Default
F	ndarray	Fluorescence trace, shape (n_frames,).	required
v	ndarray	Pool values buffer, shape (n_frames,).	required
w	ndarray	Pool weights buffer, shape (n_frames,).	required
t	ndarray	Pool start times buffer, shape (n_frames,).	required
l	ndarray	Pool lengths buffer, shape (n_frames,).	required
s	ndarray	Output spike trace, shape (n_frames,). Modified in place.	required
tau	float	Timescale of the indicator decay in seconds.	required
fs	float	Sampling rate per plane in Hz.	required

Source code in suite2p/extraction/dcnv.py

@njit([
    "float32[:], float32[:], float32[:], int64[:], float32[:], float32[:], float32, float32"
], cache=True)
def oasis_trace(F, v, w, t, l, s, tau, fs):
    """
    Spike deconvolution on a single neuron using the OASIS algorithm.

    Parameters
    ----------
    F : numpy.ndarray
        Fluorescence trace, shape (n_frames,).
    v : numpy.ndarray
        Pool values buffer, shape (n_frames,).
    w : numpy.ndarray
        Pool weights buffer, shape (n_frames,).
    t : numpy.ndarray
        Pool start times buffer, shape (n_frames,).
    l : numpy.ndarray
        Pool lengths buffer, shape (n_frames,).
    s : numpy.ndarray
        Output spike trace, shape (n_frames,). Modified in place.
    tau : float
        Timescale of the indicator decay in seconds.
    fs : float
        Sampling rate per plane in Hz.
    """
    NT = F.shape[0]
    g = -1. / (tau * fs)

    it = 0
    ip = 0

    while it < NT:
        v[ip], w[ip], t[ip], l[ip] = F[it], 1, it, 1
        while ip > 0:
            if v[ip - 1] * np.exp(g * l[ip - 1]) > v[ip]:
                # violation of the constraint means merging pools
                f1 = np.exp(g * l[ip - 1])
                f2 = np.exp(2 * g * l[ip - 1])
                wnew = w[ip - 1] + w[ip] * f2
                v[ip - 1] = (v[ip - 1] * w[ip - 1] + v[ip] * w[ip] * f1) / wnew
                w[ip - 1] = wnew
                l[ip - 1] = l[ip - 1] + l[ip]
                ip -= 1
            else:
                break
        it += 1
        ip += 1

    s[t[1:ip]] = v[1:ip] - v[:ip - 1] * np.exp(g * l[:ip - 1])

preprocess

preprocess(F, baseline, win_baseline, sig_baseline, fs, prctile_baseline=8, batch_size=100, device=torch.device('cuda'))

Preprocess fluorescence traces for spike deconvolution.

Subtracts a baseline estimated with a rolling maximin filter, constant minimum, or percentile, depending on the baseline setting.

Parameters:

Name	Type	Description	Default
F	ndarray	Fluorescence traces, shape (n_neurons, n_frames). In the pipeline, uses neuropil-subtracted fluorescence.	required
baseline	str	Method for computing the baseline of each trace. One of "maximin", "constant", or "prctile".	required
win_baseline	float	Window in seconds for the max/min filters.	required
sig_baseline	float	Standard deviation of the Gaussian filter in frames.	required
fs	float	Sampling rate per plane in Hz.	required
prctile_baseline	float, optional (default 8)	Percentile of trace to use as baseline when baseline is "prctile".	8

Returns:

Name	Type	Description
F	ndarray	Baseline-corrected fluorescence, shape (n_neurons, n_frames).

Source code in suite2p/extraction/dcnv.py

def preprocess(F, baseline, win_baseline, sig_baseline,
               fs, prctile_baseline=8, batch_size=100, 
               device=torch.device('cuda')):
    """
    Preprocess fluorescence traces for spike deconvolution.

    Subtracts a baseline estimated with a rolling maximin filter, constant
    minimum, or percentile, depending on the baseline setting.

    Parameters
    ----------
    F : numpy.ndarray
        Fluorescence traces, shape (n_neurons, n_frames). In the pipeline,
        uses neuropil-subtracted fluorescence.
    baseline : str
        Method for computing the baseline of each trace. One of "maximin",
        "constant", or "prctile".
    win_baseline : float
        Window in seconds for the max/min filters.
    sig_baseline : float
        Standard deviation of the Gaussian filter in frames.
    fs : float
        Sampling rate per plane in Hz.
    prctile_baseline : float, optional (default 8)
        Percentile of trace to use as baseline when baseline is "prctile".

    Returns
    -------
    F : numpy.ndarray
        Baseline-corrected fluorescence, shape (n_neurons, n_frames).
    """
    win = int(win_baseline * fs)
    if baseline == "maximin":
        # Flow = gaussian_filter(F, [0., sig_baseline])
        # Flow = minimum_filter1d(Flow, win)
        # Flow = maximum_filter1d(Flow, win)
        F = baseline_maximin(F, win_baseline, sig_baseline, fs, 
                             batch_size=batch_size, device=device)
    elif baseline == "constant":
        Flow = gaussian_filter(F, [0., sig_baseline])
        Flow = np.amin(Flow)
        F -= Flow
    elif baseline == "prctile":
        Flow = np.percentile(F, prctile_baseline, axis=1)
        Flow = np.expand_dims(Flow, axis=1)
        F -= Flow
    else:
        Flow = 0.

    return F

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

extract_traces

extract_traces(f_in, cell_masks, neuropil_masks, batch_size=500, device=torch.device('cuda'))

Extract fluorescence traces using cell and neuropil masks.

Performs sparse matrix multiplication on the GPU to extract ROI and neuropil traces from registered frames.

Parameters:

Name	Type	Description	Default
f_in	ndarray or BinaryFile	Registered frames, shape (n_frames, Ly, Lx).	required
cell_masks	list of tuple	Each element is a tuple (pixel_indices, weights) for one ROI, where pixel_indices are flattened and weights are normalized to sum to 1.	required
neuropil_masks	list of numpy.ndarray	Each element is an array of flattened pixel indices for the neuropil mask of one ROI.	required
batch_size	int, optional (default 500)	Number of frames processed per batch.	500
device	torch.device, optional (default torch.device("cuda"))	Torch device for performing operations.	device('cuda')

Returns:

Name	Type	Description
F	ndarray	ROI fluorescence traces, shape (n_rois, n_frames).
Fneu	ndarray	Neuropil fluorescence traces, shape (n_rois, n_frames).

Source code in suite2p/extraction/extract.py

def extract_traces(f_in, cell_masks, neuropil_masks, batch_size=500, 
                    device = torch.device("cuda")):
    """
    Extract fluorescence traces using cell and neuropil masks.

    Performs sparse matrix multiplication on the GPU to extract ROI and
    neuropil traces from registered frames.

    Parameters
    ----------
    f_in : numpy.ndarray or BinaryFile
        Registered frames, shape (n_frames, Ly, Lx).
    cell_masks : list of tuple
        Each element is a tuple (pixel_indices, weights) for one ROI, where
        pixel_indices are flattened and weights are normalized to sum to 1.
    neuropil_masks : list of numpy.ndarray
        Each element is an array of flattened pixel indices for the neuropil
        mask of one ROI.
    batch_size : int, optional (default 500)
        Number of frames processed per batch.
    device : torch.device, optional (default torch.device("cuda"))
        Torch device for performing operations.

    Returns
    -------
    F : numpy.ndarray
        ROI fluorescence traces, shape (n_rois, n_frames).
    Fneu : numpy.ndarray
        Neuropil fluorescence traces, shape (n_rois, n_frames).
    """
    n_frames, Ly, Lx = f_in.shape
    t0 = time.time()
    batch_size = min(batch_size, 1000)
    ncells = len(cell_masks)

    # Check if mps and force to be on the CPU so that extraction can be run
    # TODO: Once, sparse_coo_tensor works on sparseMPS backend, we should remove this check
    if device.type == 'mps':
        device = torch.device('cpu')

    if neuropil_masks is not None:
        npix_neuropil = torch.Tensor([len(nm) for nm in neuropil_masks]).to(device)
        # create coo tensor of neuropil masks
        ccol_indices = [m for nm in neuropil_masks for m in nm]
        row_indices = [k for k in range(len(neuropil_masks)) for m in neuropil_masks[k]]
        inds = torch.Tensor([ccol_indices, row_indices]).to(device)
        # convert to csc (tried creating csc directly but it was slow)
        nmasks = torch.sparse_coo_tensor(inds, torch.ones(len(row_indices), device=device),
                                         size=(Ly*Lx, ncells))
        nmasks = nmasks.to_sparse_csc()

    ccol_indices = [m for cm in cell_masks for m in cm[0]]
    row_indices = [k for k in range(len(cell_masks)) for m in cell_masks[k][0]]
    cell_lam = torch.Tensor([l for cm in cell_masks for l in cm[1]]).to(device)
    inds = torch.Tensor([ccol_indices, row_indices]).to(device)
    cmasks = torch.sparse_coo_tensor(inds, cell_lam,
                                     size=(Ly*Lx, ncells))
    cmasks = cmasks.to_sparse_csc()

    F = np.zeros((ncells, n_frames), np.float32)
    Fneu = np.zeros((ncells, n_frames), np.float32)

    batch_size = int(batch_size)

    n_batches = int(np.ceil(n_frames / batch_size))
    logger.info(f"Extracting fluorescence from {n_frames} frames in {n_batches} batches")
    tqdm_out = TqdmToLogger(logger, level=logging.INFO)
    for n in trange(n_batches, mininterval=10, file=tqdm_out):
        tstart, tend = n * batch_size, min((n+1) * batch_size, n_frames)
        data = torch.from_numpy(f_in[tstart : tend]).to(device)
        data = data.reshape(-1, Ly*Lx).float()

        if neuropil_masks is not None:
            Fneu_batch = (data @ nmasks) / npix_neuropil
            Fneu[:, tstart : tend] = Fneu_batch.T.cpu().numpy()

        F_batch = data @ cmasks
        F[:, tstart : tend] = F_batch.T.cpu().numpy()

    return F, Fneu

extraction_wrapper

extraction_wrapper(stat, f_reg, f_reg_chan2=None, cell_masks=None, neuropil_masks=None, settings=default_settings()['extraction'], device=torch.device('cuda'))

Main fluorescence extraction function.

Creates cell and neuropil masks (if not provided) and extracts fluorescence traces for both functional and anatomical channels.

Parameters:

Name	Type	Description	Default
stat	ndarray	Array of ROI statistics dictionaries, each containing "ypix", "xpix", "lam", "radius", and "overlap".	required
f_reg	ndarray or BinaryFile	Registered functional frames, shape (n_frames, Ly, Lx).	required
f_reg_chan2	numpy.ndarray or BinaryFile, optional (default None)	Registered anatomical channel frames, shape (n_frames, Ly, Lx).	None
cell_masks	list of tuple, optional (default None)	Pre-computed cell masks. If None, masks are created from stat.	None
neuropil_masks	list of numpy.ndarray, optional (default None)	Pre-computed neuropil masks. If None, masks are created from stat.	None
settings	dict	Extraction settings dictionary.	default_settings()['extraction']
device	torch.device, optional (default torch.device("cuda"))	Torch device for performing operations.	device('cuda')

Returns:

Name	Type	Description
F	ndarray	ROI fluorescence traces for functional channel, shape (n_rois, n_frames).
Fneu	ndarray	Neuropil fluorescence traces for functional channel, shape (n_rois, n_frames).
F_chan2	ndarray or None	ROI fluorescence traces for anatomical channel, or None if f_reg_chan2 is not provided.
Fneu_chan2	ndarray or None	Neuropil fluorescence traces for anatomical channel, or None if f_reg_chan2 is not provided.

Source code in suite2p/extraction/extract.py

def extraction_wrapper(stat, f_reg, f_reg_chan2=None, cell_masks=None,
                       neuropil_masks=None, settings=default_settings()["extraction"],
                        device = torch.device("cuda")):
    """
    Main fluorescence extraction function.

    Creates cell and neuropil masks (if not provided) and extracts
    fluorescence traces for both functional and anatomical channels.

    Parameters
    ----------
    stat : numpy.ndarray
        Array of ROI statistics dictionaries, each containing "ypix", "xpix",
        "lam", "radius", and "overlap".
    f_reg : numpy.ndarray or BinaryFile
        Registered functional frames, shape (n_frames, Ly, Lx).
    f_reg_chan2 : numpy.ndarray or BinaryFile, optional (default None)
        Registered anatomical channel frames, shape (n_frames, Ly, Lx).
    cell_masks : list of tuple, optional (default None)
        Pre-computed cell masks. If None, masks are created from stat.
    neuropil_masks : list of numpy.ndarray, optional (default None)
        Pre-computed neuropil masks. If None, masks are created from stat.
    settings : dict
        Extraction settings dictionary.
    device : torch.device, optional (default torch.device("cuda"))
        Torch device for performing operations.

    Returns
    -------
    F : numpy.ndarray
        ROI fluorescence traces for functional channel, shape (n_rois, n_frames).
    Fneu : numpy.ndarray
        Neuropil fluorescence traces for functional channel, shape (n_rois, n_frames).
    F_chan2 : numpy.ndarray or None
        ROI fluorescence traces for anatomical channel, or None if f_reg_chan2
        is not provided.
    Fneu_chan2 : numpy.ndarray or None
        Neuropil fluorescence traces for anatomical channel, or None if
        f_reg_chan2 is not provided.
    """
    n_frames, Ly, Lx = f_reg.shape
    batch_size = settings["batch_size"]
    if cell_masks is None:
        t10 = time.time()
        cell_masks, neuropil_masks0 = create_masks(stat, Ly, Lx, 
                                                   lam_percentile=settings["lam_percentile"], 
                                                allow_overlap=settings["allow_overlap"], 
                                                neuropil_extract=settings["neuropil_extract"], 
                                                inner_neuropil_radius=settings["inner_neuropil_radius"],
                                                min_neuropil_pixels=settings["min_neuropil_pixels"], 
                                                circular_neuropil=settings["circular_neuropil"])
        if neuropil_masks is None:
            neuropil_masks = neuropil_masks0
        logger.info("Masks created, %0.2f sec." % (time.time() - t10))

    t0 = time.time()
    ncells = len(cell_masks)
    logger.info("functional channel:")
    F, Fneu = extract_traces(f_reg, cell_masks, neuropil_masks, 
                             batch_size=batch_size, device=device)
    if f_reg_chan2 is not None:
        logger.info("anatomical channel:")
        F_chan2, Fneu_chan2 = extract_traces(f_reg_chan2, cell_masks, neuropil_masks,
                                             batch_size=batch_size, device=device)
    else:
        F_chan2, Fneu_chan2 = None, None

    logger.info("Extracted fluorescence from %d ROIs in %d frames, %0.2f sec." %
          (ncells, n_frames, time.time() - t0))


    return F, Fneu, F_chan2, Fneu_chan2

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

create_cell_mask

create_cell_mask(stat, Ly, Lx, allow_overlap=False)

Create the cell mask for a single ROI.

Parameters:

Name	Type	Description	Default
stat	dict	ROI statistics dictionary containing "ypix", "xpix", "lam", and "overlap".	required
Ly	int	Height of the image in pixels.	required
Lx	int	Width of the image in pixels.	required
allow_overlap	bool, optional (default False)	If True, include overlapping pixels in the cell mask.	False

Returns:

Name	Type	Description
cell_mask	ndarray	Flattened pixel indices for the ROI.
lam_normed	ndarray	Normalized pixel weights summing to 1.

Source code in suite2p/extraction/masks.py

def create_cell_mask(stat, Ly, Lx, allow_overlap = False):
    """
    Create the cell mask for a single ROI.

    Parameters
    ----------
    stat : dict
        ROI statistics dictionary containing "ypix", "xpix", "lam", and
        "overlap".
    Ly : int
        Height of the image in pixels.
    Lx : int
        Width of the image in pixels.
    allow_overlap : bool, optional (default False)
        If True, include overlapping pixels in the cell mask.

    Returns
    -------
    cell_mask : numpy.ndarray
        Flattened pixel indices for the ROI.
    lam_normed : numpy.ndarray
        Normalized pixel weights summing to 1.
    """
    mask = ... if allow_overlap else ~stat["overlap"]
    cell_mask = np.ravel_multi_index((stat["ypix"], stat["xpix"]), (Ly, Lx))
    cell_mask = cell_mask[mask]
    lam = stat["lam"][mask]
    lam_normed = lam / lam.sum() if lam.size > 0 else np.empty(0)
    return cell_mask, lam_normed

create_cell_pix

create_cell_pix(stats, Ly, Lx, lam_percentile=50.0)

Create a binary image indicating which pixels contain a cell.

Pixels with low cell weights can be excluded using lam_percentile; set lam_percentile to 0 to disable filtering.

Parameters:

Name	Type	Description	Default
stats	list of dict	List of ROI statistics dictionaries, each containing "ypix", "xpix", "lam", and "radius".	required
Ly	int	Height of the image in pixels.	required
Lx	int	Width of the image in pixels.	required
lam_percentile	float, optional (default 50.0)	Percentile filter threshold for excluding low-weight pixels.	50.0

Returns:

Name	Type	Description
cell_pix	ndarray	Boolean array of shape (Ly, Lx), True where a cell pixel exists.

Source code in suite2p/extraction/masks.py

def create_cell_pix(stats, Ly, Lx, lam_percentile = 50.0):
    """
    Create a binary image indicating which pixels contain a cell.

    Pixels with low cell weights can be excluded using lam_percentile;
    set lam_percentile to 0 to disable filtering.

    Parameters
    ----------
    stats : list of dict
        List of ROI statistics dictionaries, each containing "ypix", "xpix",
        "lam", and "radius".
    Ly : int
        Height of the image in pixels.
    Lx : int
        Width of the image in pixels.
    lam_percentile : float, optional (default 50.0)
        Percentile filter threshold for excluding low-weight pixels.

    Returns
    -------
    cell_pix : numpy.ndarray
        Boolean array of shape (Ly, Lx), True where a cell pixel exists.
    """
    cell_pix = np.zeros((Ly, Lx))
    lammap = np.zeros((Ly, Lx))
    radii = np.zeros(len(stats))
    for ni, stat in enumerate(stats):
        radii[ni] = stat["radius"]
        ypix = stat["ypix"]
        xpix = stat["xpix"]
        lam = stat["lam"]
        lammap[ypix, xpix] = np.maximum(lammap[ypix, xpix], lam)
    radius = np.median(radii)
    if lam_percentile > 0.0:
        filt = percentile_filter(lammap, percentile=lam_percentile,
                                 size=int(radius * 5))
        cell_pix = ~np.logical_or(lammap < filt, lammap == 0)
    else:
        cell_pix = lammap > 0.0

    return cell_pix

create_masks

create_masks(stats, Ly, Lx, lam_percentile=50.0, allow_overlap=False, neuropil_extract=True, inner_neuropil_radius=2, min_neuropil_pixels=350, circular_neuropil=False)

Create cell and neuropil masks for all ROIs.

Parameters:

Name	Type	Description	Default
stats	list of dict	List of ROI statistics dictionaries, each containing "ypix", "xpix", "lam", "radius", and "overlap".	required
Ly	int	Height of the image in pixels.	required
Lx	int	Width of the image in pixels.	required
lam_percentile	float, optional (default 50.)	Percentile filter threshold for excluding low-weight cell pixels from neuropil masks. Set to 0 to disable.	50.0
allow_overlap	bool, optional (default False)	If True, include overlapping pixels in cell masks.	False
neuropil_extract	bool, optional (default True)	If True, compute neuropil masks.	True
inner_neuropil_radius	int, optional (default 2)	Number of pixels to extend around each ROI as an exclusion zone for the neuropil mask.	2
min_neuropil_pixels	int, optional (default 350)	Minimum number of pixels in the neuropil mask.	350
circular_neuropil	bool, optional (default False)	If True, extend neuropil masks circularly rather than as rectangles (SLOW).	False

Returns:

Name	Type	Description
cell_masks	list of tuple	Each element is a tuple (pixel_indices, weights) for one ROI.
neuropil_masks	list of numpy.ndarray or None	Each element is an array of flattened pixel indices for the neuropil mask. None if neuropil_extract is False.

Source code in suite2p/extraction/masks.py

def create_masks(stats, Ly, Lx, lam_percentile=50., 
                 allow_overlap=False, neuropil_extract=True, inner_neuropil_radius=2,
                 min_neuropil_pixels=350, circular_neuropil=False):
    """
    Create cell and neuropil masks for all ROIs.

    Parameters
    ----------
    stats : list of dict
        List of ROI statistics dictionaries, each containing "ypix", "xpix",
        "lam", "radius", and "overlap".
    Ly : int
        Height of the image in pixels.
    Lx : int
        Width of the image in pixels.
    lam_percentile : float, optional (default 50.)
        Percentile filter threshold for excluding low-weight cell pixels from
        neuropil masks. Set to 0 to disable.
    allow_overlap : bool, optional (default False)
        If True, include overlapping pixels in cell masks.
    neuropil_extract : bool, optional (default True)
        If True, compute neuropil masks.
    inner_neuropil_radius : int, optional (default 2)
        Number of pixels to extend around each ROI as an exclusion zone for
        the neuropil mask.
    min_neuropil_pixels : int, optional (default 350)
        Minimum number of pixels in the neuropil mask.
    circular_neuropil : bool, optional (default False)
        If True, extend neuropil masks circularly rather than as rectangles (SLOW).

    Returns
    -------
    cell_masks : list of tuple
        Each element is a tuple (pixel_indices, weights) for one ROI.
    neuropil_masks : list of numpy.ndarray or None
        Each element is an array of flattened pixel indices for the neuropil
        mask. None if neuropil_extract is False.
    """

    cell_pix = create_cell_pix(stats, Ly=Ly, Lx=Lx,
                               lam_percentile=lam_percentile)
    cell_masks = [
        create_cell_mask(stat, Ly=Ly, Lx=Lx, allow_overlap=allow_overlap)
        for stat in stats
    ]
    if neuropil_extract:
        neuropil_masks = create_neuropil_masks(
            ypixs=[stat["ypix"] for stat in stats],
            xpixs=[stat["xpix"] for stat in stats], cell_pix=cell_pix,
            inner_neuropil_radius=inner_neuropil_radius,
            min_neuropil_pixels=min_neuropil_pixels,
            circular=circular_neuropil)
    else:
        neuropil_masks = None
    return cell_masks, neuropil_masks

create_neuropil_masks

create_neuropil_masks(ypixs, xpixs, cell_pix, inner_neuropil_radius=2, min_neuropil_pixels=350, circular=False)

Create surround neuropil masks for ROIs by extending each ROI outward.

Parameters:

Name	Type	Description	Default
ypixs	list of numpy.ndarray	Y-coordinates of the pixels for each ROI.	required
xpixs	list of numpy.ndarray	X-coordinates of the pixels for each ROI.	required
cell_pix	ndarray	Binary array of shape (Ly, Lx), True where a cell pixel exists. These pixels are excluded from neuropil masks.	required
inner_neuropil_radius	int, optional (default 2)	Number of pixels to extend around each ROI as an exclusion zone.	2
min_neuropil_pixels	int, optional (default 350)	Minimum number of pixels in each neuropil mask.	350
circular	bool, optional (default False)	If True, extend neuropil masks circularly (slower). If False, extend as rectangles.	False

Returns:

Name	Type	Description
neuropil_ipix	list of numpy.ndarray	Each element is an array of flattened pixel indices for the neuropil mask of one ROI.

Source code in suite2p/extraction/masks.py

def create_neuropil_masks(ypixs, xpixs, cell_pix, inner_neuropil_radius=2,
                          min_neuropil_pixels=350, circular=False):
    """
    Create surround neuropil masks for ROIs by extending each ROI outward.

    Parameters
    ----------
    ypixs : list of numpy.ndarray
        Y-coordinates of the pixels for each ROI.
    xpixs : list of numpy.ndarray
        X-coordinates of the pixels for each ROI.
    cell_pix : numpy.ndarray
        Binary array of shape (Ly, Lx), True where a cell pixel exists.
        These pixels are excluded from neuropil masks.
    inner_neuropil_radius : int, optional (default 2)
        Number of pixels to extend around each ROI as an exclusion zone.
    min_neuropil_pixels : int, optional (default 350)
        Minimum number of pixels in each neuropil mask.
    circular : bool, optional (default False)
        If True, extend neuropil masks circularly (slower). If False, extend
        as rectangles.

    Returns
    -------
    neuropil_ipix : list of numpy.ndarray
        Each element is an array of flattened pixel indices for the neuropil
        mask of one ROI.
    """
    valid_pixels = lambda cell_pix, ypix, xpix: cell_pix[ypix, xpix] < .5
    extend_by = 5

    Ly, Lx = cell_pix.shape
    assert len(xpixs) == len(ypixs)
    neuropil_ipix = []
    idx = 0
    for ypix, xpix in zip(ypixs, xpixs):
        neuropil_mask = np.zeros((Ly, Lx), bool)
        # extend to get ring of dis-allowed pixels
        ypix, xpix = extendROI(ypix, xpix, Ly, Lx, niter=inner_neuropil_radius)
        nring = np.sum(valid_pixels(cell_pix, ypix,
                                    xpix))  # count how many pixels are valid

        nreps = count()
        ypix1, xpix1 = ypix.copy(), xpix.copy()
        while next(nreps) < 100 and np.sum(valid_pixels(
                cell_pix, ypix1, xpix1)) - nring <= min_neuropil_pixels:
            if circular:
                ypix1, xpix1 = extendROI(ypix1, xpix1, Ly, Lx,
                                         extend_by)  # keep extending
            else:
                ypix1, xpix1 = np.meshgrid(
                    np.arange(max(0,
                                  ypix1.min() - extend_by),
                              min(Ly,
                                  ypix1.max() + extend_by + 1), 1, int),
                    np.arange(max(0,
                                  xpix1.min() - extend_by),
                              min(Lx,
                                  xpix1.max() + extend_by + 1), 1, int), indexing="ij")

        ix = valid_pixels(cell_pix, ypix1, xpix1)
        neuropil_mask[ypix1[ix], xpix1[ix]] = True
        neuropil_mask[ypix, xpix] = False
        neuropil_ipix.append(np.ravel_multi_index(np.nonzero(neuropil_mask), (Ly, Lx)))
        idx += 1

    return neuropil_ipix