Merge pull request #1379 from flatironinstitute/dev

dev to main for release 1.11.2
flatironinstitute · Jul 25, 2024 · c887cb9 · c887cb9
2 parents e7e8641 + e8abf3b
commit c887cb9
Show file tree

Hide file tree

Showing 60 changed files with 1,004 additions and 692 deletions.
diff --git a/.gitignore b/.gitignore
@@ -52,7 +52,7 @@ __pycache__/
 
 # C extensions
 *.so
-*.c
+caiman/source_extraction/cnmf/oasis.cpp
 
 #movie tif
 *.tif

diff --git a/README.md b/README.md
@@ -13,7 +13,7 @@ Caiman Central
 --------------
 - [Caiman Central](https://github.com/flatironinstitute/caiman_central) is the hub for sharing information about CaImAn. Information on quarterly community meetings, workshops, other events, and any other communications between the developers and the user community can be found there.
 
-# Quick start :rocket:
+# Quick start
 Follow these three steps to get started quickly, from installation to working through a demo notebook. If you do not already have conda installed, [you can find it here](https://docs.conda.io/en/latest/miniconda.html). There is a video walkthrough of the following steps [here](https://youtu.be/b63zAmKihIY?si=m7WleTwdU0rJup_2).
 
 ### Step 1: Install caiman
@@ -41,7 +41,7 @@ Jupyter will open. Navigate to demos/notebooks/ and click on `demo_pipeline.ipyn
 ## For installation help
 Caiman should install easily on Linux, Mac, and Windows. If you run into problems, we have a dedicated [installation page](./docs/source/Installation.rst): the details there should help you troubleshoot. If you don't find what you need there, *please* [create an issue](https://github.com/flatironinstitute/CaImAn/issues) at GitHub, and we will help you get it sorted out. 
 
-# Demo notebooks :page_with_curl:
+# Demo notebooks
 Caiman provides demo notebooks to showcase each of our main features, from motion correction to online CNMF. We recommend starting with the CNMF notebook (`demo_pipeline.ipynb`), which contains more explanation and details than the other notebooks: it covers many concepts that will be used without explanation in the other notebooks. The CNMFE notebook (`demo_pipeline_cnmfE.ipynb`), is also more detailed. Once you've gotten things set up and worked through those "anchor" notebooks, the best way to get started is to work through the demo notebook that most closely matches your use case; you should be able to adapt it for your particular needs.
 
 The main use cases and notebooks are listed in the following table:
@@ -72,7 +72,7 @@ Caiman also provides commandline demos, similar to the notebooks, demonstrating
 - If you have found a bug, we recommend searching the [issues at github](https://github.com/flatironinstitute/CaImAn/issues) and opening a new issue if you can't find the solution there. 
 - If there is a feature you would like to see implemented, feel free to come chat at the above forums or open an issue at Github.
 
-# How to contribute :hammer:
+# How to contribute
  Caiman is an open-source project and improves because of contributions from users all over the world. If there is something about Caiman that you would like to work on, then please reach out. We are always looking for more contributors, so please come read the [contributors page](./CONTRIBUTING.md) for more details about how. 
 
 # Videos 
@@ -89,7 +89,7 @@ The following talks are more in depth:
 * https://www.youtube.com/watch?v=z6TlH28MLRo
 
 
-# Related repositories :pushpin:
+# Related repositories
 There are many repositories that use Caiman, or help make using Caiman easier.
 
 * [use\_cases repo](https://github.com/flatironinstitute/caiman_use_cases):  additional code (unmaintained) demonstrating how to reproduce results in some Caiman-related papers, and how to use/extend Caiman.

diff --git a/bin/caiman_gui.py b/bin/caiman_gui.py
@@ -53,6 +53,8 @@ def make_color_img(img, gain=255, min_max=None, out_type=np.uint8):
 # load object saved by CNMF
 fpath = F.getOpenFileName(caption='Load CNMF Object',
                           filter='HDF5 (*.h5 *.hdf5);;NWB (*.nwb)')[0]
+if fpath == '':
+    raise Exception("You must provide a filename")
 cnm_obj = load_CNMF(fpath)
 
 # movie

diff --git a/caiman/base/movies.py b/caiman/base/movies.py
diff --git a/caiman/base/rois.py b/caiman/base/rois.py
@@ -28,46 +28,36 @@
     pass
 
 
-def com(A: np.ndarray, d1: int, d2: int, d3: Optional[int] = None) -> np.array:
+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
-              matrix of spatial components (d x K)
+         A: np.ndarray or scipy.sparse array or matrix
+              matrix of spatial components (d x K).
 
-         d1:  int
-              number of pixels in x-direction
-
-         d2:  int
-              number of pixels in y-direction
-
-         d3:  int
-              number of pixels in z-direction
+         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 the given dimensions.
 
      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
+    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])
+
+    # take weighted sum of pixel positions along each coordinate
+    cm = (coor @ A / A.sum(axis=0)).T
     return np.array(cm)
 
 
@@ -222,6 +212,7 @@ def nf_match_neurons_in_binary_masks(masks_gt,
             indices false pos
 
     """
+    logger = logging.getLogger("caiman")
 
     _, d1, d2 = np.shape(masks_gt)
     dims = d1, d2
@@ -262,7 +253,7 @@ def nf_match_neurons_in_binary_masks(masks_gt,
     performance['precision'] = TP / (TP + FP)
     performance['accuracy'] = (TP + TN) / (TP + FP + FN + TN)
     performance['f1_score'] = 2 * TP / (2 * TP + FP + FN)
-    logging.debug(performance)
+    logger.debug(performance)
 
     idx_tp = np.where(np.array(costs) < thresh_cost)[0]
     idx_tp_ben = matches[0][idx_tp]    # ground truth
@@ -307,8 +298,8 @@ def nf_match_neurons_in_binary_masks(masks_gt,
             pl.show()
             pl.axis('off')
         except Exception as e:
-            logging.warning("not able to plot precision recall: graphics failure")
-            logging.warning(e)
+            logger.warning("not able to plot precision recall: graphics failure")
+            logger.warning(e)
     return idx_tp_gt, idx_tp_comp, idx_fn_gt, idx_fp_comp, performance
 
 
@@ -402,11 +393,7 @@ def register_ROIs(A1,
             ROIs from session 2 aligned to session 1
 
     """
-
-    #    if 'csc_matrix' not in str(type(A1)):
-    #        A1 = scipy.sparse.csc_matrix(A1)
-    #    if 'csc_matrix' not in str(type(A2)):
-    #        A2 = scipy.sparse.csc_matrix(A2)
+    logger = logging.getLogger("caiman")
 
     if 'ndarray' not in str(type(A1)):
         A1 = A1.toarray()
@@ -498,7 +485,7 @@ def register_ROIs(A1,
     performance['precision'] = TP / (TP + FP)
     performance['accuracy'] = (TP + TN) / (TP + FP + FN + TN)
     performance['f1_score'] = 2 * TP / (2 * TP + FP + FN)
-    logging.info(performance)
+    logger.info(performance)
 
     if plot_results:
         if Cn is None:
@@ -531,11 +518,6 @@ def register_ROIs(A1,
         pl.title('Mismatches')
         pl.axis('off')
 
-
-#        except Exception as e:
-#            logging.warning("not able to plot precision recall usually because we are on travis")
-#            logging.warning(e)
-
     return matched_ROIs1, matched_ROIs2, non_matched1, non_matched2, performance, A2
 
 
@@ -596,6 +578,7 @@ def register_multisession(A,
             by component k in A_union
 
     """
+    logger = logging.getLogger("caiman")
 
     n_sessions = len(A)
     templates = list(templates)
@@ -625,7 +608,7 @@ def register_multisession(A,
                                     enclosed_thr=enclosed_thr)
 
         mat_sess, mat_un, nm_sess, nm_un, _, A2 = reg_results
-        logging.info(len(mat_sess))
+        logger.info(len(mat_sess))
         A_union = A2.copy()
         A_union[:, mat_un] = A[sess][:, mat_sess]
         A_union = np.concatenate((A_union.toarray(), A[sess][:, nm_sess]), axis=1)
@@ -773,6 +756,7 @@ def distance_masks(M_s:list, cm_s: list[list], max_dist: float, enclosed_thr: Op
 
 def find_matches(D_s, print_assignment: bool = False) -> tuple[list, list]:
     # todo todocument
+    logger = logging.getLogger("caiman")
 
     matches = []
     costs = []
@@ -781,7 +765,7 @@ def find_matches(D_s, print_assignment: bool = False) -> tuple[list, list]:
         # we make a copy not to set changes in the original
         DD = D.copy()
         if np.sum(np.where(np.isnan(DD))) > 0:
-            logging.error('Exception: Distance Matrix contains invalid value NaN')
+            logger.error('Exception: Distance Matrix contains invalid value NaN')
             raise Exception('Distance Matrix contains invalid value NaN')
 
         # we do the hungarian
@@ -794,10 +778,10 @@ def find_matches(D_s, print_assignment: bool = False) -> tuple[list, list]:
         for row, column in indexes2:
             value = DD[row, column]
             if print_assignment:
-                logging.debug(('(%d, %d) -> %f' % (row, column, value)))
+                logger.debug(f'({row}, {column}) -> {value}')
             total.append(value)
-        logging.debug(('FOV: %d, shape: %d,%d total cost: %f' % (ii, DD.shape[0], DD.shape[1], np.sum(total))))
-        logging.debug((time.time() - t_start))
+        logger.debug(f'FOV: {ii}, shape: {DD.shape[0]},{DD.shape[1]} total cost: {np.sum(total)}')
+        logger.debug(time.time() - t_start)
         costs.append(total)
         # send back the results in the format we want
     return matches, costs
@@ -828,6 +812,8 @@ def link_neurons(matches: list[list[tuple]],
         neurons: list of arrays representing the indices of neurons in each FOV
 
     """
+    logger = logging.getLogger("caiman")
+
     if min_FOV_present is None:
         min_FOV_present = len(matches)
 
@@ -852,7 +838,7 @@ def link_neurons(matches: list[list[tuple]],
             neurons.append(neuron)
 
     neurons = np.array(neurons).T
-    logging.info(f'num_neurons: {num_neurons}')
+    logger.info(f'num_neurons: {num_neurons}')
     return neurons
 
 
@@ -927,6 +913,8 @@ def nf_read_roi(fileobj) -> np.ndarray:
 
     Adapted from https://gist.github.com/luispedro/3437255
     '''
+    logger = logging.getLogger("caiman")
+
     # This is based on:
     # http://rsbweb.nih.gov/ij/developer/source/ij/io/RoiDecoder.java.html
     # http://rsbweb.nih.gov/ij/developer/source/ij/io/RoiEncoder.java.html
@@ -967,8 +955,7 @@ def getfloat():
 
     magic = fileobj.read(4)
     if magic != 'Iout':
-        #        raise IOError('Magic number not found')
-        logging.warning('Magic number not found')
+        logger.warning('Magic number not found')
     version = get16()
 
     # It seems that the roi type field occupies 2 Bytes, but only one is used
@@ -977,13 +964,6 @@ def getfloat():
     # Discard second Byte:
     get8()
 
-    #    if not (0 <= roi_type < 11):
-    #        logging.error(('roireader: ROI type %s not supported' % roi_type))
-    #
-    #    if roi_type != 7:
-    #
-    #        logging.error(('roireader: ROI type %s not supported (!= 7)' % roi_type))
-
     top = get16()
     left = get16()
     bottom = get16()
@@ -1060,6 +1040,8 @@ def nf_merge_roi_zip(fnames: list[str], idx_to_keep: list[list], new_fold: str):
             name of the output zip file (without .zip extension)
 
     """
+    logger = logging.getLogger("caiman")
+
     folders_rois = []
     files_to_keep = []
     # unzip the files and keep only the ones that are requested
@@ -1068,7 +1050,7 @@ def nf_merge_roi_zip(fnames: list[str], idx_to_keep: list[list], new_fold: str):
         folders_rois.append(dirpath)
         with zipfile.ZipFile(fn) as zf:
             name_rois = zf.namelist()
-            logging.debug(len(name_rois))
+            logger.debug(len(name_rois))
         zip_ref = zipfile.ZipFile(fn, 'r')
         zip_ref.extractall(dirpath)
         files_to_keep.append([os.path.join(dirpath, ff) for ff in np.array(name_rois)[idx]])
@@ -1119,6 +1101,8 @@ def extract_binary_masks_blob(A,
         neg_examples:
 
     """
+    logger = logging.getLogger("caiman")
+
     params = cv2.SimpleBlobDetector_Params()
     params.minCircularity = minCircularity
     params.minInertiaRatio = minInertiaRatio
@@ -1146,7 +1130,7 @@ def extract_binary_masks_blob(A,
     neg_examples = []
 
     for count, comp in enumerate(A.tocsc()[:].T):
-        logging.debug(count)
+        logger.debug(count)
         comp_d = np.array(comp.todense())
         gray_image = np.reshape(comp_d, dims, order='F')
         gray_image = (gray_image - np.min(gray_image)) / \
@@ -1171,7 +1155,7 @@ def extract_binary_masks_blob(A,
             edges = (label_objects == (1 + idx_largest))
             edges = scipy.ndimage.binary_fill_holes(edges)
         else:
-            logging.warning('empty component')
+            logger.warning('empty component')
             edges = np.zeros_like(edges)
 
         masks_ws.append(edges)
@@ -1275,14 +1259,15 @@ def detect_duplicates_and_subsets(binary_masks,
                                   dist_thr: float = 0.1,
                                   min_dist=10,
                                   thresh_subset: float = 0.8):
+    logger = logging.getLogger("caiman")
 
     cm = [scipy.ndimage.center_of_mass(mm) for mm in binary_masks]
     sp_rois = scipy.sparse.csc_matrix(np.reshape(binary_masks, (binary_masks.shape[0], -1)).T)
     D = distance_masks([sp_rois, sp_rois], [cm, cm], min_dist)[0]
     np.fill_diagonal(D, 1)
     overlap = sp_rois.T.dot(sp_rois).toarray()
     sz = np.array(sp_rois.sum(0))
-    logging.info(sz.shape)
+    logger.info(sz.shape)
     overlap = overlap / sz.T
     np.fill_diagonal(overlap, 0)
     # pairs of duplicate indices
@@ -1297,7 +1282,7 @@ def detect_duplicates_and_subsets(binary_masks,
         metric = r_values.squeeze()
     else:
         metric = sz.squeeze()
-        logging.debug('***** USING MAX AREA BY DEFAULT')
+        logger.debug('***** USING MAX AREA BY DEFAULT')
 
     overlap_tmp = overlap.copy() >= thresh_subset
     overlap_tmp = overlap_tmp * metric[:, None]