From 522d456be8eeca6ac12f73b8de872521a92695d1 Mon Sep 17 00:00:00 2001 From: Ethan Blackwood Date: Fri, 28 Jun 2024 11:50:53 -0400 Subject: [PATCH 1/7] Fix CoM calculation when swap_dim is true --- caiman/utils/visualization.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/caiman/utils/visualization.py b/caiman/utils/visualization.py index a82833507..825757a1b 100644 --- a/caiman/utils/visualization.py +++ b/caiman/utils/visualization.py @@ -392,18 +392,16 @@ def get_contours(A, dims, thr=0.9, thr_method='nrg', swap_dim=False): if 'csc_matrix' not in str(type(A)): A = csc_matrix(A) d, nr = np.shape(A) - # if we are on a 3D video - if len(dims) == 3: - d1, d2, d3 = dims - x, y = np.mgrid[0:d2:1, 0:d3:1] - else: - d1, d2 = dims - x, y = np.mgrid[0:d1:1, 0:d2:1] + d1, d2 = dims[:2] coordinates = [] # get the center of mass of neurons( patches ) - cm = caiman.base.rois.com(A, *dims) + # com assumes F-order + if swap_dim: + cm = caiman.base.rois.com(A, *dims[::-1])[:, ::-1] + else: + cm = caiman.base.rois.com(A, *dims) # for each patches for i in range(nr): From 1f88e493a6756fd0bb6dca580c930c32110a192d Mon Sep 17 00:00:00 2001 From: Ethan Blackwood Date: Tue, 2 Jul 2024 12:13:26 -0400 Subject: [PATCH 2/7] Fix plot_contours to use CoMs from updated get_contours --- caiman/utils/visualization.py | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/caiman/utils/visualization.py b/caiman/utils/visualization.py index 825757a1b..7f03e4652 100644 --- a/caiman/utils/visualization.py +++ b/caiman/utils/visualization.py @@ -1096,16 +1096,11 @@ def plot_contours(A, Cn, thr=None, thr_method='max', maxthr=0.2, nrgthr=0.9, dis plt.plot(*v.T, c=colors, **contour_args) if display_numbers: - d1, d2 = np.shape(Cn) - d, nr = np.shape(A) - cm = caiman.base.rois.com(A, d1, d2) + nr = A.shape[1] if max_number is None: - max_number = A.shape[1] - for i in range(np.minimum(nr, max_number)): - if swap_dim: - ax.text(cm[i, 0], cm[i, 1], str(i + 1), color=colors, **number_args) - else: - ax.text(cm[i, 1], cm[i, 0], str(i + 1), color=colors, **number_args) + max_number = nr + for i, c in zip(range(np.minimum(nr, max_number)), coordinates): + ax.text(c['CoM'][1], c['CoM'][0], str(i + 1), color=colors, **number_args) return coordinates def plot_shapes(Ab, dims, num_comps=15, size=(15, 15), comps_per_row=None, From b38f0faa75859536e7c889436a86c7437d3577b1 Mon Sep 17 00:00:00 2001 From: Ethan Blackwood Date: Tue, 2 Jul 2024 16:22:12 -0400 Subject: [PATCH 3/7] Add rewritten com function with C/F order and add test --- caiman/base/rois.py | 44 ++++++++++++----------------------- caiman/tests/test_toydata.py | 17 ++++++++++++++ caiman/utils/visualization.py | 6 +---- 3 files changed, 33 insertions(+), 34 deletions(-) diff --git a/caiman/base/rois.py b/caiman/base/rois.py index b5fff51e1..86e0fd6b0 100644 --- a/caiman/base/rois.py +++ b/caiman/base/rois.py @@ -28,46 +28,32 @@ pass -def com(A: np.ndarray, d1: int, d2: int, d3: Optional[int] = None) -> np.array: +def com(A, *dims: int, order='F') -> np.ndarray: """Calculation of the center of mass for spatial components Args: - A: np.ndarray - matrix of spatial components (d x K) + A: np.ndarray or scipy.sparse array or matrix + matrix of spatial components (d x K); column indices map to dimensions according to F-order. - d1: int - number of pixels in x-direction - - d2: int - number of pixels in y-direction - - d3: int - number of pixels in z-direction + *dims: D ints + each positional argument after A defines the extent of one of the dimensions of the data. + + order: 'C' or 'F' + how each column of A should be reshaped to match dims. Returns: cm: np.ndarray - center of mass for spatial components (K x 2 or 3) + center of mass for spatial components (K x D) """ - if 'csc_matrix' not in str(type(A)): A = scipy.sparse.csc_matrix(A) - if d3 is None: - Coor = np.matrix([np.outer(np.ones(d2), np.arange(d1)).ravel(), - np.outer(np.arange(d2), np.ones(d1)).ravel()], - dtype=A.dtype) - else: - Coor = np.matrix([ - np.outer(np.ones(d3), - np.outer(np.ones(d2), np.arange(d1)).ravel()).ravel(), - np.outer(np.ones(d3), - np.outer(np.arange(d2), np.ones(d1)).ravel()).ravel(), - np.outer(np.arange(d3), - np.outer(np.ones(d2), np.ones(d1)).ravel()).ravel() - ], - dtype=A.dtype) - - cm = (Coor * A / A.sum(axis=0)).T + # make coordinate arrays where coor[d] increases from 0 to npixels[d]-1 along the dth axis + coors = np.meshgrid(*[range(d) for d in dims], indexing='ij') + coor = np.stack([c.ravel(order=order) for c in coors]) + + # take weighted sum of pixel positions along each coordinate + cm = (coor @ A / A.sum(axis=0)).T return np.array(cm) diff --git a/caiman/tests/test_toydata.py b/caiman/tests/test_toydata.py index 2a0f2a8e0..bf61fd7d3 100644 --- a/caiman/tests/test_toydata.py +++ b/caiman/tests/test_toydata.py @@ -6,6 +6,7 @@ import caiman.source_extraction.cnmf.params from caiman.source_extraction import cnmf as cnmf +from caiman.utils.visualization import get_contours #%% @@ -64,6 +65,22 @@ def pipeline(D): ] npt.assert_allclose(corr, 1, .05) + # Check that get_contours works regardless of swap_dim + if D == 2: # Can't expect swap_dim to work for 3D data in the same way + coor_normal = get_contours(cnm.estimates.A, dims, swap_dim=False) + coor_swapped = get_contours(cnm.estimates.A, dims[::-1], swap_dim=True) + for c_normal, c_swapped in zip(coor_normal, coor_swapped): + # have to sort coordinates b/c starting point is unimportant & depend on orientation + # also, the first point is repeated and this may be a different point depending on orientation. + # to get around this, compare differences instead (have to take absolute value b/c direction may be opposite) + def normalize_coords(coords): + abs_diffs = np.abs(np.diff(coords, axis=0)) + sort_order = np.lexsort(abs_diffs.T) + return abs_diffs[sort_order, :] + + npt.assert_allclose(normalize_coords(c_normal['coordinates']), normalize_coords(c_swapped['coordinates'][:, ::-1])) + npt.assert_allclose(c_normal['CoM'], c_swapped['CoM'][::-1]) + def test_2D(): pipeline(2) diff --git a/caiman/utils/visualization.py b/caiman/utils/visualization.py index 7f03e4652..0b6bbb652 100644 --- a/caiman/utils/visualization.py +++ b/caiman/utils/visualization.py @@ -397,11 +397,7 @@ def get_contours(A, dims, thr=0.9, thr_method='nrg', swap_dim=False): coordinates = [] # get the center of mass of neurons( patches ) - # com assumes F-order - if swap_dim: - cm = caiman.base.rois.com(A, *dims[::-1])[:, ::-1] - else: - cm = caiman.base.rois.com(A, *dims) + cm = caiman.base.rois.com(A, *dims, order='C' if swap_dim else 'F') # for each patches for i in range(nr): From 467c591af6ff749318afbef19b7163704848eaf4 Mon Sep 17 00:00:00 2001 From: Ethan Blackwood Date: Fri, 5 Jul 2024 10:10:43 -0400 Subject: [PATCH 4/7] Restore individual dimension arguments for com --- caiman/base/rois.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/caiman/base/rois.py b/caiman/base/rois.py index 86e0fd6b0..417405828 100644 --- a/caiman/base/rois.py +++ b/caiman/base/rois.py @@ -28,18 +28,18 @@ pass -def com(A, *dims: int, order='F') -> np.ndarray: +def com(A, d1: int, d2: int, d3: Optional[int] = None, order: str = 'F') -> np.ndarray: """Calculation of the center of mass for spatial components Args: A: np.ndarray or scipy.sparse array or matrix - matrix of spatial components (d x K); column indices map to dimensions according to F-order. + matrix of spatial components (d x K). - *dims: D ints - each positional argument after A defines the extent of one of the dimensions of the data. + d1, d2, d3: ints + d1, d2, and (optionally) d3 are the original dimensions of the data. order: 'C' or 'F' - how each column of A should be reshaped to match dims. + how each column of A should be reshaped to match the given dimensions. Returns: cm: np.ndarray @@ -48,6 +48,10 @@ def com(A, *dims: int, order='F') -> np.ndarray: if 'csc_matrix' not in str(type(A)): A = scipy.sparse.csc_matrix(A) + dims = [d1, d2] + if d3 is not None: + dims.append(d3) + # make coordinate arrays where coor[d] increases from 0 to npixels[d]-1 along the dth axis coors = np.meshgrid(*[range(d) for d in dims], indexing='ij') coor = np.stack([c.ravel(order=order) for c in coors]) From d8835987801667615cdb9ef31782032c6263a7e8 Mon Sep 17 00:00:00 2001 From: Ethan Blackwood Date: Fri, 5 Jul 2024 18:58:45 -0400 Subject: [PATCH 5/7] Add slice_dim option for get_contours and 3D test for swap_dim --- caiman/tests/test_toydata.py | 38 +++++++++++++++++++++++------------ caiman/utils/visualization.py | 37 +++++++++++++++++++++++++--------- 2 files changed, 53 insertions(+), 22 deletions(-) diff --git a/caiman/tests/test_toydata.py b/caiman/tests/test_toydata.py index bf61fd7d3..beb51a786 100644 --- a/caiman/tests/test_toydata.py +++ b/caiman/tests/test_toydata.py @@ -66,20 +66,32 @@ def pipeline(D): npt.assert_allclose(corr, 1, .05) # Check that get_contours works regardless of swap_dim - if D == 2: # Can't expect swap_dim to work for 3D data in the same way - coor_normal = get_contours(cnm.estimates.A, dims, swap_dim=False) - coor_swapped = get_contours(cnm.estimates.A, dims[::-1], swap_dim=True) - for c_normal, c_swapped in zip(coor_normal, coor_swapped): - # have to sort coordinates b/c starting point is unimportant & depend on orientation - # also, the first point is repeated and this may be a different point depending on orientation. - # to get around this, compare differences instead (have to take absolute value b/c direction may be opposite) - def normalize_coords(coords): - abs_diffs = np.abs(np.diff(coords, axis=0)) - sort_order = np.lexsort(abs_diffs.T) - return abs_diffs[sort_order, :] + coor_normal = get_contours(cnm.estimates.A, dims, swap_dim=False) + coor_swapped = get_contours(cnm.estimates.A, dims[::-1], swap_dim=True) + for c_normal, c_swapped in zip(coor_normal, coor_swapped): + if D == 3: + for plane_coor_normal, plane_coor_swapped in zip(c_normal['coordinates'], c_swapped['coordinates']): + compare_contour_coords(plane_coor_normal, plane_coor_swapped[:, ::-1]) + else: + compare_contour_coords(c_normal['coordinates'], c_swapped['coordinates'][:, ::-1]) - npt.assert_allclose(normalize_coords(c_normal['coordinates']), normalize_coords(c_swapped['coordinates'][:, ::-1])) - npt.assert_allclose(c_normal['CoM'], c_swapped['CoM'][::-1]) + npt.assert_allclose(c_normal['CoM'], c_swapped['CoM'][::-1]) + +def compare_contour_coords(coords1: np.ndarray, coords2: np.ndarray): + """ + Compare 2 matrices of contour coordinates that should be the same, but may be calculated in a different order/ + from different starting points. + + The first point of each contour component is repeated, and this may be a different point depending on orientation. + To get around this, compare differences instead (have to take absolute value b/c direction may be opposite). + Also sort coordinates b/c starting point is unimportant & depends on orientation + """ + diffs_sorted = [] + for coords in [coords1, coords2]: + abs_diffs = np.abs(np.diff(coords, axis=0)) + sort_order = np.lexsort(abs_diffs.T) + diffs_sorted.append(abs_diffs[sort_order, :]) + npt.assert_allclose(diffs_sorted[0], diffs_sorted[1]) def test_2D(): diff --git a/caiman/utils/visualization.py b/caiman/utils/visualization.py index 0b6bbb652..cceb58f59 100644 --- a/caiman/utils/visualization.py +++ b/caiman/utils/visualization.py @@ -20,7 +20,7 @@ from skimage.measure import find_contours import sys from tempfile import NamedTemporaryFile -from typing import Any, Optional +from typing import Any, Optional, Literal, Union from warnings import warn import caiman.base.rois @@ -366,7 +366,7 @@ def plot_unit(uid, scl): .redim.range(unit_id=(0, nr-1), scale=(0.0, 1.0))) -def get_contours(A, dims, thr=0.9, thr_method='nrg', swap_dim=False): +def get_contours(A, dims, thr=0.9, thr_method='nrg', swap_dim=False, slice_dim: Union[int, Literal['auto']] = 'auto'): """Gets contour of spatial components and returns their coordinates Args: @@ -374,15 +374,24 @@ def get_contours(A, dims, thr=0.9, thr_method='nrg', swap_dim=False): Matrix of Spatial components (d x K) dims: tuple of ints - Spatial dimensions of movie (x, y[, z]) + Spatial dimensions of movie thr: scalar between 0 and 1 Energy threshold for computing contours (default 0.9) - thr_method: [optional] string + thr_method: string Method of thresholding: 'max' sets to zero pixels that have value less than a fraction of the max value 'nrg' keeps the pixels that contribute up to a specified fraction of the energy + + swap_dim: bool + If False (default), each column of A should be reshaped in F-order to recover the mask; + this is correct if the dimensions have not been reordered from (y, x[, z]). + If True, each column should be reshaped in C-order; this is correct for dims = ([z, ]x, y). + + slice_dim: int or 'auto' + Which dimension to slice along if we have 3D data. (i.e., get contours on each plane along this axis). + The default ('auto') is 0 if swap_dim is True, else -1. Returns: Coor: list of coordinates with center of mass and @@ -431,9 +440,9 @@ def get_contours(A, dims, thr=0.9, thr_method='nrg', swap_dim=False): Bmat = np.reshape(Bvec, dims, order='C') else: Bmat = np.reshape(Bvec, dims, order='F') - pars['coordinates'] = [] - # for each dimensions we draw the contour - for B in (Bmat if len(dims) == 3 else [Bmat]): + + def get_slice_coords(B: np.ndarray) -> np.ndarray: + """Get contour coordinates for a 2D slice""" vertices = find_contours(B.T, thr) # this fix is necessary for having disjoint figures and borders plotted correctly v = np.atleast_2d([np.nan, np.nan]) @@ -449,9 +458,19 @@ def get_contours(A, dims, thr=0.9, thr_method='nrg', swap_dim=False): vtx = np.concatenate((vtx, vtx[0, np.newaxis]), axis=0) v = np.concatenate( (v, vtx, np.atleast_2d([np.nan, np.nan])), axis=0) + return v + + if len(dims) == 2: + pars['coordinates'] = get_slice_coords(Bmat) + else: + # make a list of the contour coordinates for each 2D slice + pars['coordinates'] = [] + if slice_dim == 'auto': + slice_dim = 0 if swap_dim else -1 + for s in range(dims[slice_dim]): + B = Bmat.take(s, axis=slice_dim) + pars['coordinates'].append(get_slice_coords(B)) - pars['coordinates'] = v if len( - dims) == 2 else (pars['coordinates'] + [v]) pars['CoM'] = np.squeeze(cm[i, :]) pars['neuron_id'] = i + 1 coordinates.append(pars) From 830967c747334fc3bd2c620d96a3714590415ac0 Mon Sep 17 00:00:00 2001 From: Ethan Blackwood Date: Fri, 5 Jul 2024 19:48:02 -0400 Subject: [PATCH 6/7] Minor changes to make contours more consistent/predictable --- caiman/utils/visualization.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/caiman/utils/visualization.py b/caiman/utils/visualization.py index cceb58f59..707781d0e 100644 --- a/caiman/utils/visualization.py +++ b/caiman/utils/visualization.py @@ -401,7 +401,6 @@ def get_contours(A, dims, thr=0.9, thr_method='nrg', swap_dim=False, slice_dim: if 'csc_matrix' not in str(type(A)): A = csc_matrix(A) d, nr = np.shape(A) - d1, d2 = dims[:2] coordinates = [] @@ -443,6 +442,7 @@ def get_contours(A, dims, thr=0.9, thr_method='nrg', swap_dim=False, slice_dim: def get_slice_coords(B: np.ndarray) -> np.ndarray: """Get contour coordinates for a 2D slice""" + d1, d2 = B.shape vertices = find_contours(B.T, thr) # this fix is necessary for having disjoint figures and borders plotted correctly v = np.atleast_2d([np.nan, np.nan]) @@ -451,8 +451,8 @@ def get_slice_coords(B: np.ndarray) -> np.ndarray: if num_close_coords < 2: if num_close_coords == 0: # case angle - newpt = np.round(vtx[-1, :] / [d2, d1]) * [d2, d1] - vtx = np.concatenate((vtx, newpt[np.newaxis, :]), axis=0) + newpt = np.round(np.mean(vtx[[0, -1], :], axis=0) / [d2, d1]) * [d2, d1] + vtx = np.concatenate((newpt[np.newaxis, :], vtx, newpt[np.newaxis, :]), axis=0) else: # case one is border vtx = np.concatenate((vtx, vtx[0, np.newaxis]), axis=0) From 04d1ca20da4a41322ac4d4530654cff28b733632 Mon Sep 17 00:00:00 2001 From: Ethan Blackwood Date: Mon, 8 Jul 2024 10:20:24 -0400 Subject: [PATCH 7/7] Use None instead of 'auto' --- caiman/utils/visualization.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/caiman/utils/visualization.py b/caiman/utils/visualization.py index 707781d0e..9f9a6f2e1 100644 --- a/caiman/utils/visualization.py +++ b/caiman/utils/visualization.py @@ -20,7 +20,7 @@ from skimage.measure import find_contours import sys from tempfile import NamedTemporaryFile -from typing import Any, Optional, Literal, Union +from typing import Any, Optional from warnings import warn import caiman.base.rois @@ -366,7 +366,7 @@ def plot_unit(uid, scl): .redim.range(unit_id=(0, nr-1), scale=(0.0, 1.0))) -def get_contours(A, dims, thr=0.9, thr_method='nrg', swap_dim=False, slice_dim: Union[int, Literal['auto']] = 'auto'): +def get_contours(A, dims, thr=0.9, thr_method='nrg', swap_dim=False, slice_dim: Optional[int] = None): """Gets contour of spatial components and returns their coordinates Args: @@ -389,9 +389,9 @@ def get_contours(A, dims, thr=0.9, thr_method='nrg', swap_dim=False, slice_dim: this is correct if the dimensions have not been reordered from (y, x[, z]). If True, each column should be reshaped in C-order; this is correct for dims = ([z, ]x, y). - slice_dim: int or 'auto' + slice_dim: int or None Which dimension to slice along if we have 3D data. (i.e., get contours on each plane along this axis). - The default ('auto') is 0 if swap_dim is True, else -1. + The default (None) is 0 if swap_dim is True, else -1. Returns: Coor: list of coordinates with center of mass and @@ -465,7 +465,7 @@ def get_slice_coords(B: np.ndarray) -> np.ndarray: else: # make a list of the contour coordinates for each 2D slice pars['coordinates'] = [] - if slice_dim == 'auto': + if slice_dim is None: slice_dim = 0 if swap_dim else -1 for s in range(dims[slice_dim]): B = Bmat.take(s, axis=slice_dim)