brain_loop_search.brain_utils

  1"""
  2
  3This module contains all utilities associated with brain projects, so other modules can be used to do without coupling
  4the brain ontology.
  5
  6The `CCFv3Ontology` is an example of how to derive an ontology.
  7
  8`draw_loop_in_ccf` is a visualization based on the brainrender project, also a good example.
  9What is tricky is that this lib needs internet connection to download the atlas, and used some deprecated functions
 10by other libs, so it may not work well for first-time. You might need to modify the lib's code.
 11running. (The usage of brainrender will block other processes from using the package, so it's dumped.)
 12
 13"""
 14# import itertools
 15# import os
 16import typing
 17# import networkx as nx
 18
 19import pandas as pd
 20import numpy as np
 21from importlib_resources import files, as_file
 22from .packing import Ontology
 23# from brainrender import Scene, settings
 24# from myterial import grey, white, grey_dark
 25# from vedo.shapes import Spline, Arrow, Tube, Line
 26
 27import colorsys
 28import random
 29# from brainrender._colors import map_color
 30
 31
 32with as_file(files('brain_loop_search') / 'structures.csv') as path:
 33    ccfv3 = pd.read_csv(path, index_col=1)
 34
 35
 36class CCFv3Ontology(Ontology):
 37    def __init__(self):
 38        """
 39        Use the already stashed ccfv3 structure table, containing the ID and acronym of each brain region.
 40        """
 41        super(Ontology, self).__init__()
 42        self._tree = ccfv3.copy()
 43        self._tree['parents'] = self._tree['structure_id_path']. \
 44            str.removeprefix('/').str.removesuffix('/').str.split('/')  # create ancestors as lists
 45        self._tree['parents'] = self._tree['parents'].apply(lambda x:[int(i) for i in x[:-1]])    # remove the last one(self)
 46        self._tree['depth'] = self._tree['parents'].apply(len)  # create depths
 47        self._find_children()
 48        self._leveling()
 49
 50    def _find_children(self):
 51        self._tree['children'] = [[] for i in range(len(self._tree))]
 52        for ind, row in self._tree.iterrows():
 53            id = row['parent_structure_id']
 54            if not np.isnan(id):
 55                self._tree.at[id,'children'].append(ind)
 56
 57    def _leveling(self):
 58        self._tree['level'] = 0
 59        for st in self._tree.index[self._tree['children'].apply(len) == 0]:
 60            par = self._tree.at[st, 'parent_structure_id']
 61            count = 1
 62            while par in self._tree.index and self._tree.at[par, 'level'] < count:
 63                self._tree.at[par, 'level'] = count
 64                par = self._tree.at[par, 'parent_structure_id']
 65                count += 1
 66
 67    def check_include(self, vert: typing.Iterable):
 68        return set(vert).issubset(self._tree.index)
 69
 70    def levels_of(self, vert: typing.Iterable) -> pd.Series:
 71        return self._tree.loc[vert, 'level']
 72
 73    def depths_of(self, vert: typing.Iterable) -> pd.Series:
 74        return self._tree.loc[vert, 'depth']
 75
 76    def ancestors_of(self, vert: typing.Iterable) -> pd.Series:
 77        return self._tree.loc[vert, 'parents']
 78
 79    def immediate_children_of(self, vert: typing.Iterable) -> pd.Series:
 80        return self._tree.loc[vert, 'children']
 81
 82
 83# def draw_brain_graph(graph: nx.DiGraph, path: str | os.PathLike, thr: float = 0, render_ops: dict = None, cmap='jet'):
 84#     """
 85#     Plot a directed graph among ccfv3 brain structures.
 86#
 87#     Note: It will reset some of the brainrender global settings.
 88#     Do not use this in interactive mode like jupyter, where some render options may not work.
 89#
 90#     Using brainrender can cause some problem, as it
 91#     will attempt downloading the brain atlas from the internet and ping google.com beforehand.
 92#
 93#     If you are not connected, it will fail.
 94#
 95#     One solution to do this is to remove the ping in the package source. If the downloading is too slow,
 96#     you can manually download it from their website and decompress it into the package's storage directory, usually set
 97#     in `$HOME/.brainglobe`. Anyway, you can hack into their code to see it for yourself.
 98#
 99#     Another problem I found with brainrender is that it uses numpy's deprecated features, you might need to
