diff --git a/CHANGELOG.md b/CHANGELOG.md index 2762187..ebd49a7 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/README.md b/README.md index 576101c..0cc59d5 100644 --- a/README.md +++ b/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. diff --git a/iblatlas/__init__.py b/iblatlas/__init__.py index e19725c..b522ae2 100644 --- a/iblatlas/__init__.py +++ b/iblatlas/__init__.py @@ -194,4 +194,4 @@ http://help.brain-map.org/download/attachments/8323525/Mouse_Common_Coordinate_Framework.pdf """ -__version__ = '0.2.1' \ No newline at end of file +__version__ = '0.3.0' \ No newline at end of file diff --git a/iblatlas/atlas.py b/iblatlas/atlas.py index 8068f0a..2828c31 100644 --- a/iblatlas/atlas.py +++ b/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,7 +1255,7 @@ 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) @@ -1244,10 +1263,10 @@ def get_brain_exit(traj, 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) @@ -1255,7 +1274,7 @@ def get_brain_entry(traj, 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): diff --git a/iblatlas/regions.py b/iblatlas/regions.py index 77460cd..0907779 100644 --- a/iblatlas/regions.py +++ b/iblatlas/regions.py @@ -675,10 +675,14 @@ 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 ------- @@ -686,15 +690,17 @@ def remap(self, region_ids, source_map='Allen', target_map='Beryl'): 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 diff --git a/iblatlas/tests/test_atlas.py b/iblatlas/tests/test_atlas.py index 5bdd7d5..d8a242d 100644 --- a/iblatlas/tests/test_atlas.py +++ b/iblatlas/tests/test_atlas.py @@ -94,7 +94,7 @@ 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]) @@ -102,6 +102,10 @@ def test_remap(self): 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):