Surfaces

CHANGELOG.md

@@ -1,3 +1,8 @@
+## [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
+- BrainAtlas.mask(): method that computes a mask of the convex brain volume shape
+
 ## [0.2.1]
 
 ### Modified

README.md

@@ -1,2 +1,20 @@
 # iblatlas
-IBL atlas module
+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.
+
+One of the requirment of this repository is to use minimal requirements, based of standard matplotlib, numpy and scipy libraries, and excludes visualisation tool.
+
+The documentation can be found here: 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
+
+Once merged the package is uploaded to PyPI using GitHub Actions.

iblatlas/__init__.py

@@ -194,4 +194,4 @@
    http://help.brain-map.org/download/attachments/8323525/Mouse_Common_Coordinate_Framework.pdf
 
 """
-__version__ = '0.2.1'
+__version__ = '0.3.0'

iblatlas/atlas.py

@@ -496,6 +496,18 @@ def _get_cache_dir():
         path_atlas = Path(par.CACHE_DIR).joinpath('histology', 'ATLAS', 'Needles', 'Allen', 'flatmaps')
         return path_atlas
 
+    def mask(self):
+        """
+        Returns a Boolean mask volume of the brain atlas, where 1 is inside the convex brain and 0 is outside
+        This returns an ovoid volume shaped like the brain and this will contain void values in the ventricules.
+        :return: np.array Bool (nap, nml, ndv)
+        """
+        self.compute_surface()
+        mask = np.logical_and(np.cumsum(self.label != 0, axis=-1) > 0,
+                              np.flip(np.cumsum(np.flip(self.label, axis=-1) != 0, axis=-1), axis=-1) > 0)
+        mask[np.isnan(self.top)] = 0
+        return mask
+
     def compute_surface(self):
         """
         Get the volume top, bottom, left and right surfaces, and from these the outer surface of
@@ -1197,28 +1209,35 @@ def tip(self):
         return sph2cart(- self.depth, self.theta, self.phi) + np.array((self.x, self.y, self.z))
 
     @staticmethod
-    def _get_surface_intersection(traj, brain_atlas, surface='top'):
+    def _get_surface_intersection(traj, brain_atlas, surface='top', mode='raise'):
         """
-        TODO Document!
+        Computes the intersection of a trajectory with either the top or the bottom surface of an atlas.
 
         Parameters
         ----------
-        traj
-        brain_atlas
-        surface
+        traj: iblatlas.atlas.Trajectory object
+        brain_atlas: iblatlas.atlas.BrainAtlas (or descendant) object
+        surface: str, optional (defaults to 'top') 'top' or 'bottom'
+        mode: str, optional (defaults to 'raise') 'raise' or 'none': raise an error if no intersection
+         with the brain surface is found otherwise returns None
 
         Returns
         -------
-
+        xyz: np.array, 3 elements, x, y, z coordinates of the intersection point with the surface
+             None if no intersection is found and mode is not set to 'raise'
         """
         brain_atlas.compute_surface()
-
         distance = traj.mindist(brain_atlas.srf_xyz)
         dist_sort = np.argsort(distance)
         # In some cases the nearest two intersection points are not the top and bottom of brain
         # So we find all intersection points that fall within one voxel and take the one with
         # highest dV to be entry and lowest dV to be exit
         idx_lim = np.sum(distance[dist_sort] * 1e6 < np.max(brain_atlas.res_um))
+        if idx_lim == 0:  # no intersection found
+            if mode == 'raise':
+                raise ValueError('No intersection found with brain surface')
+            else:
+                return
         dist_lim = dist_sort[0:idx_lim]
         z_val = brain_atlas.srf_xyz[dist_lim, 2]
         if surface == 'top':
@@ -1236,26 +1255,26 @@ def _get_surface_intersection(traj, brain_atlas, surface='top'):
         return xyz
 
     @staticmethod
-    def get_brain_exit(traj, brain_atlas):
+    def get_brain_exit(traj, brain_atlas, mode='raise'):
         """
         Given a Trajectory and a BrainAtlas object, computes the brain exit coordinate as the
         intersection of the trajectory and the brain surface (brain_atlas.surface)
         :param brain_atlas:
         :return: 3 element array x,y,z
         """
         # Find point where trajectory intersects with bottom of brain
-        return Insertion._get_surface_intersection(traj, brain_atlas, surface='bottom')
+        return Insertion._get_surface_intersection(traj, brain_atlas, surface='bottom', mode=mode)
 
     @staticmethod
-    def get_brain_entry(traj, brain_atlas):
+    def get_brain_entry(traj, brain_atlas, mode='raise'):
         """
         Given a Trajectory and a BrainAtlas object, computes the brain entry coordinate as the
         intersection of the trajectory and the brain surface (brain_atlas.surface)
         :param brain_atlas:
         :return: 3 element array x,y,z
         """
         # Find point where trajectory intersects with top of brain