100#     refactor the lib source until it passes.
101#
102#     :param graph: a list of shortest paths consisting of brain structure IDs. Each sublist will be assigned a different
103#     random color. You need to make sure the heads and tails are repeated in adjacent lists.
104#     :param path: screenshot save path.
105#     :param thr: only edge weights over this will be plotted.
106#     :param render_ops: render options. Default is None to use the default options, see the code.
107#     :param cmap: matplotlib colormap.
108#     """
109#     if render_ops is None:
110#         render_ops = {
111#             'interactive': False,
112#             'camera': {
113#                 'pos': (4811, 3225, -42167),
114#                 'viewup': (0, -1, 0),
115#                 'clippingRange': (24770, 51413),
116#                 'focalPoint': (7252, 4096, -5657),
117#                 'distance': 36602
118#             },
119#             'zoom': 2
120#         }
121#     settings.SHOW_AXES = False
122#     settings.SHADER_STYLE = "cartoon"
123#     settings.ROOT_ALPHA = .05
124#     settings.ROOT_COLOR = grey
125#     settings.BACKGROUND_COLOR = white
126#
127#     scene = Scene(atlas_name='allen_mouse_100um')
128#
129#     regions = list(graph.nodes)
130#     regions = ccfv3.loc[regions, 'acronym']
131#     scene.add_brain_region(*list(regions), alpha=.02, hemisphere='left', silhouette=True)
132#     e = graph.edges(data=True)
133#     data = [d['weight'] for u, v, d in e if u != v]
134#     vmax, vmin = max(data), min(data)
135#
136#     # change silhouette
137#     for i in scene.get_actors(br_class="brain region"):
138#         i._silhouette_kwargs['lw'] = 1
139#         i._silhouette_kwargs['color'] = grey_dark
140#
141#     rt = scene.get_actors(br_class="brain region", name="root")[0]
142#     rt._silhouette_kwargs['lw'] = 1
143#     rt._silhouette_kwargs['color'] = grey
144#
145#     for u, v, d in e:
146#         if u == v or d['weight'] < thr:
147#             continue
148#         # tube
149#         # get a proper center of each region (this is difficult, for brain structures can be very twisted)
150#         # then connect them to make a spline for a tube, which will envelop arrows
151#         run = list(regions.loc[[u, v]])
152#         actors = scene.get_actors(br_class="brain region", name=run)
153#         sorted_actors = [None] * len(actors)
154#         run_map = dict(zip(run, range(len(run))))
155#         for a in actors:
156#             sorted_actors[run_map[a.name]] = a
157#         z_mean = [np.mean(m.points()[:, 2]) for m in sorted_actors]
158#         centers = [np.mean(m.points()[m.points()[:, 2] - z < 10], axis=0) * (1, 1, -1) for m, z in zip(sorted_actors, z_mean)]
159#         spl = Line(*centers, res=3)
160#         pts = spl.points()
161#         c = map_color(d['weight'], cmap, vmin, vmax)
162#         feint = list(colorsys.rgb_to_hsv(*c))
163#         feint[1] /= 2
164#         feint = colorsys.hsv_to_rgb(*feint)
165#
166#         rr = 1 + (d['weight'] - vmin) / (vmax - vmin)
167#         radius = [(np.linalg.norm(i - centers[0]) + np.linalg.norm(i - centers[-1])) / 20 * rr for i in pts]
168#         scene.add(Tube(pts, radius, c=feint, alpha=0.2))
169#
170#         # arrows
171#         scene.add(*[Arrow(*i, c=c, s=rr*3) for i in itertools.pairwise(pts)])
172#
173#     scene.render(**render_ops)
174#     scene.screenshot(str(path))
175#     scene.close()
176#
177#
178# def draw_single_loop(loop: list[list], path: str | os.PathLike, render_ops: dict = None):
179#     """
180#     Plot one loop in the ccfv3 atlas using brainrender.
181#
182#     Note: It will reset some of the brainrender global settings.
183#     Do not use this in interactive mode like jupyter, where some render options may not work.
184#
185#     Using brainrender can cause some problem, as it
186#     will attempt downloading the brain atlas from the internet and ping google.com beforehand.
187#
188#     If you are not connected, it will fail.
189#
190#     One solution to do this is to remove the ping in the package source. If the downloading is too slow,
191#     you can manually download it from their website and decompress it into the package's storage directory, usually set
192#     in `$HOME/.brainglobe`. Anyway, you can hack into their code to see it for yourself.
193#
194#     Another problem I found with brainrender is that it uses numpy's deprecated features, you might need to
195#     refactor the lib source until it passes.
196#
197#     :param loop: a list of shortest paths consisting of brain structure IDs. Each sublist will be assigned a different
198#     random color. You need to make sure the heads and tails are repeated in adjacent lists.
199#     :param path: screenshot save path.
200#     :param render_ops: render options. Default is None to use the default options, see the code.
201#     """
202#     if render_ops is None:
203#         render_ops = {
204#             'interactive': False,
205#             'camera': {
206#                 'pos': (4811, 3225, -42167),
207#                 'viewup': (0, -1, 0),
208#                 'clippingRange': (24770, 51413),
209#                 'focalPoint': (7252, 4096, -5657),
210#                 'distance': 36602
211#             },
212#             'zoom': 2
213#         }
214#     settings.SHOW_AXES = False
215#     settings.SHADER_STYLE = "cartoon"
216#     settings.ROOT_ALPHA = .05
217#     settings.ROOT_COLOR = grey
218#     settings.BACKGROUND_COLOR = white
219#
220#     scene = Scene(atlas_name='allen_mouse_100um')
221#
222#     # the root brain (usually this is just a background and not used to plot loops)
223#     rt = scene.get_actors(br_class="brain region", name="root")[0]
224#     rt._silhouette_kwargs['lw'] = 1
225#     rt._silhouette_kwargs['color'] = grey
226#
227#     text_map = {}
228#     count = 0
229#     for run in loop:
230#         run = list(ccfv3.loc[run, 'acronym'])
231#         for i in run[:-1]:
232#             if i not in text_map:
233#                 text_map[i] = []
234#             count += 1
235#             text_map[i].append(str(count))
236#
237#     for run in loop:    # each run in a loop is one sssp, will be marked by different colors
238#         # map to ccf acronym
239#         run = list(ccfv3.loc[run, 'acronym'])
240#
241#         # all traversed brain structures but the first and last one (axes)
242#         scene.add_brain_region(*run[1:-1], alpha=.2, hemisphere='left', silhouette=False)
243#
244#         # axis regions will add silhouette, and bigger alpha
245#         scene.add_brain_region(run[0], run[-1], alpha=.5, hemisphere='left', silhouette=False)
246#         scene.add_silhouette(*scene.get_actors(br_class="brain region", name=[run[0], run[-1]]), lw=2)
247#
248#         # random hue
249#         hue = random.random()
250#
251#         # tube
252#         # get a proper center of each region (this is difficult, for brain structures can be very twisted)
253#         # then connect them to make a spline for a tube, which will envelop arrows
254#         actors = scene.get_actors(br_class="brain region", name=run)
255#         sorted_actors = [None] * len(actors)
256#         run_map = dict(zip(run, range(len(run))))
257#         for a in actors:
258#             sorted_actors[run_map[a.name]] = a
259#         z_mean = [np.mean(m.points()[:, 2]) for m in sorted_actors]
260#         centers = [np.mean(m.points()[m.points()[:, 2] - z < 10], axis=0) * (1, 1, -1) for m, z in zip(sorted_actors, z_mean)]
261#         if len(centers) < 3:
262#             spl = Line(*centers, res=20)
263#         else:
264#             spl = Spline(centers)
265#         pts = spl.points()
266#         radius = [(np.linalg.norm(i - centers[0]) + np.linalg.norm(i - centers[-1])) / 100 for i in pts]
267#         scene.add(Tube(pts, radius, c=colorsys.hsv_to_rgb(hue, .5, .9), alpha=0.2))
268#
269#         # arrows
270#         scene.add(*[Arrow(*i, c=colorsys.hsv_to_rgb(hue, .9, .9)) for i in itertools.pairwise(pts)])
271#
272#         # text
273#         for i in range(len(run) - 1):
274#             if run[i] in text_map:
275#                 sorted_actors[i].caption(f'{"/".join(text_map.pop(run[i]))}. {run[i]}',
276#                                          centers[i] * (1, 1, -1), (.04, .04))
277#
278#     scene.render(**render_ops)
279#     scene.screenshot(str(path))
280#     scene.close()