From f9c0d13930dccf073f8d0c41617341cd3cbd03cb Mon Sep 17 00:00:00 2001 From: jinhan Date: Wed, 11 Sep 2024 14:15:33 -0400 Subject: [PATCH 1/3] add h01 interface --- .../3_interfaces/plot_04_interface_h01.py | 148 +++++++ navis/interfaces/h01.py | 380 ++++++++++++++++++ 2 files changed, 528 insertions(+) create mode 100644 docs/examples/3_interfaces/plot_04_interface_h01.py create mode 100644 navis/interfaces/h01.py diff --git a/docs/examples/3_interfaces/plot_04_interface_h01.py b/docs/examples/3_interfaces/plot_04_interface_h01.py new file mode 100644 index 00000000..373ae9fd --- /dev/null +++ b/docs/examples/3_interfaces/plot_04_interface_h01.py @@ -0,0 +1,148 @@ +""" +H01 Dataset +============ + +In this notebook, you can learn how to work with the [H01 dataset](https://www.science.org/doi/10.1126/science.adk4858) using Navis. The H01 dataset contains 57,000 cells and 150 million synapses from a cubic millimeter of the human temporal cortex, which is [proofreading](https://h01-release.storage.googleapis.com/proofreading.html) in the CAVE. With this interface, you can access both a snapshot of the proofread dataset and the latest dataset using CAVEclient. +""" + + +# %% +## For Dev, clear cache and reload module +# import navis +# import importlib +# import navis.interfaces.h01 as h01 +# importlib.reload(h01) +import navis +import navis.interfaces.h01 as h01 + +# %% +## Examples in using Navis interface for H01 dataset + + +# %% +auth = h01.Authentication() + +# %% +## Scenario 1: Completely new user of H01 dataset + +# %% +auth.setup_token(make_new=True) + +# %% +auth.save_token(token=, overwrite=True) + +# %% +# %% +client = h01.get_cave_client() + +# %% +## Scenario 2: Existing user, new computer + +# %% +auth.setup_token(make_new=False) + +# %% +auth.save_token(token=, overwrite=True) + +# %% +client = h01.get_cave_client() + +# %% +## Get CAVE client after setting up your token in your computer + +# %% +client = h01.get_cave_client() + +# %% +## Query Tables + +# %% +client.materialize.get_versions() +print("Mat version: ", client.materialize.version) +client.materialize.get_tables() + +# %% +### Query Materialized Synapse Table + + +# %% +client.materialize.synapse_query(limit=10) + +# %% +syn = client.materialize.synapse_query( + post_ids=[864691131861340864], + # pre_ids=[ADD YOUR ROOT ID], +) +syn.head() +print(len(syn)) + +# %% +### Live Synapse Queries + + +# %% +import datetime +# check if root ID is the most recent root ID +root_id = 864691131861340864 +now = datetime.datetime.now(datetime.timezone.utc) +is_latest = client.chunkedgraph.is_latest_roots([root_id], timestamp=now) +latest_id = client.chunkedgraph.get_latest_roots(root_id, timestamp=now) +print(is_latest) +print(latest_id) + +# %% +synapse_table = client.info.get_datastack_info()['synapse_table'] +df=client.materialize.query_table(synapse_table, + timestamp = datetime.datetime.now(datetime.timezone.utc), + filter_equal_dict = {'post_pt_root_id': latest_id[0]}) +print(len(df)) +df.head() + +# %% +### Query Cells Table + +# %% +ct = client.materialize.query_table(table="cells") +ct.head() + +# %% +ct.cell_type.unique() + +# %% +### Filter by cell type + +# %% +interneuron_ids = ct[ct.cell_type == "INTERNEURON"].pt_root_id.values[:50] + +# %% +interneuron_ids = interneuron_ids[interneuron_ids != 0 ] +interneuron_ids + +# %% +len(interneuron_ids) + +# %% +interneurons = h01.fetch_neurons(interneuron_ids, lod=2, with_synapses=False) + +# %% +interneurons_ds = navis.simplify_mesh(interneurons, F=1/3) + +# %% +interneuron_ds + +# %% +import seaborn as sns +colors = {n.id: sns.color_palette('Reds', 7)[i] for i, n in enumerate(interneurons_ds)} +fig = navis.plot3d([interneurons_ds], color=colors) + +# %% +## Get Skeletons + +# %% +interneurons_sk = navis.skeletonize(interneurons, parallel=True) + +# %% +interneurons_sk + +# %% +fig = navis.plot3d([interneurons_sk[0], interneurons[0]], color=[(1, 0, 0), (1, 1, 1, .5)]) diff --git a/navis/interfaces/h01.py b/navis/interfaces/h01.py new file mode 100644 index 00000000..4024241b --- /dev/null +++ b/navis/interfaces/h01.py @@ -0,0 +1,380 @@ +from ..core import MeshNeuron, NeuronList +from .. import config, utils +import pandas as pd +import numpy as np +import warnings + +from concurrent.futures import ThreadPoolExecutor, as_completed +from functools import lru_cache +from textwrap import dedent + +import caveclient + +err_msg = dedent(""" + Failed to import `caveclient` library. Please install using pip: + + pip install caveclient -U + + """) + +try: + from caveclient import CAVEclient + import cloudvolume as cv +except ImportError: + logger.error(err_msg) + CAVEclient = None + cv = None +except BaseException: + raise + + +logger = config.get_logger(__name__) +dataset = None + +CAVE_DATASTACKS = { + 'h01_c3_flat': 'h01_c3_flat', +} + +SEG_URLS = { + 'h01_c3_flat': 'graphene://https://local.brain-wire-test.org/segmentation/table/h01_full0_v2' +} + +class Authentication: + def __init__(self): + self.url = "https://global.brain-wire-test.org/" + self.auth = caveclient.auth.AuthClient(server_address=self.url) + + def setup_token(self, make_new=True): + self.auth.setup_token(make_new) + + def save_token(self, token, overwrite): + self.auth.save_token(token, overwrite) + + +@lru_cache(None) +def get_cave_client(datastack='h01_c3_flat'): + """Get caveclient for given datastack. + + Parameters + ---------- + datastack : "h01_c3_flat" + Name of the dataset to use. + + """ + if not CAVEclient: + raise ImportError(err_msg) + + # Try mapping, else pass-through + datastack = CAVE_DATASTACKS.get(datastack, datastack) + print(CAVEclient(datastack)) + + return CAVEclient(datastack) + + +@lru_cache(None) +def get_cloudvol(url, cache=True): + """Get (cached) CloudVolume for given segmentation. + + Parameters + ---------- + url : str + + """ + if not cv: + raise ImportError(err_msg) + + return cv.CloudVolume(url, cache=cache, use_https=True, + progress=False, fill_missing=True) + + +def get_somas(root_ids, table='nucleus', datastack='h01_c3_flat'): + """Fetch somas based on nuclei segmentation for given neuron(s). + + Since this is a nucleus detection you will find that some neurons do + not have an entry despite clearly having a soma. This is due to the + "avocado problem" where the nucleus is separate from the rest of the + soma. + + Important + --------- + This data currently only exists for the 'cortex65' datastack (i.e. + "minnie65_public_v117"). + + Parameters + ---------- + root_ids : int | list of ints | None + Root ID(s) for which to fetch soma infos. Use + ``None`` to fetch complete list of annotated nuclei. + table : str + Which table to use for nucleus annotations. + datastack : "cortex65" | "cortex35" | "layer 2/3" + Which dataset to use. Internally these are mapped to the + corresponding sources (e.g. "minnie65_public_v117" for + "cortex65"). + + Returns + ------- + DataFrame + Pandas DataFrame with nuclei (see Examples). Root IDs + without a nucleus will simply not have an entry in the + table. + + """ + if datastack != 'h01_c3_flat': + warnings.warn('To our knowledge there is no nucleus segmentation ' + f'for "{datastack}". If that has changed please ' + 'get in touch on navis Github.') + + # Get/Initialize the CAVE client + client = get_cave_client(datastack) + + filter_in_dict = None + if not isinstance(root_ids, type(None)): + root_ids = utils.make_iterable(root_ids) + filter_in_dict = {'pt_root_id': root_ids} + + return client.materialize.query_table(table, filter_in_dict=filter_in_dict) + + +def fetch_neurons(x, *, lod=2, + with_synapses=True, + datastack='h01_c3_flat', + parallel=True, + max_threads=4, + **kwargs): + """Fetch neuron meshes. + + Notes + ----- + Synapses will be attached to the closest vertex on the mesh. + + Parameters + ---------- + x : str | int | list-like + Segment ID(s). Multiple Ids can be provided as list-like. + lod : int + Level of detail. Higher ``lod`` = coarser. This parameter + is ignored if the data source does not support multi-level + meshes. + with_synapses : bool, optional + If True will also attach synapses as ``.connectors``. + datastack : "cortex65" | "cortex35" | "layer 2/3" + Which dataset to use. Internally these are mapped to the + corresponding sources (e.g. "minnie65_public_v117" for + "cortex65"). + parallel : bool + If True, will use parallel threads to fetch data. + max_threads : int + Max number of parallel threads to use. + **kwargs + Keyword arguments are passed through to the initialization + of the ``navis.MeshNeurons``. + + Returns + ------- + navis.Neuronlist + Containing :class:`navis.MeshNeuron`. + + """ + x = utils.make_iterable(x, force_type=int) + client = get_cave_client(datastack) + print(client) + + if datastack in SEG_URLS: + url = SEG_URLS[datastack] + else: + url = client.info.get_datastack_info()['segmentation_source'] + vol = get_cloudvol(url) + + soma_pos = {} + + nl = [] + with ThreadPoolExecutor(max_workers=1 if not parallel else max_threads) as executor: + futures = {} + for id in x: + f = executor.submit(_fetch_single_neuron, + id, + vol=vol, + lod=lod, + client=client, + with_synapses=with_synapses, + source=datastack, + **kwargs + ) + futures[f] = id + + with config.tqdm(desc='Fetching', + total=len(x), + leave=config.pbar_leave, + disable=len(x) == 1 or config.pbar_hide) as pbar: + for f in as_completed(futures): + id = futures[f] + pbar.update(1) + try: + nl.append(f.result()) + except Exception as exc: + print(f'{id} generated an exception:', exc) + + nl = NeuronList(nl) + + for n in nl: + if n.id in soma_pos: + n.soma_pos = np.array(soma_pos[n.id]) * [8, 8, 33] + else: + n.soma_pos = None + + return nl + + +def _fetch_single_neuron(id, lod, vol, client, with_synapses=False, **kwargs): + """Fetch a single neuron.""" + # Make sure we only use `lod` if that's actually supported by the source + if 'MultiLevel' in str(type(vol.mesh)): + # mesh = vol.mesh.get(id, lod=lod, progress=False)[id] + mesh = vol.mesh.get(id, lod=lod)[id] + else: + # mesh = vol.mesh.get(id, progress=False)[id] + mesh = vol.mesh.get(id)[id] + + n = MeshNeuron(mesh, id=id, units='nm', **kwargs) + + if with_synapses: + pre = client.materialize.synapse_query(pre_ids=id) + post = client.materialize.synapse_query(post_ids=id) + + syn_table = client.materialize.synapse_table + syn_info = client.materialize.get_table_metadata(syn_table) + vxl_size = np.array(syn_info['voxel_resolution']).astype(int) + + to_concat = [] + if not pre.empty: + pre['type'] = 'pre' + locs = np.vstack(pre['pre_pt_position'].values) + locs = locs * vxl_size + pre['x'], pre['y'], pre['z'] = locs[:, 0], locs[:, 1], locs[:, 2] + pre = pre[['id', 'x', 'y', 'z', 'type', 'size']].copy() + to_concat.append(pre) + if not post.empty: + post['type'] = 'post' + locs = np.vstack(post['post_pt_position'].values) + locs = locs * vxl_size + post['x'], post['y'], post['z'] = locs[:, 0], locs[:, 1], locs[:, 2] + post = post[['id', 'x', 'y', 'z', 'type', 'size']].copy() + to_concat.append(post) + + if len(to_concat) == 1: + n.connectors = to_concat[0] + elif len(to_concat) == 2: + n.connectors = pd.concat(to_concat, axis=0).reset_index(drop=True) + + return n + + +def get_voxels(x, mip=0, bounds=None, datastack='h01_c3_flat'): + """Fetch voxels making a up given root ID. + + Parameters + ---------- + x : int + A single root ID. + mip : int + Scale at which to fetch voxels. + bounds : list, optional + Bounding box [xmin, xmax, ymin, ymax, zmin, zmax] in voxel + space. For example, the voxel resolution for mip 0 + segmentation is 8 x 8 x 40 nm. + datastack : "h01_c3_flat" + Which dataset to use. Internally these are mapped to the + corresponding sources (e.g. "h01_c3_flat" for + "h01_c3_flat"). + + Returns + ------- + voxels : (N, 3) np.ndarray + In voxel space according to `mip`. + + """ + client = get_cave_client(datastack) + # Need to get the graphene (not the precomputed) version of the data + vol_graphene = cv.CloudVolume(client.chunkedgraph.cloudvolume_path, + use_https=True, progress=False) + if datastack in SEG_URLS: + url = SEG_URLS[datastack] + else: + url = client.info.get_datastack_info()['segmentation_source'] + vol_prec = get_cloudvol(url) + + # Get L2 chunks making up this neuron + l2_ids = client.chunkedgraph.get_leaves(x, stop_layer=2) + + # Turn l2_ids into chunk indices + l2_ix = [np.array(vol_graphene.mesh.meta.meta.decode_chunk_position(l)) for l in l2_ids] + l2_ix = np.unique(l2_ix, axis=0) + + # Convert to nm + l2_nm = np.asarray(_chunks_to_nm(l2_ix, vol=vol_graphene)) + + # Convert back to voxel space (according to mip) + l2_vxl = l2_nm // vol_prec.meta.scales[mip]["resolution"] + + voxels = [] + ch_size = np.array(vol_graphene.mesh.meta.meta.graph_chunk_size) + ch_size = ch_size // (vol_prec.mip_resolution(mip) / vol_prec.mip_resolution(0)) + ch_size = np.asarray(ch_size).astype(int) + old_mip = vol_prec.mip + + if not isinstance(bounds, type(None)): + bounds = np.asarray(bounds) + if not bounds.ndim == 1 or len(bounds) != 6: + raise ValueError('`bounds` must be [xmin, xmax, ymin, ymax, zmin, zmax]') + l2_vxl = l2_vxl[np.all(l2_vxl >= bounds[::2], axis=1)] + l2_vxl = l2_vxl[np.all(l2_vxl < bounds[1::2] + ch_size, axis=1)] + + try: + vol_prec.mip = mip + for ch in config.tqdm(l2_vxl, desc='Loading'): + ct = vol_prec[ch[0]:ch[0] + ch_size[0], + ch[1]:ch[1] + ch_size[1], + ch[2]:ch[2] + ch_size[2]][:, :, :, 0] + this_vxl = np.dstack(np.where(ct == x))[0] + this_vxl = this_vxl + ch + voxels.append(this_vxl) + except BaseException: + raise + finally: + vol_prec.mip = old_mip + voxels = np.vstack(voxels) + + if not isinstance(bounds, type(None)): + voxels = voxels[np.all(voxels >= bounds[::2], axis=1)] + voxels = voxels[np.all(voxels < bounds[1::2], axis=1)] + + return voxels + + +def _chunks_to_nm(xyz_ch, vol, voxel_resolution=[4, 4, 40]): + """Map a chunk location to Euclidean space. + + Parameters + ---------- + xyz_ch : array-like + (N, 3) array of chunk indices. + vol : cloudvolume.CloudVolume + CloudVolume object associated with the chunked space. + voxel_resolution : list, optional + Voxel resolution. + + Returns + ------- + np.array + (N, 3) array of spatial points. + + """ + mip_scaling = vol.mip_resolution(0) // np.array(voxel_resolution, dtype=int) + + x_vox = np.atleast_2d(xyz_ch) * vol.mesh.meta.meta.graph_chunk_size + return ( + (x_vox + np.array(vol.mesh.meta.meta.voxel_offset(0))) + * voxel_resolution + * mip_scaling + ) From 25f8fded89c4c922aaf4b088d1fff6a92cfadcc4 Mon Sep 17 00:00:00 2001 From: jinhan Date: Tue, 17 Sep 2024 11:27:33 -0400 Subject: [PATCH 2/3] refactoring h01.py and microns.py / add example notebook for h01 interface --- .../plot_04_remote_h01.py} | 36 -- navis/interfaces/cave_utils.py | 356 ++++++++++++++++ navis/interfaces/h01.py | 384 +----------------- navis/interfaces/microns.py | 361 ++-------------- 4 files changed, 405 insertions(+), 732 deletions(-) rename docs/examples/{3_interfaces/plot_04_interface_h01.py => 4_remote/plot_04_remote_h01.py} (81%) create mode 100644 navis/interfaces/cave_utils.py diff --git a/docs/examples/3_interfaces/plot_04_interface_h01.py b/docs/examples/4_remote/plot_04_remote_h01.py similarity index 81% rename from docs/examples/3_interfaces/plot_04_interface_h01.py rename to docs/examples/4_remote/plot_04_remote_h01.py index 373ae9fd..9f8c15f1 100644 --- a/docs/examples/3_interfaces/plot_04_interface_h01.py +++ b/docs/examples/4_remote/plot_04_remote_h01.py @@ -7,51 +7,15 @@ # %% -## For Dev, clear cache and reload module -# import navis -# import importlib -# import navis.interfaces.h01 as h01 -# importlib.reload(h01) import navis import navis.interfaces.h01 as h01 # %% ## Examples in using Navis interface for H01 dataset - -# %% -auth = h01.Authentication() - -# %% -## Scenario 1: Completely new user of H01 dataset - -# %% -auth.setup_token(make_new=True) - -# %% -auth.save_token(token=, overwrite=True) - -# %% # %% client = h01.get_cave_client() -# %% -## Scenario 2: Existing user, new computer - -# %% -auth.setup_token(make_new=False) - -# %% -auth.save_token(token=, overwrite=True) - -# %% -client = h01.get_cave_client() - -# %% -## Get CAVE client after setting up your token in your computer - -# %% -client = h01.get_cave_client() # %% ## Query Tables diff --git a/navis/interfaces/cave_utils.py b/navis/interfaces/cave_utils.py new file mode 100644 index 00000000..9a4534d9 --- /dev/null +++ b/navis/interfaces/cave_utils.py @@ -0,0 +1,356 @@ +from ..core import MeshNeuron, NeuronList +from .. import config, utils +import pandas as pd +import numpy as np +import warnings + +from concurrent.futures import ThreadPoolExecutor, as_completed +from functools import lru_cache +from textwrap import dedent + +err_msg = dedent(""" + Failed to import `caveclient` library. Please install using pip: + + pip install caveclient -U + + """) + +try: + from caveclient import CAVEclient + import cloudvolume as cv +except ImportError: + config.logger.error(err_msg) + CAVEclient = None + cv = None +except BaseException: + raise + + +logger = config.get_logger(__name__) +dataset = None + +@lru_cache(None) +def get_cave_client(datastack="cortex65", server_address=None): + """Get caveclient for given datastack. + + Parameters + ---------- + datastack : "cortex65" | "cortex35" | "layer 2/3" | str + Which dataset to query. "cortex65", "cortex35" and "layer 2/3" + are internally mapped to the corresponding sources: for example, + "minnie65_public_vXXX" for "cortex65" where XXX is always the + most recent version). + + """ + if not CAVEclient: + raise ImportError(err_msg) + + return CAVEclient(datastack, server_address) + +@lru_cache(None) +def _get_cloudvol(url, cache=True): + """Get (cached) CloudVolume for given segmentation. + + Parameters + ---------- + url : str + + """ + if not cv: + raise ImportError(err_msg) + + return cv.CloudVolume( + url, cache=cache, use_https=True, parallel=10, progress=False, fill_missing=True + ) + +def _get_somas(root_ids, table="nucleus_detection_v0", datastack="cortex65"): + """Fetch somas based on nuclei segmentation for given neuron(s). + + Since this is a nucleus detection you will find that some neurons do + not have an entry despite clearly having a soma. This is due to the + "avocado problem" where the nucleus is separate from the rest of the + soma. + + Important + --------- + This data currently only exists for the 'cortex65' datastack (i.e. + "minnie65_public_vXXX"). + + Parameters + ---------- + root_ids : int | list of ints | None + Root ID(s) for which to fetch soma infos. Use + `None` to fetch complete list of annotated nuclei. + table : str + Which table to use for nucleus annotations. Also + see the `microns.list_annotation_tables()` function. + datastack : "cortex65" | "cortex35" | "layer 2/3" | str + Which dataset to query. "cortex65", "cortex35" and "layer 2/3" + are internally mapped to the corresponding sources: for example, + "minnie65_public_vXXX" for "cortex65" where XXX is always the + most recent version). + + Returns + ------- + DataFrame + Pandas DataFrame with nuclei (see Examples). Root IDs + without a nucleus will simply not have an entry in the + table. + + """ + if datastack != "cortex65": + warnings.warn( + "To our knowledge there is no nucleus segmentation " + f'for "{datastack}". If that has changed please ' + "get in touch on navis' Github." + ) + + # Get/Initialize the CAVE client + client = get_cave_client(datastack) + + filter_in_dict = None + if not isinstance(root_ids, type(None)): + root_ids = utils.make_iterable(root_ids) + filter_in_dict = {"pt_root_id": root_ids} + + return client.materialize.query_table(table, filter_in_dict=filter_in_dict) + +def fetch_neurons(x, lod, + with_synapses, + datastack, + parallel, + max_threads, + **kwargs): + """Fetch neuron meshes. + + Notes + ----- + Synapses will be attached to the closest vertex on the mesh. + + Parameters + ---------- + x : str | int | list-like + Segment ID(s). Multiple Ids can be provided as list-like. + lod : int + Level of detail. Higher ``lod`` = coarser. This parameter + is ignored if the data source does not support multi-level + meshes. + with_synapses : bool, optional + If True will also attach synapses as ``.connectors``. + datastack : "cortex65" | "cortex35" | "layer 2/3" + Which dataset to use. Internally these are mapped to the + corresponding sources (e.g. "minnie65_public_v117" for + "cortex65"). + parallel : bool + If True, will use parallel threads to fetch data. + max_threads : int + Max number of parallel threads to use. + **kwargs + Keyword arguments are passed through to the initialization + of the ``navis.MeshNeurons``. + + Returns + ------- + navis.Neuronlist + Containing :class:`navis.MeshNeuron`. + + """ + x = utils.make_iterable(x, force_type=int) + client = get_cave_client(datastack) + + vol = _get_cloudvol(client.info.segmentation_source()) # this is cached + + if datastack == "cortex65": + try: + somas = _get_somas(x, datastack=datastack) + soma_pos = somas.set_index("pt_root_id").pt_position.to_dict() + except BaseException as e: + logger.warning("Failed to fetch somas via nucleus segmentation" f"(){e})") + soma_pos = {} + else: + soma_pos = {} + + nl = [] + with ThreadPoolExecutor(max_workers=1 if not parallel else max_threads) as executor: + futures = {} + for id in x: + f = executor.submit(_fetch_single_neuron, + id, + vol=vol, + lod=lod, + client=client, + with_synapses=with_synapses, + source=datastack, + **kwargs + ) + futures[f] = id + + with config.tqdm(desc='Fetching', + total=len(x), + leave=config.pbar_leave, + disable=len(x) == 1 or config.pbar_hide) as pbar: + for f in as_completed(futures): + id = futures[f] + pbar.update(1) + try: + nl.append(f.result()) + except Exception as exc: + print(f'{id} generated an exception:', exc) + + nl = NeuronList(nl) + + for n in nl: + if n.id in soma_pos: + n.soma_pos = np.array(soma_pos[n.id]) * [8, 8, 33] + else: + n.soma_pos = None + + return nl + +def _fetch_single_neuron(id, lod, vol, client, with_synapses=False, **kwargs): + """Fetch a single neuron.""" + # Make sure we only use `lod` if that's actually supported by the source + if "MultiLevel" in str(type(vol.mesh)): + mesh = vol.mesh.get(id, lod=lod)[id] + else: + mesh = vol.mesh.get(id, deduplicate_chunk_boundaries=False)[id] + + n = MeshNeuron(mesh, id=id, units="nm", **kwargs) + + if with_synapses: + pre = client.materialize.synapse_query(pre_ids=id) + post = client.materialize.synapse_query(post_ids=id) + + syn_table = client.materialize.synapse_table + syn_info = client.materialize.get_table_metadata(syn_table) + vxl_size = np.array(syn_info["voxel_resolution"]).astype(int) + + to_concat = [] + if not pre.empty: + pre["type"] = "pre" + locs = np.vstack(pre["pre_pt_position"].values) + locs = locs * vxl_size + pre["x"], pre["y"], pre["z"] = locs[:, 0], locs[:, 1], locs[:, 2] + pre = pre[["id", "x", "y", "z", "type", "size"]].copy() + to_concat.append(pre) + if not post.empty: + post["type"] = "post" + locs = np.vstack(post["post_pt_position"].values) + locs = locs * vxl_size + post["x"], post["y"], post["z"] = locs[:, 0], locs[:, 1], locs[:, 2] + post = post[["id", "x", "y", "z", "type", "size"]].copy() + to_concat.append(post) + + if len(to_concat) == 1: + n.connectors = to_concat[0] + elif len(to_concat) == 2: + n.connectors = pd.concat(to_concat, axis=0).reset_index(drop=True) + + return n + + +def get_voxels(x, mip, bounds, datastack): + """Fetch voxels making a up given root ID. + + Parameters + ---------- + x : int + A single root ID. + mip : int + Scale at which to fetch voxels. + bounds : list, optional + Bounding box [xmin, xmax, ymin, ymax, zmin, zmax] in voxel + space. For example, the voxel resolution for mip 0 + segmentation is 8 x 8 x 40 nm. + datastack : "h01_c3_flat" + Which dataset to use. Internally these are mapped to the + corresponding sources (e.g. "h01_c3_flat" for + "h01_c3_flat"). + + Returns + ------- + voxels : (N, 3) np.ndarray + In voxel space according to `mip`. + + """ + client = get_cave_client(datastack) + # Need to get the graphene (not the precomputed) version of the data + vol_graphene = cv.CloudVolume(client.chunkedgraph.cloudvolume_path, + use_https=True, progress=False) + url = client.info.get_datastack_info()['segmentation_source'] + vol_prec = _get_cloudvol(url) + + # Get L2 chunks making up this neuron + l2_ids = client.chunkedgraph.get_leaves(x, stop_layer=2) + + # Turn l2_ids into chunk indices + l2_ix = [np.array(vol_graphene.mesh.meta.meta.decode_chunk_position(l)) for l in l2_ids] + l2_ix = np.unique(l2_ix, axis=0) + + # Convert to nm + l2_nm = np.asarray(_chunks_to_nm(l2_ix, vol=vol_graphene)) + + # Convert back to voxel space (according to mip) + l2_vxl = l2_nm // vol_prec.meta.scales[mip]["resolution"] + + voxels = [] + ch_size = np.array(vol_graphene.mesh.meta.meta.graph_chunk_size) + ch_size = ch_size // (vol_prec.mip_resolution(mip) / vol_prec.mip_resolution(0)) + ch_size = np.asarray(ch_size).astype(int) + old_mip = vol_prec.mip + + if not isinstance(bounds, type(None)): + bounds = np.asarray(bounds) + if not bounds.ndim == 1 or len(bounds) != 6: + raise ValueError('`bounds` must be [xmin, xmax, ymin, ymax, zmin, zmax]') + l2_vxl = l2_vxl[np.all(l2_vxl >= bounds[::2], axis=1)] + l2_vxl = l2_vxl[np.all(l2_vxl < bounds[1::2] + ch_size, axis=1)] + + try: + vol_prec.mip = mip + for ch in config.tqdm(l2_vxl, desc='Loading'): + ct = vol_prec[ch[0]:ch[0] + ch_size[0], + ch[1]:ch[1] + ch_size[1], + ch[2]:ch[2] + ch_size[2]][:, :, :, 0] + this_vxl = np.dstack(np.where(ct == x))[0] + this_vxl = this_vxl + ch + voxels.append(this_vxl) + except BaseException: + raise + finally: + vol_prec.mip = old_mip + voxels = np.vstack(voxels) + + if not isinstance(bounds, type(None)): + voxels = voxels[np.all(voxels >= bounds[::2], axis=1)] + voxels = voxels[np.all(voxels < bounds[1::2], axis=1)] + + return voxels + + +def _chunks_to_nm(xyz_ch, vol, voxel_resolution=[4, 4, 40]): + """Map a chunk location to Euclidean space. + + Parameters + ---------- + xyz_ch : array-like + (N, 3) array of chunk indices. + vol : cloudvolume.CloudVolume + CloudVolume object associated with the chunked space. + voxel_resolution : list, optional + Voxel resolution. + + Returns + ------- + np.array + (N, 3) array of spatial points. + + """ + mip_scaling = vol.mip_resolution(0) // np.array(voxel_resolution, dtype=int) + + x_vox = np.atleast_2d(xyz_ch) * vol.mesh.meta.meta.graph_chunk_size + return ( + (x_vox + np.array(vol.mesh.meta.meta.voxel_offset(0))) + * voxel_resolution + * mip_scaling + ) diff --git a/navis/interfaces/h01.py b/navis/interfaces/h01.py index 4024241b..06715d59 100644 --- a/navis/interfaces/h01.py +++ b/navis/interfaces/h01.py @@ -1,140 +1,14 @@ -from ..core import MeshNeuron, NeuronList -from .. import config, utils -import pandas as pd -import numpy as np -import warnings - -from concurrent.futures import ThreadPoolExecutor, as_completed from functools import lru_cache -from textwrap import dedent - -import caveclient - -err_msg = dedent(""" - Failed to import `caveclient` library. Please install using pip: - - pip install caveclient -U - - """) - -try: - from caveclient import CAVEclient - import cloudvolume as cv -except ImportError: - logger.error(err_msg) - CAVEclient = None - cv = None -except BaseException: - raise - - -logger = config.get_logger(__name__) -dataset = None - -CAVE_DATASTACKS = { - 'h01_c3_flat': 'h01_c3_flat', -} - -SEG_URLS = { - 'h01_c3_flat': 'graphene://https://local.brain-wire-test.org/segmentation/table/h01_full0_v2' -} - -class Authentication: - def __init__(self): - self.url = "https://global.brain-wire-test.org/" - self.auth = caveclient.auth.AuthClient(server_address=self.url) - - def setup_token(self, make_new=True): - self.auth.setup_token(make_new) - - def save_token(self, token, overwrite): - self.auth.save_token(token, overwrite) - +from . import cave_utils @lru_cache(None) -def get_cave_client(datastack='h01_c3_flat'): - """Get caveclient for given datastack. - - Parameters - ---------- - datastack : "h01_c3_flat" - Name of the dataset to use. - +def get_cave_client(): + """Get caveclient for H01 dataset. """ - if not CAVEclient: - raise ImportError(err_msg) - - # Try mapping, else pass-through - datastack = CAVE_DATASTACKS.get(datastack, datastack) - print(CAVEclient(datastack)) - - return CAVEclient(datastack) - - -@lru_cache(None) -def get_cloudvol(url, cache=True): - """Get (cached) CloudVolume for given segmentation. - - Parameters - ---------- - url : str - - """ - if not cv: - raise ImportError(err_msg) - - return cv.CloudVolume(url, cache=cache, use_https=True, - progress=False, fill_missing=True) - - -def get_somas(root_ids, table='nucleus', datastack='h01_c3_flat'): - """Fetch somas based on nuclei segmentation for given neuron(s). - - Since this is a nucleus detection you will find that some neurons do - not have an entry despite clearly having a soma. This is due to the - "avocado problem" where the nucleus is separate from the rest of the - soma. - - Important - --------- - This data currently only exists for the 'cortex65' datastack (i.e. - "minnie65_public_v117"). - - Parameters - ---------- - root_ids : int | list of ints | None - Root ID(s) for which to fetch soma infos. Use - ``None`` to fetch complete list of annotated nuclei. - table : str - Which table to use for nucleus annotations. - datastack : "cortex65" | "cortex35" | "layer 2/3" - Which dataset to use. Internally these are mapped to the - corresponding sources (e.g. "minnie65_public_v117" for - "cortex65"). - - Returns - ------- - DataFrame - Pandas DataFrame with nuclei (see Examples). Root IDs - without a nucleus will simply not have an entry in the - table. - - """ - if datastack != 'h01_c3_flat': - warnings.warn('To our knowledge there is no nucleus segmentation ' - f'for "{datastack}". If that has changed please ' - 'get in touch on navis Github.') - - # Get/Initialize the CAVE client - client = get_cave_client(datastack) - - filter_in_dict = None - if not isinstance(root_ids, type(None)): - root_ids = utils.make_iterable(root_ids) - filter_in_dict = {'pt_root_id': root_ids} - - return client.materialize.query_table(table, filter_in_dict=filter_in_dict) - + datastack = "h01_c3_flat" + server_address = "https://global.brain-wire-test.org/" + + return cave_utils.get_cave_client(datastack, server_address) def fetch_neurons(x, *, lod=2, with_synapses=True, @@ -143,238 +17,24 @@ def fetch_neurons(x, *, lod=2, max_threads=4, **kwargs): """Fetch neuron meshes. - - Notes - ----- - Synapses will be attached to the closest vertex on the mesh. - - Parameters - ---------- - x : str | int | list-like - Segment ID(s). Multiple Ids can be provided as list-like. - lod : int - Level of detail. Higher ``lod`` = coarser. This parameter - is ignored if the data source does not support multi-level - meshes. - with_synapses : bool, optional - If True will also attach synapses as ``.connectors``. - datastack : "cortex65" | "cortex35" | "layer 2/3" - Which dataset to use. Internally these are mapped to the - corresponding sources (e.g. "minnie65_public_v117" for - "cortex65"). - parallel : bool - If True, will use parallel threads to fetch data. - max_threads : int - Max number of parallel threads to use. - **kwargs - Keyword arguments are passed through to the initialization - of the ``navis.MeshNeurons``. - - Returns - ------- - navis.Neuronlist - Containing :class:`navis.MeshNeuron`. - """ - x = utils.make_iterable(x, force_type=int) - client = get_cave_client(datastack) - print(client) - - if datastack in SEG_URLS: - url = SEG_URLS[datastack] - else: - url = client.info.get_datastack_info()['segmentation_source'] - vol = get_cloudvol(url) - - soma_pos = {} - - nl = [] - with ThreadPoolExecutor(max_workers=1 if not parallel else max_threads) as executor: - futures = {} - for id in x: - f = executor.submit(_fetch_single_neuron, - id, - vol=vol, - lod=lod, - client=client, - with_synapses=with_synapses, - source=datastack, - **kwargs - ) - futures[f] = id - - with config.tqdm(desc='Fetching', - total=len(x), - leave=config.pbar_leave, - disable=len(x) == 1 or config.pbar_hide) as pbar: - for f in as_completed(futures): - id = futures[f] - pbar.update(1) - try: - nl.append(f.result()) - except Exception as exc: - print(f'{id} generated an exception:', exc) - - nl = NeuronList(nl) - - for n in nl: - if n.id in soma_pos: - n.soma_pos = np.array(soma_pos[n.id]) * [8, 8, 33] - else: - n.soma_pos = None - - return nl - - -def _fetch_single_neuron(id, lod, vol, client, with_synapses=False, **kwargs): - """Fetch a single neuron.""" - # Make sure we only use `lod` if that's actually supported by the source - if 'MultiLevel' in str(type(vol.mesh)): - # mesh = vol.mesh.get(id, lod=lod, progress=False)[id] - mesh = vol.mesh.get(id, lod=lod)[id] - else: - # mesh = vol.mesh.get(id, progress=False)[id] - mesh = vol.mesh.get(id)[id] - - n = MeshNeuron(mesh, id=id, units='nm', **kwargs) - - if with_synapses: - pre = client.materialize.synapse_query(pre_ids=id) - post = client.materialize.synapse_query(post_ids=id) - - syn_table = client.materialize.synapse_table - syn_info = client.materialize.get_table_metadata(syn_table) - vxl_size = np.array(syn_info['voxel_resolution']).astype(int) - - to_concat = [] - if not pre.empty: - pre['type'] = 'pre' - locs = np.vstack(pre['pre_pt_position'].values) - locs = locs * vxl_size - pre['x'], pre['y'], pre['z'] = locs[:, 0], locs[:, 1], locs[:, 2] - pre = pre[['id', 'x', 'y', 'z', 'type', 'size']].copy() - to_concat.append(pre) - if not post.empty: - post['type'] = 'post' - locs = np.vstack(post['post_pt_position'].values) - locs = locs * vxl_size - post['x'], post['y'], post['z'] = locs[:, 0], locs[:, 1], locs[:, 2] - post = post[['id', 'x', 'y', 'z', 'type', 'size']].copy() - to_concat.append(post) - - if len(to_concat) == 1: - n.connectors = to_concat[0] - elif len(to_concat) == 2: - n.connectors = pd.concat(to_concat, axis=0).reset_index(drop=True) - - return n + return cave_utils.fetch_neurons( + x, + lod=lod, + with_synapses=with_synapses, + datastack=datastack, + parallel=parallel, + max_threads=max_threads, + **kwargs + ) def get_voxels(x, mip=0, bounds=None, datastack='h01_c3_flat'): """Fetch voxels making a up given root ID. - - Parameters - ---------- - x : int - A single root ID. - mip : int - Scale at which to fetch voxels. - bounds : list, optional - Bounding box [xmin, xmax, ymin, ymax, zmin, zmax] in voxel - space. For example, the voxel resolution for mip 0 - segmentation is 8 x 8 x 40 nm. - datastack : "h01_c3_flat" - Which dataset to use. Internally these are mapped to the - corresponding sources (e.g. "h01_c3_flat" for - "h01_c3_flat"). - - Returns - ------- - voxels : (N, 3) np.ndarray - In voxel space according to `mip`. - """ - client = get_cave_client(datastack) - # Need to get the graphene (not the precomputed) version of the data - vol_graphene = cv.CloudVolume(client.chunkedgraph.cloudvolume_path, - use_https=True, progress=False) - if datastack in SEG_URLS: - url = SEG_URLS[datastack] - else: - url = client.info.get_datastack_info()['segmentation_source'] - vol_prec = get_cloudvol(url) - - # Get L2 chunks making up this neuron - l2_ids = client.chunkedgraph.get_leaves(x, stop_layer=2) - - # Turn l2_ids into chunk indices - l2_ix = [np.array(vol_graphene.mesh.meta.meta.decode_chunk_position(l)) for l in l2_ids] - l2_ix = np.unique(l2_ix, axis=0) - - # Convert to nm - l2_nm = np.asarray(_chunks_to_nm(l2_ix, vol=vol_graphene)) - - # Convert back to voxel space (according to mip) - l2_vxl = l2_nm // vol_prec.meta.scales[mip]["resolution"] - - voxels = [] - ch_size = np.array(vol_graphene.mesh.meta.meta.graph_chunk_size) - ch_size = ch_size // (vol_prec.mip_resolution(mip) / vol_prec.mip_resolution(0)) - ch_size = np.asarray(ch_size).astype(int) - old_mip = vol_prec.mip - - if not isinstance(bounds, type(None)): - bounds = np.asarray(bounds) - if not bounds.ndim == 1 or len(bounds) != 6: - raise ValueError('`bounds` must be [xmin, xmax, ymin, ymax, zmin, zmax]') - l2_vxl = l2_vxl[np.all(l2_vxl >= bounds[::2], axis=1)] - l2_vxl = l2_vxl[np.all(l2_vxl < bounds[1::2] + ch_size, axis=1)] - - try: - vol_prec.mip = mip - for ch in config.tqdm(l2_vxl, desc='Loading'): - ct = vol_prec[ch[0]:ch[0] + ch_size[0], - ch[1]:ch[1] + ch_size[1], - ch[2]:ch[2] + ch_size[2]][:, :, :, 0] - this_vxl = np.dstack(np.where(ct == x))[0] - this_vxl = this_vxl + ch - voxels.append(this_vxl) - except BaseException: - raise - finally: - vol_prec.mip = old_mip - voxels = np.vstack(voxels) - - if not isinstance(bounds, type(None)): - voxels = voxels[np.all(voxels >= bounds[::2], axis=1)] - voxels = voxels[np.all(voxels < bounds[1::2], axis=1)] - - return voxels - - -def _chunks_to_nm(xyz_ch, vol, voxel_resolution=[4, 4, 40]): - """Map a chunk location to Euclidean space. - - Parameters - ---------- - xyz_ch : array-like - (N, 3) array of chunk indices. - vol : cloudvolume.CloudVolume - CloudVolume object associated with the chunked space. - voxel_resolution : list, optional - Voxel resolution. - - Returns - ------- - np.array - (N, 3) array of spatial points. - - """ - mip_scaling = vol.mip_resolution(0) // np.array(voxel_resolution, dtype=int) - - x_vox = np.atleast_2d(xyz_ch) * vol.mesh.meta.meta.graph_chunk_size - return ( - (x_vox + np.array(vol.mesh.meta.meta.voxel_offset(0))) - * voxel_resolution - * mip_scaling - ) + return cave_utils.get_voxels( + x, + mip=mip, + bounds=bounds, + datastack=datastack + ) \ No newline at end of file diff --git a/navis/interfaces/microns.py b/navis/interfaces/microns.py index 5f489ed3..8b180691 100644 --- a/navis/interfaces/microns.py +++ b/navis/interfaces/microns.py @@ -13,15 +13,10 @@ """Interface with MICrONS datasets: https://www.microns-explorer.org/.""" -from ..core import MeshNeuron, NeuronList -from .. import config, utils -import pandas as pd -import numpy as np -import warnings - -from concurrent.futures import ThreadPoolExecutor, as_completed +from .. import config from functools import lru_cache from textwrap import dedent +from . import cave_utils err_msg = dedent(""" Failed to import `caveclient` library. Please install using pip: @@ -40,11 +35,6 @@ except BaseException: raise - -logger = config.get_logger(__name__) -dataset = None - - @lru_cache(None) def _translate_datastack(datastack): """Translate datastack to source.""" @@ -89,6 +79,7 @@ def get_datastacks(microns_only=True): raise ImportError(err_msg) stacks = CAVEclient().info.get_datastacks() + if microns_only: stacks = [s for s in stacks if "minnie" in s or "pinky" in s] return stacks @@ -112,7 +103,7 @@ def get_cave_client(datastack="cortex65"): # Try mapping, else pass-through datastack = _translate_datastack(datastack) - return CAVEclient(datastack) + return cave_utils.get_cave_client(datastack) @lru_cache(None) @@ -133,332 +124,34 @@ def list_annotation_tables(datastack="cortex65"): List of available annotation tables. """ - return get_cave_client(datastack).materialize.get_tables() - - -@lru_cache(None) -def get_cloudvol(url, cache=True): - """Get (cached) CloudVolume for given segmentation. - - Parameters - ---------- - url : str - - """ - if not cv: - raise ImportError(err_msg) - - return cv.CloudVolume( - url, cache=cache, use_https=True, parallel=10, progress=False, fill_missing=True - ) - - -def get_somas(root_ids, table="nucleus_detection_v0", datastack="cortex65"): - """Fetch somas based on nuclei segmentation for given neuron(s). - - Since this is a nucleus detection you will find that some neurons do - not have an entry despite clearly having a soma. This is due to the - "avocado problem" where the nucleus is separate from the rest of the - soma. - - Important - --------- - This data currently only exists for the 'cortex65' datastack (i.e. - "minnie65_public_vXXX"). - - Parameters - ---------- - root_ids : int | list of ints | None - Root ID(s) for which to fetch soma infos. Use - `None` to fetch complete list of annotated nuclei. - table : str - Which table to use for nucleus annotations. Also - see the `microns.list_annotation_tables()` function. - datastack : "cortex65" | "cortex35" | "layer 2/3" | str - Which dataset to query. "cortex65", "cortex35" and "layer 2/3" - are internally mapped to the corresponding sources: for example, - "minnie65_public_vXXX" for "cortex65" where XXX is always the - most recent version). - - Returns - ------- - DataFrame - Pandas DataFrame with nuclei (see Examples). Root IDs - without a nucleus will simply not have an entry in the - table. - - """ - if datastack != "cortex65": - warnings.warn( - "To our knowledge there is no nucleus segmentation " - f'for "{datastack}". If that has changed please ' - "get in touch on navis' Github." - ) - - # Get/Initialize the CAVE client - client = get_cave_client(datastack) - - filter_in_dict = None - if not isinstance(root_ids, type(None)): - root_ids = utils.make_iterable(root_ids) - filter_in_dict = {"pt_root_id": root_ids} + return cave_utils.get_cave_client(datastack).materialize.get_tables() - return client.materialize.query_table(table, filter_in_dict=filter_in_dict) - -def fetch_neurons( - x, - *, - lod=2, - with_synapses=False, - datastack="cortex65", - parallel=True, - max_threads=4, - **kwargs, -): +def fetch_neurons(x, *, lod=2, + with_synapses=True, + datastack='h01_c3_flat', + parallel=True, + max_threads=4, + **kwargs): """Fetch neuron meshes. - - Notes - ----- - Synapses will be attached to the closest vertex on the mesh. - - Parameters - ---------- - x : str | int | list-like - Segment ID(s). Multiple Ids can be provided as list-like. - lod : int - Level of detail. Higher `lod` = coarser. This parameter - is ignored if the data source does not support multi-level - meshes. - with_synapses : bool, optional - If True will also attach synapses as `.connectors`. - datastack : "cortex65" | "cortex35" | "layer 2/3" | str - Which dataset to query. "cortex65", "cortex35" and "layer 2/3" - are internally mapped to the corresponding sources: for example, - "minnie65_public_vXXX" for "cortex65" where XXX is always the - most recent version). - parallel : bool - If True, will use parallel threads to fetch data. - max_threads : int - Max number of parallel threads to use. - **kwargs - Keyword arguments are passed through to the initialization - of the `navis.MeshNeurons`. - - Returns - ------- - navis.Neuronlist - Containing [`navis.MeshNeuron`][]. - """ - x = utils.make_iterable(x, force_type=int) - client = get_cave_client(datastack) - - vol = get_cloudvol(client.info.segmentation_source()) # this is cached - if datastack == "cortex65": - try: - somas = get_somas(x, datastack=datastack) - soma_pos = somas.set_index("pt_root_id").pt_position.to_dict() - except BaseException as e: - logger.warning("Failed to fetch somas via nucleus segmentation" f"(){e})") - soma_pos = {} - else: - soma_pos = {} - - nl = [] - with ThreadPoolExecutor(max_workers=1 if not parallel else max_threads) as executor: - futures = {} - for id in x: - f = executor.submit( - _fetch_single_neuron, - id, - vol=vol, - lod=lod, - client=client, - with_synapses=with_synapses, - source=datastack, - **kwargs, - ) - futures[f] = id - - with config.tqdm( - desc="Fetching", - total=len(x), - leave=config.pbar_leave, - disable=len(x) == 1 or config.pbar_hide, - ) as pbar: - for f in as_completed(futures): - id = futures[f] - pbar.update(1) - try: - nl.append(f.result()) - except Exception as exc: - print(f"{id} generated an exception:", exc) - - nl = NeuronList(nl) - - for n in nl: - if n.id in soma_pos: - n.soma_pos = np.array(soma_pos[n.id]) * [4, 4, 40] - else: - n.soma_pos = None - - return nl - - -def _fetch_single_neuron(id, lod, vol, client, with_synapses=False, **kwargs): - """Fetch a single neuron.""" - # Make sure we only use `lod` if that's actually supported by the source - if "MultiLevel" in str(type(vol.mesh)): - mesh = vol.mesh.get(id, lod=lod)[id] - else: - mesh = vol.mesh.get(id, deduplicate_chunk_boundaries=False)[id] - - n = MeshNeuron(mesh, id=id, units="nm", **kwargs) - - if with_synapses: - pre = client.materialize.synapse_query(pre_ids=id) - post = client.materialize.synapse_query(post_ids=id) - - syn_table = client.materialize.synapse_table - syn_info = client.materialize.get_table_metadata(syn_table) - vxl_size = np.array(syn_info["voxel_resolution"]).astype(int) - - to_concat = [] - if not pre.empty: - pre["type"] = "pre" - locs = np.vstack(pre["pre_pt_position"].values) - locs = locs * vxl_size - pre["x"], pre["y"], pre["z"] = locs[:, 0], locs[:, 1], locs[:, 2] - pre = pre[["id", "x", "y", "z", "type", "size"]].copy() - to_concat.append(pre) - if not post.empty: - post["type"] = "post" - locs = np.vstack(post["post_pt_position"].values) - locs = locs * vxl_size - post["x"], post["y"], post["z"] = locs[:, 0], locs[:, 1], locs[:, 2] - post = post[["id", "x", "y", "z", "type", "size"]].copy() - to_concat.append(post) - - if len(to_concat) == 1: - n.connectors = to_concat[0] - elif len(to_concat) == 2: - n.connectors = pd.concat(to_concat, axis=0).reset_index(drop=True) - - return n - - -def get_voxels(x, mip=0, bounds=None, datastack="cortex65"): - """Fetch voxels making a up given root ID. - - This is useful if you want to generate a high-resolution mesh for a neuron. - - Parameters - ---------- - x : int - A single root ID. - mip : int - Scale at which to fetch voxels. - bounds : list, optional - Bounding box [xmin, xmax, ymin, ymax, zmin, zmax] in voxel - space. For example, the voxel resolution for mip 0 - segmentation is 8 x 8 x 40 nm. - datastack : "cortex65" | "cortex35" | "layer 2/3" | str - Which dataset to query. "cortex65", "cortex35" and "layer 2/3" - are internally mapped to the corresponding sources: for example, - "minnie65_public_vXXX" for "cortex65" where XXX is always the - most recent version). - - Returns - ------- - voxels : (N, 3) np.ndarray - In voxel space according to `mip`. - - """ - client = get_cave_client(datastack) - # Need to get the graphene (not the precomputed) version of the data - vol_graphene = cv.CloudVolume( - client.chunkedgraph.cloudvolume_path, use_https=True, progress=False + return cave_utils.fetch_neurons( + x, + lod=lod, + with_synapses=with_synapses, + datastack=datastack, + parallel=parallel, + max_threads=max_threads, + **kwargs ) - url = client.info.get_datastack_info()["segmentation_source"] - vol_prec = get_cloudvol(url) # this is cached - - # Get L2 chunks making up this neuron - l2_ids = client.chunkedgraph.get_leaves(x, stop_layer=2) - - # Turn l2_ids into chunk indices - l2_ix = [ - np.array(vol_graphene.mesh.meta.meta.decode_chunk_position(l)) for l in l2_ids - ] - l2_ix = np.unique(l2_ix, axis=0) - - # Convert to nm - l2_nm = np.asarray(_chunks_to_nm(l2_ix, vol=vol_graphene)) - - # Convert back to voxel space (according to mip) - l2_vxl = l2_nm // vol_prec.meta.scales[mip]["resolution"] - - voxels = [] - ch_size = np.array(vol_graphene.mesh.meta.meta.graph_chunk_size) - ch_size = ch_size // (vol_prec.mip_resolution(mip) / vol_prec.mip_resolution(0)) - ch_size = np.asarray(ch_size).astype(int) - - if not isinstance(bounds, type(None)): - bounds = np.asarray(bounds) - if not bounds.ndim == 1 or len(bounds) != 6: - raise ValueError("`bounds` must be [xmin, xmax, ymin, ymax, zmin, zmax]") - l2_vxl = l2_vxl[np.all(l2_vxl >= (bounds[::2] - ch_size), axis=1)] - l2_vxl = l2_vxl[np.all(l2_vxl <= (bounds[1::2] + ch_size), axis=1)] - - old_mip = vol_prec.mip - try: - vol_prec.mip = mip - for ch in config.tqdm(l2_vxl, desc="Loading"): - ct = vol_prec[ - ch[0] : ch[0] + ch_size[0], - ch[1] : ch[1] + ch_size[1], - ch[2] : ch[2] + ch_size[2], - ][:, :, :, 0] - this_vxl = np.dstack(np.where(ct == x))[0] - this_vxl = this_vxl + ch - voxels.append(this_vxl) - except BaseException: - raise - finally: - vol_prec.mip = old_mip - voxels = np.vstack(voxels) - - if not isinstance(bounds, type(None)): - voxels = voxels[np.all(voxels >= bounds[::2], axis=1)] - voxels = voxels[np.all(voxels < bounds[1::2], axis=1)] - - return voxels - - -def _chunks_to_nm(xyz_ch, vol, voxel_resolution=[4, 4, 40]): - """Map a chunk location to Euclidean space. - - Parameters - ---------- - xyz_ch : array-like - (N, 3) array of chunk indices. - vol : cloudvolume.CloudVolume - CloudVolume object associated with the chunked space. - voxel_resolution : list, optional - Voxel resolution. - - Returns - ------- - np.array - (N, 3) array of spatial points. +def get_voxels(x, mip=0, bounds=None, datastack='h01_c3_flat'): + """Fetch voxels making a up given root ID. """ - mip_scaling = vol.mip_resolution(0) // np.array(voxel_resolution, dtype=int) - - x_vox = np.atleast_2d(xyz_ch) * vol.mesh.meta.meta.graph_chunk_size - return ( - (x_vox + np.array(vol.mesh.meta.meta.voxel_offset(0))) - * voxel_resolution - * mip_scaling - ) + return cave_utils.get_voxels( + x, + mip=mip, + bounds=bounds, + datastack=datastack + ) \ No newline at end of file From f9d8c578719755b1398c052f9162f842dce50f70 Mon Sep 17 00:00:00 2001 From: jinhan Date: Fri, 11 Oct 2024 17:22:33 -0400 Subject: [PATCH 3/3] resolve comments --- docs/examples/4_remote/plot_04_remote_h01.py | 15 +++-- navis/interfaces/cave_utils.py | 22 +++---- navis/interfaces/h01.py | 62 +++++++++++++++++--- navis/interfaces/microns.py | 57 +++++++++++++++++- 4 files changed, 130 insertions(+), 26 deletions(-) diff --git a/docs/examples/4_remote/plot_04_remote_h01.py b/docs/examples/4_remote/plot_04_remote_h01.py index 9f8c15f1..f21c81b7 100644 --- a/docs/examples/4_remote/plot_04_remote_h01.py +++ b/docs/examples/4_remote/plot_04_remote_h01.py @@ -2,7 +2,10 @@ H01 Dataset ============ -In this notebook, you can learn how to work with the [H01 dataset](https://www.science.org/doi/10.1126/science.adk4858) using Navis. The H01 dataset contains 57,000 cells and 150 million synapses from a cubic millimeter of the human temporal cortex, which is [proofreading](https://h01-release.storage.googleapis.com/proofreading.html) in the CAVE. With this interface, you can access both a snapshot of the proofread dataset and the latest dataset using CAVEclient. +In this notebook, you can learn how to work with the H01 dataset using Navis. + +The [H01 dataset](https://www.science.org/doi/10.1126/science.adk4858) contains 57,000 cells and 150 million synapses from a cubic millimeter of the human temporal cortex, which is [proofreading](https://h01-release.storage.googleapis.com/proofreading.html) in the CAVE. +With this interface, you can access both a snapshot of the proofread dataset and the latest dataset using CAVEclient. """ @@ -16,19 +19,18 @@ # %% client = h01.get_cave_client() - # %% ## Query Tables # %% client.materialize.get_versions() -print("Mat version: ", client.materialize.version) + +# %% client.materialize.get_tables() # %% ### Query Materialized Synapse Table - # %% client.materialize.synapse_query(limit=10) @@ -38,12 +40,13 @@ # pre_ids=[ADD YOUR ROOT ID], ) syn.head() + +# %% print(len(syn)) # %% ### Live Synapse Queries - # %% import datetime # check if root ID is the most recent root ID @@ -92,7 +95,7 @@ interneurons_ds = navis.simplify_mesh(interneurons, F=1/3) # %% -interneuron_ds +interneurons_ds # %% import seaborn as sns diff --git a/navis/interfaces/cave_utils.py b/navis/interfaces/cave_utils.py index 9a4534d9..4d3ae7ec 100644 --- a/navis/interfaces/cave_utils.py +++ b/navis/interfaces/cave_utils.py @@ -35,7 +35,7 @@ def get_cave_client(datastack="cortex65", server_address=None): Parameters ---------- - datastack : "cortex65" | "cortex35" | "layer 2/3" | str + datastack : "cortex65" | "cortex35" | "layer 2/3" | "h01_c3_flat" | str Which dataset to query. "cortex65", "cortex35" and "layer 2/3" are internally mapped to the corresponding sources: for example, "minnie65_public_vXXX" for "cortex65" where XXX is always the @@ -84,7 +84,7 @@ def _get_somas(root_ids, table="nucleus_detection_v0", datastack="cortex65"): table : str Which table to use for nucleus annotations. Also see the `microns.list_annotation_tables()` function. - datastack : "cortex65" | "cortex35" | "layer 2/3" | str + datastack : "cortex65" | "cortex35" | "layer 2/3" | "h01_c3_flat" | str Which dataset to query. "cortex65", "cortex35" and "layer 2/3" are internally mapped to the corresponding sources: for example, "minnie65_public_vXXX" for "cortex65" where XXX is always the @@ -137,10 +137,11 @@ def fetch_neurons(x, lod, meshes. with_synapses : bool, optional If True will also attach synapses as ``.connectors``. - datastack : "cortex65" | "cortex35" | "layer 2/3" - Which dataset to use. Internally these are mapped to the - corresponding sources (e.g. "minnie65_public_v117" for - "cortex65"). + datastack : "cortex65" | "cortex35" | "layer 2/3" | "h01_c3_flat" | str + Which dataset to query. "cortex65", "cortex35" and "layer 2/3" + are internally mapped to the corresponding sources: for example, + "minnie65_public_vXXX" for "cortex65" where XXX is always the + most recent version). parallel : bool If True, will use parallel threads to fetch data. max_threads : int @@ -262,10 +263,11 @@ def get_voxels(x, mip, bounds, datastack): Bounding box [xmin, xmax, ymin, ymax, zmin, zmax] in voxel space. For example, the voxel resolution for mip 0 segmentation is 8 x 8 x 40 nm. - datastack : "h01_c3_flat" - Which dataset to use. Internally these are mapped to the - corresponding sources (e.g. "h01_c3_flat" for - "h01_c3_flat"). + datastack : "cortex65" | "cortex35" | "layer 2/3" | "h01_c3_flat" | str + Which dataset to query. "cortex65", "cortex35" and "layer 2/3" + are internally mapped to the corresponding sources: for example, + "minnie65_public_vXXX" for "cortex65" where XXX is always the + most recent version). Returns ------- diff --git a/navis/interfaces/h01.py b/navis/interfaces/h01.py index 06715d59..92534b80 100644 --- a/navis/interfaces/h01.py +++ b/navis/interfaces/h01.py @@ -1,24 +1,50 @@ -from functools import lru_cache from . import cave_utils -@lru_cache(None) +DATASTACK = "h01_c3_flat" +SERVER_ADDRESS = "https://global.brain-wire-test.org/" + def get_cave_client(): """Get caveclient for H01 dataset. """ - datastack = "h01_c3_flat" - server_address = "https://global.brain-wire-test.org/" - - return cave_utils.get_cave_client(datastack, server_address) + return cave_utils.get_cave_client(DATASTACK, SERVER_ADDRESS) def fetch_neurons(x, *, lod=2, with_synapses=True, - datastack='h01_c3_flat', + datastack=DATASTACK, parallel=True, max_threads=4, **kwargs): """Fetch neuron meshes. - """ + Notes + ----- + Synapses will be attached to the closest vertex on the mesh. + + Parameters + ---------- + x : str | int | list-like + Segment ID(s). Multiple Ids can be provided as list-like. + lod : int + Level of detail. Higher ``lod`` = coarser. This parameter + is ignored if the data source does not support multi-level + meshes. + with_synapses : bool, optional + If True will also attach synapses as ``.connectors``. + datastack : str + DATASTACK. + parallel : bool + If True, will use parallel threads to fetch data. + max_threads : int + Max number of parallel threads to use. + **kwargs + Keyword arguments are passed through to the initialization + of the ``navis.MeshNeurons``. + + Returns + ------- + navis.Neuronlist + Containing :class:`navis.MeshNeuron`. + """ return cave_utils.fetch_neurons( x, lod=lod, @@ -29,8 +55,26 @@ def fetch_neurons(x, *, lod=2, **kwargs ) -def get_voxels(x, mip=0, bounds=None, datastack='h01_c3_flat'): +def get_voxels(x, mip=0, bounds=None, datastack=DATASTACK): """Fetch voxels making a up given root ID. + + Parameters + ---------- + x : int + A single root ID. + mip : int + Scale at which to fetch voxels. + bounds : list, optional + Bounding box [xmin, xmax, ymin, ymax, zmin, zmax] in voxel + space. For example, the voxel resolution for mip 0 + segmentation is 8 x 8 x 40 nm. + datastack : str + DATASTACK. + + Returns + ------- + voxels : (N, 3) np.ndarray + In voxel space according to `mip`. """ return cave_utils.get_voxels( x, diff --git a/navis/interfaces/microns.py b/navis/interfaces/microns.py index 8b180691..4488e54f 100644 --- a/navis/interfaces/microns.py +++ b/navis/interfaces/microns.py @@ -85,7 +85,6 @@ def get_datastacks(microns_only=True): return stacks -@lru_cache(None) def get_cave_client(datastack="cortex65"): """Get caveclient for given datastack. @@ -134,6 +133,39 @@ def fetch_neurons(x, *, lod=2, max_threads=4, **kwargs): """Fetch neuron meshes. + + Notes + ----- + Synapses will be attached to the closest vertex on the mesh. + + Parameters + ---------- + x : str | int | list-like + Segment ID(s). Multiple Ids can be provided as list-like. + lod : int + Level of detail. Higher ``lod`` = coarser. This parameter + is ignored if the data source does not support multi-level + meshes. + with_synapses : bool, optional + If True will also attach synapses as ``.connectors``. + datastack : "cortex65" | "cortex35" | "layer 2/3" | str + Which dataset to query. "cortex65", "cortex35" and "layer 2/3" + are internally mapped to the corresponding sources: for example, + "minnie65_public_vXXX" for "cortex65" where XXX is always the + most recent version). + parallel : bool + If True, will use parallel threads to fetch data. + max_threads : int + Max number of parallel threads to use. + **kwargs + Keyword arguments are passed through to the initialization + of the ``navis.MeshNeurons``. + + Returns + ------- + navis.Neuronlist + Containing :class:`navis.MeshNeuron`. + """ return cave_utils.fetch_neurons( @@ -148,6 +180,29 @@ def fetch_neurons(x, *, lod=2, def get_voxels(x, mip=0, bounds=None, datastack='h01_c3_flat'): """Fetch voxels making a up given root ID. + + + Parameters + ---------- + x : int + A single root ID. + mip : int + Scale at which to fetch voxels. + bounds : list, optional + Bounding box [xmin, xmax, ymin, ymax, zmin, zmax] in voxel + space. For example, the voxel resolution for mip 0 + segmentation is 8 x 8 x 40 nm. + datastack : "cortex65" | "cortex35" | "layer 2/3" | str + Which dataset to query. "cortex65", "cortex35" and "layer 2/3" + are internally mapped to the corresponding sources: for example, + "minnie65_public_vXXX" for "cortex65" where XXX is always the + most recent version). + + Returns + ------- + voxels : (N, 3) np.ndarray + In voxel space according to `mip`. + """ return cave_utils.get_voxels( x,