Cellatlas (#5)

* refactor get cache function to return higher level instead of `flatmaps` folder

* wip plenty of stuff

* make sure slice dimension are consistent regardless of the volume memory layout

* add the agea atlas backend with geometry

* bugfixes

* windows CI tests require memmap to be closed to remove tempdir

* review comments modifications

.flake8

@@ -1,6 +1,6 @@
 [flake8]
 max-line-length = 130
-ignore = W504, W503, E266
+ignore = W504, W503, E266, D, BLK
 exclude =
     .git,
     __pycache__,

CHANGELOG.md

@@ -1,3 +1,13 @@
+## [1.0.0]
+
+### Added
+- `iblatlas.genomics` module for working with genomics data from Allen contains
+  - the Allen gene expression atlas in the `iblatlas.genomics.agea` module
+  - the Allen cell types atlas in the `iblatlas.genomics.merfish` module
+### Modified
+- slices of the atlas are now always returned with consistent sizes regardless of the volume layout on disk
+- atlases now can have an extra dimension in the image volume, to allow for multiple layers
+
 ## [0.3.0]
 ### Modified
 - Insertion._get_surface_intersection: can return None if and when mode is not set to raise and insertion is not intersecting with any surface

README.md

@@ -1,20 +1,24 @@
 # iblatlas
-This repository contains tools to manipulate 3D representations of the mouse brain anatomy for electrophysiolgy experiments.
-The tools are mainly using the the Allen Brain Atlas although some other atlases can be used and are written in Python.
+Tools to manipulate hierarchical 3D representations of the mouse brain anatomy for electrophysiolgy experiments.
+The tools are mainly using the the Allen CCF although some other atlases can be used.
 
-One of the requirment of this repository is to use minimal requirements, based of standard matplotlib, numpy and scipy libraries, and excludes visualisation tool.
+**This repository uses minimal requirements, based of standard `matplotlib`, `numpy` and `scipy` libraries, to the exclusion
+of more complex visualization tools such as `pyqt`.**
 
-The documentation can be found here: https://int-brain-lab.github.io/iblenv/notebooks_external/atlas_working_with_ibllib_atlas.html
+## Documentation
+[Full package documentation](https://int-brain-lab.github.io/iblenv/_autosummary/iblatlas.html#module-iblatlas)
+
+[Notebooks and examples](https://int-brain-lab.github.io/iblenv/notebooks_external/atlas_working_with_ibllib_atlas.html)
 
 ## Installation
 `pip install iblatlas`
 
-
 ## Contributing
 Changes are merged by pull requests.
 Release checklist:
-- [ ] Update version in `iblatlas/__init__.py`
-- [ ] Update `CHANGELOG.md`
-- [ ] Create a pull request to the `main` branch on GitHub
+- [x] Update version in `iblatlas/__init__.py`
+- [x] Update `CHANGELOG.md`
+- [x] Create a pull request to the `main` branch on GitHub
+- [x] Once the PR is merged, create a new tag and push the tag
 
-Once merged the package is uploaded to PyPI using GitHub Actions.
+Once a tag is pushed on main the package is uploaded to PyPI using GitHub Actions.

examples/atlas_genomics_load_allen_gene_expression_atlas.py

@@ -0,0 +1,14 @@
+import matplotlib.pyplot as plt
+from iblatlas.genomics import agea
+
+df_genes, gene_expression_volumes, atlas_agea = agea.load()
+
+igenes = (0,)
+fig, axs = plt.subplots(3, 2, sharex=True, sharey=True)
+atlas_agea.plot_cslice(0, ax=axs[0, 0])
+atlas_agea.plot_cslice(0, ax=axs[1, 0], volume='annotation')
+atlas_agea.plot_cslice(0, ax=axs[2, 0], volume=gene_expression_volumes[igenes[0]], cmap='viridis')
+atlas_agea.plot_sslice(0, ax=axs[0, 1])
+atlas_agea.plot_sslice(0, ax=axs[1, 1], volume='annotation')
+atlas_agea.plot_sslice(0, ax=axs[2, 1], volume=gene_expression_volumes[igenes[0]], cmap='viridis')
+fig.tight_layout()

iblatlas/__init__.py

@@ -3,8 +3,7 @@
 For examples and tutorials on using the IBL atlas package, see
 https://docs.internationalbrainlab.org/atlas_examples.html
 
-.. TODO Explain differences between this package and the Allen SDK.
-Much of this was adapted from the `cortexlab allenCCF repository <https://github.com/cortex-lab/allenCCF>`_.
+
 
 Terminology
 -----------
@@ -61,11 +60,15 @@
 bregma and lambda are at the same DV height) [6]_, however this is not currently accounted for.
 
 
-Mappings
---------
-In addition to the atlases there are also multiple brain region mappings that serve one of two purposes: 1. control the
-granularity particular brain regions; 2. support differing anatomical sub-devisions or nomenclature.  The two Allen atlas mappings
-below were created somewhat arbirarily by Nick Steinmetz to aid in analysis:
+Brain Regions Mappings
+----------------------
+The Allen CCF ontology regroups brain regions in a hierarchical manner, with up to 9 levels of parcellation. The
+labels volumes provided by the Allen Institute are at the finest level of parcellation, ie. leaf nodes.
+However it is often useful to regroup regions at higher levels of the hierarchy for aggregation purposes. The mappings
+represent the correspondence between the leaf nodes and their grouping for given analysis purposes. This is useful for:
+ 1. control the granularity particular brain regions;
+ 2. support differing anatomical sub-devisions or nomenclature.
+The two Allen atlas mappings below were created by Nick Steinmetz to aid in neural analysis:
 
 1. Beryl - brain atlas annotations without layer sub-divisions or certain ganglial/nucleus sub-devisisions (e.g. the core/shell
    sub-division of the lateral geniculate nucleus). Fibre tracts, pia, etc. are also absent.  The choice of which areas to combine
@@ -90,7 +93,7 @@
 Notes
 -----
 The Allen atlas and the CCF annotations have different release dates and versions [10]_. The annotations used by IBL are the 2017
-version.
+version (CCFv2).
 
 The IBL uses the following conventions:
 
@@ -107,7 +110,7 @@
 
 Examples
 --------
-Below are some breif API examples. For in depth tutorials on using the IBL atlas package, see
+Below are some brief API examples. For in depth tutorials on using the IBL atlas package, see
 https://docs.internationalbrainlab.org/atlas_examples.html.
 
 Find bregma position in indices * resolution in um
@@ -138,7 +141,7 @@
   `Mappings`_.
 * **swanson_regions.npy** - A 1D array of length 323 containing the Allen CCF brain region IDs
 * **mappings.pqt** - A table of mappings.  Each column defines a mapping, with the '-lr' suffix indicating a lateralized version.
-  The rows contain the correspondence of each mapping to the int64 index of the lateralized Allen structure tree.  The table is
+  The rows contain the correspondence of each mapping to the int64 index of the lateralized Allen structure tree. The table is
   generated by iblatlas.regions.BrainRegions._compute_mappings.
 
 Remote files
@@ -159,8 +162,6 @@
   region were then simplified using the `Ramer Douglas Peucker algorithm <https://rdp.readthedocs.io/en/latest/>`_
 * **swanson2allen.npz** - TODO Document who made this, its contents, purpose and data type
 * **<flatmap_name>_<res_um>.nrrd** - TODO Document who made this, its contents, purpose and data type
-* **gene-expression.pqt** - TODO Document who made this, its contents, purpose and data type
-* **gene-expression.bin** - TODO Document who made this, its contents, purpose and data type.
 
 .. [*] The annotation and average template volumes were created from the images provided in Supplemtary Data 4 of Chon et al. [3]_
    and stitched together as a single volume using SimpleITK.
@@ -192,6 +193,5 @@
    ontology and flatmaps. J Comp Neurol. [doi 10.1002/cne.24381]
 .. [10] Allen Mouse Common Coordinate Framework Technical White Paper (October 2017 v3)
    http://help.brain-map.org/download/attachments/8323525/Mouse_Common_Coordinate_Framework.pdf
-
 """
-__version__ = '0.3.0'
+__version__ = '0.4.0'

iblatlas/atlas.py

@@ -148,19 +148,6 @@ def nxyz(self):
         """numpy.array: Coordinates of the element volume[0, 0, 0] in the coordinate space."""
         return np.array([self.nx, self.ny, self.nz])
 
-    """Methods ratios to indices"""
-    def r2ix(self, r):
-        # FIXME Document
-        return int((self.nx - 1) * r)
-
-    def r2iy(self, r):
-        # FIXME Document
-        return int((self.nz - 1) * r)
-
-    def r2iz(self, r):
-        # FIXME Document
-        return int((self.nz - 1) * r)
-
     """Methods distance to indices"""
     @staticmethod
     def _round(i, round=True):
@@ -459,10 +446,13 @@ class BrainAtlas:
     yet this class can be used for other atlases arises.
     """
 
-    """numpy.array: An image volume."""
+    """numpy.array: A 3D image volume."""
     image = None
-    """numpy.array: An annotation label volume."""
+    """numpy.array: A 3D annotation label volume."""
     label = None
+    """numpy.array: One or several optional data volumes, the 3 last dimensions should match
+    the image and the label volumes dimensions"""
+    volumes = None
 
     def __init__(self, image, label, dxyz, regions, iorigin=[0, 0, 0],
                  dims2xyz=[0, 1, 2], xyz2dims=[0, 1, 2]):
@@ -478,8 +468,8 @@ def __init__(self, image, label, dxyz, regions, iorigin=[0, 0, 0],
         self.image = image
         self.label = label
         self.regions = regions
-        self.dims2xyz = dims2xyz
-        self.xyz2dims = xyz2dims
+        self.dims2xyz = np.array(dims2xyz)
+        self.xyz2dims = np.array(xyz2dims)
         assert np.all(self.dims2xyz[self.xyz2dims] == np.array([0, 1, 2]))
         assert np.all(self.xyz2dims[self.dims2xyz] == np.array([0, 1, 2]))
         # create the coordinate transform object that maps volume indices to real world coordinates
@@ -492,8 +482,13 @@ def __init__(self, image, label, dxyz, regions, iorigin=[0, 0, 0],
 
     @staticmethod
     def _get_cache_dir():
+        """
+        ./histology/ATLAS/Needles/Allen
+        Where . is the main ONE cache directory
+        :return: pathlib.Path
+        """
         par = one.params.get(silent=True)
-        path_atlas = Path(par.CACHE_DIR).joinpath('histology', 'ATLAS', 'Needles', 'Allen', 'flatmaps')
+        path_atlas = Path(par.CACHE_DIR).joinpath('histology', 'ATLAS', 'Needles', 'Allen')
         return path_atlas
 
     def mask(self):
@@ -729,7 +724,7 @@ def _plot_slice(im, extent, ax=None, cmap=None, volume=None, **kwargs):
         cmap : str, matplotlib.colors.Colormap
             The Colormap instance or registered colormap name used to map scalar data to colors.
             Defaults to 'bone'.
-        volume : str
+        volume : str | np.array
             If 'boundary', assumes image is an outline of boundaries between all regions.
             FIXME How does this affect the plot?
         **kwargs
@@ -746,10 +741,11 @@ def _plot_slice(im, extent, ax=None, cmap=None, volume=None, **kwargs):
         if not cmap:
             cmap = plt.get_cmap('bone')
 
-        if volume == 'boundary':
-            imb = np.zeros((*im.shape[:2], 4), dtype=np.uint8)
-            imb[im == 1] = np.array([0, 0, 0, 255])
-            im = imb
+        if isinstance(volume, str):
+            if volume == 'boundary':
+                imb = np.zeros((*im.shape[:2], 4), dtype=np.uint8)
+                imb[im == 1] = np.array([0, 0, 0, 255])
+                im = imb
 
         ax.imshow(im, extent=extent, cmap=cmap, **kwargs)
         return ax
@@ -795,7 +791,10 @@ def slice(self, coordinate, axis, volume='image', mode='raise', region_values=No
             - if volume='volume', region_values must have shape ba.image.shape
             - if volume='value', region_values must have shape ba.regions.id
         :param mapping: mapping to use. Options can be found using ba.regions.mappings.keys()
-        :return: 2d array or 3d RGB numpy int8 array
+        :return: 2d array or 3d RGB numpy int8 array of dimensions:
+            - 0: nap x ndv (sagittal slice)
+            - 1: nml x ndv (coronal slice)
+            - 2: nap x nml (horizontal slice)
         """
         if axis == 0:
             index = self.bc.x2i(np.array(coordinate), mode=mode)
@@ -804,43 +803,53 @@ def slice(self, coordinate, axis, volume='image', mode='raise', region_values=No
         elif axis == 2:
             index = self.bc.z2i(np.array(coordinate), mode=mode)
 
-        # np.take is 50 thousand times slower than straight slicing !
         def _take(vol, ind, axis):
+            """
+            This is a 2 steps process to get the slice along the correct axis
+            1) slice the volume according to the mapped axis corresponding to the sclice
+            we do this because np.take is 50 thousand times slower than straight slicing !
+            2) reshape the output array according to the slice specifications
+            """
+            volume_axis = self.xyz2dims[axis]
             if mode == 'clip':
-                ind = np.minimum(np.maximum(ind, 0), vol.shape[axis] - 1)
-            if axis == 0:
-                return vol[ind, :, :]
-            elif axis == 1:
-                return vol[:, ind, :]
-            elif axis == 2:
-                return vol[:, :, ind]
+                ind = np.minimum(np.maximum(ind, 0), vol.shape[volume_axis] - 1)
+            if volume_axis == 0:
+                slic = vol[ind, :, :]
+            elif volume_axis == 1:
+                slic = vol[:, ind, :]
+            elif volume_axis == 2:
+                slic = vol[:, :, ind]
+            output_sizes = [[1, 2], [0, 2], [1, 0]]  # we expect those sizes where index is the axis
+            if np.diff(self.xyz2dims[output_sizes[axis]])[0] < 0:
+                slic = slic.transpose().copy()
+            return slic
 
         def _take_remap(vol, ind, axis, mapping):
             # For the labels, remap the regions indices according to the mapping
             return self._get_mapping(mapping=mapping)[_take(vol, ind, axis)]
 
         if isinstance(volume, np.ndarray):
-            return _take(volume, index, axis=self.xyz2dims[axis])
+            return _take(volume, index, axis=axis)
         elif volume in 'annotation':
-            iregion = _take_remap(self.label, index, self.xyz2dims[axis], mapping)
+            iregion = _take_remap(self.label, index, axis=axis, mapping=mapping)
             return self._label2rgb(iregion)
         elif volume == 'image':
-            return _take(self.image, index, axis=self.xyz2dims[axis])
+            return _take(self.image, index, axis=axis)
         elif volume == 'value':
-            return region_values[_take_remap(self.label, index, self.xyz2dims[axis], mapping)]
+            return region_values[_take_remap(self.label, index, axis, mapping)]
         elif volume == 'image':
-            return _take(self.image, index, axis=self.xyz2dims[axis])
+            return _take(self.image, index, axis=axis)
         elif volume in ['surface', 'edges']:
             self.compute_surface()
-            return _take(self.surface, index, axis=self.xyz2dims[axis])
+            return _take(self.surface, index, axis=axis)
         elif volume == 'boundary':
-            iregion = _take_remap(self.label, index, self.xyz2dims[axis], mapping)
+            iregion = _take_remap(self.label, index, axis, mapping)
             return self.compute_boundaries(iregion)
 
         elif volume == 'volume':
             if bc is not None:
                 index = bc.xyz2i(np.array([coordinate] * 3))[axis]
-            return _take(region_values, index, axis=self.xyz2dims[axis])
+            return _take(region_values, index, axis=axis)
 
     def compute_boundaries(self, values):
         """
@@ -896,7 +905,6 @@ def plot_cslice(self, ap_coordinate, volume='image', mapping=None, region_values
         :param **kwargs: matplotlib.pyplot.imshow kwarg arguments
         :return: matplotlib ax object
         """
-
         cslice = self.slice(ap_coordinate, axis=1, volume=volume, mapping=mapping, region_values=region_values)
         return self._plot_slice(np.moveaxis(cslice, 0, 1), extent=self.extent(axis=1), volume=volume, **kwargs)
 
@@ -964,14 +972,11 @@ def plot_top(self, volume='annotation', mapping=None, region_values=None, ax=Non
         :param kwargs:
         :return:
         """
-
         self.compute_surface()
         ix, iy = np.meshgrid(np.arange(self.bc.nx), np.arange(self.bc.ny))
         iz = self.bc.z2i(self.top)
         inds = self._lookup_inds(np.stack((ix, iy, iz), axis=-1))
-
         regions = self._get_mapping(mapping=mapping)[self.label.flat[inds]]
-
         if volume == 'annotation':
             im = self._label2rgb(regions)
         elif volume == 'image':
@@ -1290,8 +1295,10 @@ class AllenAtlas(BrainAtlas):
 
     """numpy.array: A diffusion weighted imaging (DWI) image volume.
 
-    The Allen atlas DWI average template volume has with the shape (ap, ml, dv) and contains uint16
-    values.  FIXME What do the values represent?
+    This average template brain was created from images of 1,675 young adult C57BL/6J mouse brains
+    acquired using serial two-photon tomography.
+    This volume has a C-ordered shape of (ap, ml, dv) and contains uint16
+    values.
     """
     image = None