Merge pull request #21 from int-brain-lab/depth_lookup

Depth lookup

Author: mayofaulkner

CHANGELOG.md

@@ -1,3 +1,9 @@
+## [0.7.0]
+### Added
+- `iblatlas.atlas.xyz_to_depth` a lookup from xyz coord to depth in cortex
+### Modified
+- Moved `iblatlas.plots.get_bc10` to `iblatlas.atlas.get_bc`
+
 ## [0.6.0]
 ### Added
 - `iblatlas.genomics.agea` option to load pre-processed volume

iblatlas/__init__.py

@@ -194,4 +194,4 @@
 .. [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.6.0'
+__version__ = '0.7.0'

iblatlas/atlas.py

@@ -1723,3 +1723,81 @@ def _read_volume(file_volume):
             raise ValueError(
                 f'"{file_volume.suffix}" files not supported, must be either ".nrrd" or ".npz"')
         return volume
+
+
+def get_bc(res_um=10):
+    """
+    Get BrainCoordinates object for Allen Atlas with bregma origin without reading in volumes
+
+    Parameters
+    ----------
+    res_um : float or int
+        The resolution of the brain coordinates object to return in um. Options are 50, 25 or 10
+
+    Returns
+    -------
+    BrainCoordinates object with bregma origin
+    """
+
+    dims2xyz = np.array([1, 0, 2])
+    scaling = np.array([1, 1, 1])
+    sf = 50 / res_um
+    im_shape = np.array([int(264 * sf), int(228 * sf), int(160 * sf)])
+    iorigin = (ALLEN_CCF_LANDMARKS_MLAPDV_UM['bregma'] / res_um)
+    dxyz = res_um * 1e-6 * np.array([1, -1, -1]) * scaling
+    nxyz = np.array(im_shape)[dims2xyz]
+    bc = BrainCoordinates(nxyz=nxyz, xyz0=(0, 0, 0), dxyz=dxyz)
+    bc = BrainCoordinates(nxyz=nxyz, xyz0=-bc.i2xyz(iorigin), dxyz=dxyz)
+
+    return bc
+
+
+def _download_depth_files(file_name):
+    """
+    These files have been generated using this script
+    https://github.com/int-brain-lab/ibldevtools/blob/master/Mayo/flatmaps/2025-03-20_depths_from_streamlines.py
+    :param file_name:
+    :return:
+    """
+    file_path = BrainAtlas._get_cache_dir().joinpath('depths', file_name)
+    if not file_path.exists():
+        file_path.parent.mkdir(exist_ok=True, parents=True)
+        aws.s3_download_file(f'atlas/depths/{file_path.name}', file_path)
+
+    return np.load(file_path)
+
+
+def xyz_to_depth(xyz, res_um=25):
+    """
+    For a given xyz coordinates find the depth from the surface of the cortex. Note the lookup will only
+    work for xyz cooordinates that are in the Isocortex of the Allen volume. If coordinates outside of this
+    region are given then the depth is returned as nan.
+
+    Parameters
+    ----------
+    xyz : numpy.array
+        An (n, 3) array of Cartesian coordinates. The order is ML, AP, DV and coordinates should be given in meters
+        relative to bregma.
+
+    res_um : float or int
+        The resolution of the brain atlas to do the depth lookup
+
+    Returns
+    -------
+        numpy.array
+        The depths from the surface of the cortex for each cartesian coordinate. If the coordinate does not lie within
+        the Isocortex, depth value returned is nan
+    """
+
+    ind_flat = _download_depth_files(f'depths_ind_{res_um}.npy')
+    depths = _download_depth_files(f'depths_um_{res_um}.npy')
+    bc = get_bc(res_um=res_um)
+
+    ixyz = bc.xyz2i(xyz, mode='clip')
+    iravel = np.ravel_multi_index((ixyz[:, 1], ixyz[:, 0], ixyz[:, 2]), (bc.ny, bc.nx, bc.nz))
+    a, b = ismember(ind_flat, iravel)
+
+    lookup_depths = np.full(iravel.shape, np.nan, dtype=np.float32)
+    lookup_depths[b] = depths[a]
+
+    return lookup_depths

iblatlas/genomics/gene_expression_scrapping/06-denoise-impute.py

@@ -78,7 +78,7 @@ def compute_agea_ppca(input_memmap, atlas_agea, output_memmap=None, ppca_n_compo
     gexps[gexps < 0] = np.nan
     C, ss, M, X, Ye = pyppca.ppca(gexps, d=ppca_n_components, dia=True)
     output_memmap = np.copy(input_memmap) if output_memmap is None else output_memmap
-    output_memmap[:, *np.unravel_index(inside_idx, input_memmap.shape[1:])] = Ye.T
+    output_memmap[:, *np.unravel_index(inside_idx, input_memmap.shape[1:])] = Ye.T  # noqa
     output_memmap[:, *outside_idx] = input_memmap[:, *outside_idx]
     return output_memmap
 

iblatlas/gui/atlasview.py

@@ -467,7 +467,7 @@ def cursor2xyamp(self, qpoint):
             region = None
         i = ctrl.atlas._lookup(xyz, mode='clip')
         vol = ctrl.atlas.image if isinstance(ctrl.volume, str) else ctrl.volume
-        v = vol[*np.unravel_index(i, vol.shape)]
+        v = vol[*np.unravel_index(i, vol.shape)] # noqa
         return iw, ih, w, h, v, region
 
 

iblatlas/plots.py

@@ -18,33 +18,11 @@
 
 from iblatlas.flatmaps import FlatMap, _swanson_labels_positions, swanson, swanson_json
 from iblatlas.regions import BrainRegions
-from iblatlas.atlas import AllenAtlas, BrainCoordinates, ALLEN_CCF_LANDMARKS_MLAPDV_UM
+from iblatlas.atlas import AllenAtlas, get_bc
 
 _logger = logging.getLogger(__name__)
 
 
-def get_bc_10():
-    """
-    Get BrainCoordinates object for 10um Allen Atlas
-
-    Returns
-    -------
-    BrainCoordinates object
-    """
-    dims2xyz = np.array([1, 0, 2])
-    res_um = 10
-    scaling = np.array([1, 1, 1])
-    image_10 = np.array([1320, 1140, 800])
-
-    iorigin = (ALLEN_CCF_LANDMARKS_MLAPDV_UM['bregma'] / res_um)
-    dxyz = res_um * 1e-6 * np.array([1, -1, -1]) * scaling
-    nxyz = np.array(image_10)[dims2xyz]
-    bc = BrainCoordinates(nxyz=nxyz, xyz0=(0, 0, 0), dxyz=dxyz)
-    bc = BrainCoordinates(nxyz=nxyz, xyz0=-bc.i2xyz(iorigin), dxyz=dxyz)
-
-    return bc
-
-
 def plot_polygon(ax, xy, color, reg_id, edgecolor='k', linewidth=0.3, alpha=1):
     """
     Function to plot matplotlib polygon on an axis
@@ -295,7 +273,7 @@ def _plot_slice_vector(coords, slice, values, mapping, empty_color='silver', mas
     if ba.res_um == 10:
         bc10 = ba.bc
     else:
-        bc10 = get_bc_10()
+        bc10 = get_bc(res_um=10)
 
     if ax is None:
         fig, ax = plt.subplots()

iblatlas/tests/test_atlas.py

@@ -1,10 +1,11 @@
 import unittest
+from unittest import mock
 
 import numpy as np
 import matplotlib.pyplot as plt
 
 from iblatlas.atlas import (BrainCoordinates, cart2sph, sph2cart, Trajectory, BrainAtlas,
-                            Insertion, ALLEN_CCF_LANDMARKS_MLAPDV_UM, AllenAtlas)
+                            Insertion, ALLEN_CCF_LANDMARKS_MLAPDV_UM, AllenAtlas, get_bc, xyz_to_depth)
 from iblatlas.regions import BrainRegions, FranklinPaxinosRegions
 from iblatlas.plots import prepare_lr_data, reorder_data
 from iblutil.numerical import ismember
@@ -729,5 +730,20 @@ def test_conversions(self):
         np.testing.assert_equal(self.frs.acronym2id('LO-2/3'), np.array([]))
 
 
+class TestDepths(unittest.TestCase):
+
+    @mock.patch('iblatlas.atlas._download_depth_files',
+                side_effect=[np.array([35838117, 41837813]), np.array([400, 600])])
+    def test_depth_lookup(self, mock_file):
+        bc = get_bc(25)
+        # First two coordinates in isocortex, third in thalamus
+        idx = np.array([[245, 274, 37], [286, 327, 53], [225, 287, 170]])  # ap, ml, dv
+        xyz = np.c_[bc.i2x(idx[:, 1]), bc.i2y(idx[:, 0]), bc.i2z(idx[:, 2])]
+
+        depths = xyz_to_depth(xyz)
+
+        assert np.array_equal(depths, np.array([400, 600, np.nan]), equal_nan=True)
+
+
 if __name__ == "__main__":
     unittest.main(exit=False, verbosity=2)