Find closest trajectories

Finds the probe insertions that are closest a a specified probe insertion of interest and plots their location in a 3D image. Additionally prints out the information of closeby sessions

# Author: Mayo Faulkner
# import modules
import numpy as np
from one.api import ONE
from mayavi import mlab

import ibllib.pipes.histology as histology
import ibllib.atlas as atlas
from atlaselectrophysiology import rendering

# Instantiate brain atlas and one
brain_atlas = atlas.AllenAtlas(25)
one = ONE(base_url='https://openalyx.internationalbrainlab.org')

# Find all trajectories with histology tracing
all_hist = one.alyx.rest('trajectories', 'list', provenance='Ephys aligned histology track')
# Some do not have tracing, exclude these ones
sess_with_hist = [sess for sess in all_hist if sess['x'] is not None]
traj_ids = [sess['id'] for sess in sess_with_hist]
# Compute trajectory objects for each of the trajectories
trajectories = [atlas.Insertion.from_dict(sess) for sess in sess_with_hist]

# Find the trajectory of the id that you want to find closeby probe insertions for
subject = 'CSH_ZAD_029'
date = '2020-09-19'
probe_label = 'probe01'
traj_origin_id = one.alyx.rest('trajectories', 'list', provenance='Ephys aligned histology track',
                               subject=subject, date=date, probe=probe_label)[0]['id']
# Find the index of this trajectory in the list of all trajectories
chosen_traj = traj_ids.index(traj_origin_id)

# Define active part of probe ~ 200um from tip and ~ (200 + 3900)um to top of channels
depths = np.arange(200, 4100, 20) / 1e6
traj_coords = np.empty((len(traj_ids), len(depths), 3))

# For each trajectory compute the xyz coords at positions depths along trajectory
for iT, traj in enumerate(trajectories):
    traj_coords[iT, :] = histology.interpolate_along_track(np.vstack([traj.tip, traj.entry]),
                                                           depths)

# Find the average distance between all positions compared to trjaectory of interest
avg_dist = np.mean(np.sqrt(np.sum((traj_coords - traj_coords[chosen_traj]) ** 2, axis=2)), axis=1)

# Sort according to those that are closest
closest_traj = np.argsort(avg_dist)

close_sessions = []
# Make a 3D plot showing trajectory of interest (in black) and the 10 nearest trajectories (blue)
fig = rendering.figure(grid=False)
for iSess, sess_idx in enumerate(closest_traj[0:10]):

    mlapdv = brain_atlas.xyz2ccf(traj_coords[sess_idx])
    if iSess == 0:
        mlab.plot3d(mlapdv[:, 1], mlapdv[:, 2], mlapdv[:, 0],
                    line_width=1, tube_radius=10, color=(0, 0, 0))
    else:
        mlab.plot3d(mlapdv[:, 1], mlapdv[:, 2], mlapdv[:, 0],
                    line_width=1, tube_radius=10, color=(0.0, 0.4, 0.5))

    mlab.text3d(mlapdv[0, 1], mlapdv[0, 2], mlapdv[0, 0], str(iSess),
                line_width=4, color=(0, 0, 0), figure=fig, scale=150)

    close_sessions.append((sess_with_hist[sess_idx]['session']['subject'] + ' ' +
                           sess_with_hist[sess_idx]['session']['start_time'][:10] +
                           ' ' + sess_with_hist[sess_idx]['probe_name'] + ': dist = ' +
                           str(avg_dist[closest_traj[iSess]] * 1e6)))

print(close_sessions)

['CSH_ZAD_029 2020-09-19 probe01: dist = 0.0', 'SWC_043 2020-09-21 probe00: dist = 2456.5049323924577', 'KS023 2019-12-10 probe01: dist = 3161.001547693022', 'CSHL049 2020-01-08 probe00: dist = 3214.5209187462738']