Cell Detection
Summary
The cell detection algorithm consists of reducing the dimensionality of the data (principal components computation), smoothing spatial principal components, finding peaks in these components, and extending ROIs spatially around these peaks. On each iteration of peak extraction, the neuropil is estimated from large spatial masks and subtracted from the spatial components. This is to improve cell detection and to help avoid extracting neuropil components with large spatial extents.
SVDs ( = PCs) of data
Before computing the principal components of the movie, we bin the data
such that we have at least as many frames to take the SVD of as
specified in the option ops['navg_frames_svd']
. The bin size will be
the maximum of nframes/ops['navg_frames_svd']
and
ops['tau'] * ops['fs']
(the number of samples per transient). We
then bin the movie into this bin size and subtract the mean of the
binned movie across time. Then we smooth the movie in Y and X with a
gaussian filter of standard deviation sig = ops['diameter']/10
. The
we normalize the pixels by their noise variance. The noise variance is
variance of each pixel in the movie across time (at least 1e-10). Then
we compute the covariance of the movie (mov @ mov.T
). Then we
compute the SVD of the covariance and keep the top
ops['nsvd_for_roi']
spatial components (components that are Y x X).
The function that performs this is celldetect2.getSVDdata
and it
requires the ops described above, and Ly, Lx, yrange, xrange, and a
reg_file location.
Sourcery
After the spatial components are found, we perform an iterative algorithm to find the cells in the components. Each iteration consists of the following steps:
Smoothing of spatial components: The components are smoothed with a Gaussian filter in Y and X with standard deviation
sig = ops['diameter']
(this matrix is calledus
). Note that diameter can be a list (for unequal pixel/um in Y and X). Next the mean of the squared smoothed components is computed. The mean of the squared un-smoothed components is also computed. The correlation map is defined as the element-wise division of the smoothed components by the unsmoothed components. The function that computes the correlation map iscelldetect2.getVmap
.Detection of peaks in correlation map: On each iteration, up to 200 peaks are extracted from the correlation map. These are the largest remaining peaks such that they are greater than the threshold, which is set to be proportional to the median of the peaks in the whole correlation map:
ops['threshold_scaling'] * np.median(peaks[peaks>1e-4])
. The initial activitycode
for this newly detected peak is the value ofus
(Gaussian smoothed PCs) at this peak. This is a vector of values across the PCs (nPCs in length).ROI extension: The ROI is iteratively extended around its currently defined pixels +/- 1 in each direction. First, the new pixel weights (
lam
) of the extended ROI are computed. The weightslam
are the unsmoothed PCs projected into thecode
dimension. The pixels that are greater thanmax(lam)/5
are kept. Thelam
’s are normalized to be unit norm. The newcode
is recomputed from the new weights, and is the unsmoothed PCs projected onto thelam
weights. Then this extension procedure is repeated until no pixels are greater thanmax(lam)/5
.Neuropil computation: Now that the new codes are computed, the neuropil is estimated. We set spatial basis functions for the neuropil, which are raised cosines that tile the FOV. The parameter
ops['ratio_neuropil']
determines how big you expect the neuropil basis functions to be relative to the cell diameter (ops['diameter']
). The default is 6. This results in a tiling of 7x7 raised cosines if your FOV is 512x512 pixels and your diameter is 12 pixels. For one-photon recordings, we recommend settingops['ratio_neuropil']
to 2 or 3. Next we perform regression to compute the contribution of the neuropil on the PCs, and we subtract the estimated neuropil contribution from theU
PCs. And these steps are repeated until the stopping criterion is reached.
Stopping criterion: The number of cells detected in the first
iteration is defined as Nfirst
. The cell detection is stopped if the
number of cells detected in the current iteration is less than
Nfirst/10
or if the iteration is the last iteration (defined by
ops['max_iterations']
).
Refinement step: We remove masks which have more than a fraction
ops['max_overlap']
of their pixels overlapping with other masks.
Also, if ops['connected']=1
, then only the connected regions of ROIs
are kept. If you are looking for dendritic components, you may want to
set ops['connected']=0
.