|
| 1 | +import numpy as np |
| 2 | +from brainbox.numerical import ismember |
| 3 | +from oneibl.one import ONE |
| 4 | +from ibllib.pipes import histology |
| 5 | +from ibllib.atlas import AllenAtlas |
| 6 | +from ibllib.ephys.neuropixel import TIP_SIZE_UM, SITES_COORDINATES |
| 7 | +from ibllib.pipes.ephys_alignment import EphysAlignment |
| 8 | +import time |
| 9 | + |
| 10 | +PROV_2_VAL = { |
| 11 | + 'Resolved': 90, |
| 12 | + 'Ephys aligned histology track': 70, |
| 13 | + 'Histology track': 50, |
| 14 | + 'Micro-manipulator': 30, |
| 15 | + 'Planned': 10} |
| 16 | + |
| 17 | +VAL_2_PROV = {v: k for k, v in PROV_2_VAL.items()} |
| 18 | + |
| 19 | + |
| 20 | +class ProbeModel: |
| 21 | + def __init__(self, one=None, ba=None): |
| 22 | + |
| 23 | + self.one = one or ONE() |
| 24 | + self.ba = ba or AllenAtlas(25) |
| 25 | + self.traj = {'Planned': {}, |
| 26 | + 'Micro-manipulator': {}, |
| 27 | + 'Histology track': {}, |
| 28 | + 'Ephys aligned histology track': {}, |
| 29 | + 'Resolved': {}, |
| 30 | + 'Best': {}} |
| 31 | + |
| 32 | + self.get_traj_for_provenance(provenance='Histology track', django=['x__isnull,False']) |
| 33 | + self.get_traj_for_provenance(provenance='Ephys aligned histology track') |
| 34 | + self.get_traj_for_provenance(provenance='Ephys aligned histology track', |
| 35 | + django=['probe_insertion__json__extended_qc__' |
| 36 | + 'alignment_resolved,True'], prov_dict='Resolved') |
| 37 | + self.find_traj_is_best(provenance='Histology track') |
| 38 | + self.find_traj_is_best(provenance='Ephys aligned histology track') |
| 39 | + self.traj['Resolved']['is_best'] = np.arange(len(self.traj['Resolved']['traj'])) |
| 40 | + |
| 41 | + @staticmethod |
| 42 | + def get_traj_info(traj): |
| 43 | + return traj['probe_insertion'], traj['x'], traj['y'] |
| 44 | + |
| 45 | + def get_traj_for_provenance(self, provenance='Histology track', django=None, |
| 46 | + prov_dict=None): |
| 47 | + |
| 48 | + if django is None: |
| 49 | + django = [] |
| 50 | + |
| 51 | + if prov_dict is None: |
| 52 | + prov_dict = provenance |
| 53 | + |
| 54 | + django_base = ['probe_insertion__session__project__name__icontains,' |
| 55 | + 'ibl_neuropixel_brainwide_01'] |
| 56 | + django_str = ','.join(django_base + django) |
| 57 | + |
| 58 | + self.traj[prov_dict]['traj'] = np.array(self.one.alyx.rest('trajectories', 'list', |
| 59 | + provenance=provenance, |
| 60 | + django=django_str)) |
| 61 | + ins_ids, x, y = zip(*[self.get_traj_info(traj) for traj in self.traj[prov_dict]['traj']]) |
| 62 | + self.traj[prov_dict]['ins'] = np.array(ins_ids) |
| 63 | + self.traj[prov_dict]['x'] = np.array(x) |
| 64 | + self.traj[prov_dict]['y'] = np.array(y) |
| 65 | + |
| 66 | + def compute_best_for_provenance(self, provenance='Histology track'): |
| 67 | + val = PROV_2_VAL[provenance] |
| 68 | + prov_to_include = [] |
| 69 | + for k, v in VAL_2_PROV.items(): |
| 70 | + if k >= val: |
| 71 | + prov_to_include.append(v) |
| 72 | + |
| 73 | + for iP, prov in enumerate(prov_to_include): |
| 74 | + if not 'is_best' in self.traj[prov].keys(): |
| 75 | + self.find_traj_is_best(prov) |
| 76 | + |
| 77 | + if iP == 0: |
| 78 | + self.traj['Best']['traj'] = self.traj[prov]['traj'][self.traj[prov]['is_best']] |
| 79 | + self.traj['Best']['ins'] = self.traj[prov]['ins'][self.traj[prov]['is_best']] |
| 80 | + self.traj['Best']['x'] = self.traj[prov]['x'][self.traj[prov]['is_best']] |
| 81 | + self.traj['Best']['y'] = self.traj[prov]['y'][self.traj[prov]['is_best']] |
| 82 | + else: |
| 83 | + self.traj['Best']['traj'] = np.r_[self.traj['Best']['traj'], |
| 84 | + (self.traj[prov]['traj'] |
| 85 | + [self.traj[prov]['is_best']])] |
| 86 | + self.traj['Best']['ins'] = np.r_[self.traj['Best']['ins'], |
| 87 | + (self.traj[prov]['ins'] |
| 88 | + [self.traj[prov]['is_best']])] |
| 89 | + self.traj['Best']['x'] = np.r_[self.traj['Best']['x'], |
| 90 | + self.traj[prov]['x'][self.traj[prov]['is_best']]] |
| 91 | + self.traj['Best']['y'] = np.r_[self.traj['Best']['y'], |
| 92 | + self.traj[prov]['y'][self.traj[prov]['is_best']]] |
| 93 | + |
| 94 | + def find_traj_is_best(self, provenance='Histology track'): |
| 95 | + val = PROV_2_VAL[provenance] |
| 96 | + next_provenance = VAL_2_PROV[val + 20] |
| 97 | + |
| 98 | + if not 'traj' in self.traj[provenance].keys(): |
| 99 | + self.get_traj_for_provenance(provenance) |
| 100 | + if not 'traj' in self.traj[next_provenance].keys(): |
| 101 | + self.get_traj_for_provenance(next_provenance) |
| 102 | + |
| 103 | + isin, _ = ismember(self.traj[provenance]['ins'], |
| 104 | + self.traj[next_provenance]['ins']) |
| 105 | + self.traj[provenance]['is_best'] = np.where(np.invert(isin))[0] |
| 106 | + |
| 107 | + # Special exception for planned provenance |
| 108 | + if provenance == 'Planned': |
| 109 | + next_provenance = VAL_2_PROV[val + 40] |
| 110 | + if not 'traj' in self.traj[next_provenance].keys(): |
| 111 | + self.get_traj_for_provenance(next_provenance) |
| 112 | + isin, _ = ismember(self.traj[provenance]['ins'][self.traj[provenance]['is_best']], |
| 113 | + self.traj[next_provenance]['ins']) |
| 114 | + self.traj[provenance]['is_best'] = (self.traj[provenance]['is_best'] |
| 115 | + [np.where(np.invert(isin))[0]]) |
| 116 | + |
| 117 | + def get_channels(self, provenance): |
| 118 | + |
| 119 | + depths = SITES_COORDINATES[:, 1] |
| 120 | + start = time.time() |
| 121 | + |
| 122 | + # Need to account for case when no insertions for planned and micro can skip the insertions |
| 123 | + insertions = [] |
| 124 | + step = 150 |
| 125 | + for i in range(np.ceil(self.traj[provenance]['ins'].shape[0] / step).astype(np.int)): |
| 126 | + insertions += self.one.alyx.rest( |
| 127 | + 'insertions', 'list', django=f"id__in,{list(self.traj[provenance]['ins'][i * step:(i + 1) * step])}") |
| 128 | + |
| 129 | + ins_id = np.ravel([np.where(ins['id'] == self.traj[provenance]['ins']) |
| 130 | + for ins in insertions]) |
| 131 | + |
| 132 | + end = time.time() |
| 133 | + print(end-start) |
| 134 | + |
| 135 | + start1 = time.time() |
| 136 | + for iT, (traj, ins) in enumerate(zip(self.traj[provenance]['traj'][ins_id], insertions)): |
| 137 | + try: |
| 138 | + xyz_picks = np.array(ins['json']['xyz_picks']) / 1e6 |
| 139 | + if traj['provenance'] == 'Histology track': |
| 140 | + xyz_picks = xyz_picks[np.argsort(xyz_picks[:, 2]), :] |
| 141 | + xyz_channels = histology.interpolate_along_track(xyz_picks, (depths + TIP_SIZE_UM) / 1e6) |
| 142 | + else: |
| 143 | + align_key = ins['json']['extended_qc']['alignment_stored'] |
| 144 | + feature = traj['json'][align_key][0] |
| 145 | + track = traj['json'][align_key][1] |
| 146 | + ephysalign = EphysAlignment(xyz_picks, depths, track_prev=track, |
| 147 | + feature_prev=feature, |
| 148 | + brain_atlas=self.ba, speedy=True) |
| 149 | + xyz_channels = ephysalign.get_channel_locations(feature, track) |
| 150 | + |
| 151 | + if iT == 0: |
| 152 | + all_channels = xyz_channels |
| 153 | + else: |
| 154 | + all_channels = np.r_[all_channels, xyz_channels] |
| 155 | + except Exception as err: |
| 156 | + print(err) |
| 157 | + print(traj['id']) |
| 158 | + |
| 159 | + end = time.time() |
| 160 | + print(end-start1) |
| 161 | + |
| 162 | + return all_channels |
| 163 | + |
| 164 | +from scipy.signal import fftconvolve |
| 165 | +cvol = np.zeros(ba.image.shape, dtype=np.float) |
| 166 | +val, counts = np.unique(ba._lookup(all_channels), return_counts=True) |
| 167 | +cvol[np.unravel_index(val, cvol.shape)] = counts |
| 168 | +# |
| 169 | +from ibllib.dsp import fcn_cosine |
| 170 | +DIST_FCN = np.array([100, 150]) / 1e6 |
| 171 | +dx = ba.bc.dx |
| 172 | +template = np.arange(- np.max(DIST_FCN) - dx, np.max(DIST_FCN) + 2 * dx, dx) ** 2 |
| 173 | +kernel = sum(np.meshgrid(template, template, template)) |
| 174 | +kernel = 1 - fcn_cosine(DIST_FCN)(np.sqrt(kernel)) |
| 175 | +# |
| 176 | +cvol = fftconvolve(cvol, kernel) |
| 177 | +# |
| 178 | +# |
| 179 | +# cvol[np.unravel_index(ba._lookup(all_channels), cvol.shape)] = 1 |
| 180 | + |
| 181 | +# from ibllib.atlas import AllenAtlas |
| 182 | +# ba = AllenAtlas() |
| 183 | +# import vedo |
| 184 | +# import numpy as np |
| 185 | +# |
| 186 | +# actor = vedo.Volume(ba.image, c='bone', spacing = np.array([25]*3), mapper='smart', mode=0, alphaGradient=0.5) |
| 187 | +# plt = vedo.Plotter() |
| 188 | +# plt.add(actor) |
| 189 | +# plt.show() |
| 190 | + |
| 191 | + |
0 commit comments