-        return Insertion._get_surface_intersection(traj, brain_atlas, surface='top')
+        return Insertion._get_surface_intersection(traj, brain_atlas, surface='top', mode=mode)
 
 
 class AllenAtlas(BrainAtlas):

iblatlas/regions.py

@@ -675,26 +675,32 @@ def remap(self, region_ids, source_map='Allen', target_map='Beryl'):
         ----------
         region_ids : array_like of int
             The region IDs to remap.
-        source_map : str
-            The source map name, in `self.mappings`.
-        target_map : str
-            The target map name, in `self.mappings`.
+        source_map : str | np.array
+            if string, it should be a key existing in self.mappings such as 'Allen', 'Cosmos' or 'Beryl'
+            if np.array, it should be a [nregions,] array containing the mapped regions index such as
+            mapping[original_allen_region_index] = target_mapping_region_index
+        target_map : str | np.array
+            if string, it should be a key existing in self.mappings such as 'Allen', 'Cosmos' or 'Beryl'
+            if np.array, it should be a [nregions,] array containing the mapped regions index such as
+            mapping[original_allen_region_index] = target_mapping_region_index
 
         Returns
         -------
         numpy.array of int
             The input IDs mapped to `target_map`.
         """
         isnan = np.isnan(region_ids)
+        idx_source_map = self.mappings[source_map] if isinstance(source_map, str) else source_map
+        idx_target_map = self.mappings[target_map] if isinstance(target_map, str) else target_map
         if np.sum(isnan) > 0:
             # In case the user provides nans
             nan_loc = np.where(isnan)[0]
-            _, inds = ismember(region_ids[~isnan], self.id[self.mappings[source_map]])
-            mapped_ids = self.id[self.mappings[target_map][inds]].astype(float)
+            _, inds = ismember(region_ids[~isnan], self.id[idx_source_map])
+            mapped_ids = self.id[idx_target_map[inds]].astype(float)
             mapped_ids = np.insert(mapped_ids, nan_loc, np.full(nan_loc.shape, np.nan))
         else:
-            _, inds = ismember(region_ids, self.id[self.mappings[source_map]])
-            mapped_ids = self.id[self.mappings[target_map][inds]]
+            _, inds = ismember(region_ids, self.id[idx_source_map])
+            mapped_ids = self.id[idx_target_map[inds]]
 
         return mapped_ids
 

iblatlas/tests/test_atlas.py

@@ -94,14 +94,18 @@ def test_remap(self):
         atlas_id = np.array([463, 685])  # CA3 and PO
         cosmos_id = self.brs.remap(atlas_id, source_map='Allen', target_map='Cosmos')
         expected_cosmos_id = [1089, 549]  # HPF and TH
-        assert np.all(cosmos_id == expected_cosmos_id)
+        np.testing.assert_equal(cosmos_id, expected_cosmos_id)
 
         # Test remap when we have nans
         atlas_id = np.array([463, np.nan, 685])
         cosmos_id = self.brs.remap(atlas_id, source_map='Allen', target_map='Cosmos')
         expected_cosmos_id = np.array([1089, np.nan, 549], dtype=float)  # HPF and TH
         np.testing.assert_equal(cosmos_id, expected_cosmos_id)
 
+        # Test with providing the mapping array instead
+        cosmos_id = self.brs.remap(atlas_id, source_map=self.brs.mappings['Allen'], target_map=self.brs.mappings['Cosmos'])
+        np.testing.assert_equal(cosmos_id, expected_cosmos_id)
+
     def test_id2id(self):
         # Test remapping of atlas id to atlas id
         atlas_id = np.array([463, 685])
@@ -330,6 +334,11 @@ def test_compute_regions_volume(self):
         self.ba.compute_regions_volume()
         self.assertTrue(self.ba.regions.volume.shape == self.ba.regions.acronym.shape)
 
+    def test_compute_surfaces(self):
+        self.ba.compute_surface()
+        self.assertEqual(self.ba.surface.shape, self.ba.label.shape)
+        self.assertEqual(self.ba.surface.shape, self.ba.mask().shape)
+
 
 class TestAtlasPlots(unittest.TestCase):
 
@@ -538,6 +547,12 @@ def test_init_from_track(self):
         self.assertTrue(brain_entry[2] == brain_atlas.bc.i2z(100))
         brain_exit = insertion.get_brain_exit(insertion.trajectory, brain_atlas)
         self.assertTrue(brain_exit[2] == brain_atlas.bc.i2z(104))
+        # test getting no intersection with the brain surface
+        brain_atlas.srf_xyz *= np.NaN
+        with self.assertRaises(ValueError):
+            insertion.get_brain_entry(insertion.trajectory, brain_atlas)
+        self.assertIsNone(insertion.get_brain_entry(insertion.trajectory, brain_atlas, mode='none'))
+        self.assertIsNone(insertion.get_brain_exit(insertion.trajectory, brain_atlas, mode='none'))
 
 
 class TestTrajectory(unittest.TestCase):